K Hoffmann

J Bel

E Gaztañaga

Linear and non-linear bias: predictions vs. measurements

Keywords: methods, analytical -methods, statistical -galaxies, haloes -dark matter -large-scale structure of Universe

We study the linear and non-linear bias parameters which determine the mapping between the distributions of galaxies and the full matter density fields, comparing different measurements and predictions. Associating galaxies with dark matter haloes in the MICE Grand Challenge N-body simulation we directly measure the bias parameters by comparing the smoothed density fluctuations of haloes and matter in the same region at different positions as a function of smoothing scale. Alternatively we measure the bias parameters by matching the probability distributions of halo and matter density fluctuations, which can be applied to observations. These direct bias measurements are compared to corresponding measurements from two-point and different third-order correlations, as well as predictions from the peak-background model, which we presented in previous articles using the same data. We find an overall variation of the linear bias measurements and predictions of ∼ 5% with respect to results from two-point correlations for different halo samples with masses between ∼ 10 12 -10 15 h -1 M ⊙ at the redshifts z = 0.0 and 0.5. Variations between the second-and thirdorder bias parameters from the different methods show larger variations, but with consistent trends in mass and redshift. The various bias measurements reveal a tight relation between the linear and the quadratic bias parameters, which is consistent with results from the literature based on simulations with different cosmologies. Such a universal relation might improve constraints on cosmological models, derived from second-order clustering statistics at small scales or higher-order clustering statistics.

INTRODUCTION

The increase of data from upcoming and next generation of galaxy surveys is pulling down errors on the observed statistics of the large-scale galaxy distribution. Thus, the inferences on cosmological models from these statistics requires a modelling of cosmological fluids and their statistical properties with an accuracy of at least the same order of magnitude as the observational errors.

One of the largest uncertainties comes from the modelling of the mapping between the observed fluctuations of galaxies to the fluctuation of the underlying matter distribution (hereafter referred to as δg and δm respectively). These fluctuations are defined as normalised deviations of the density ρ, smoothed typically with a top-hat window of characteristic scale R, from the mean density of the universe ρ at the position r,

δ(r) ≡ ρ(r) - ρ ρ
.

(1)

The mapping from δg to δm is described by the so-called bias function, δg = F [δm, ∇ij Φv], where ∇ijΦv are secondorder derivatives of the velocity potential. The latter relate δg to the matter distribution beyond the smoothing scale R and are therefore referred as non-local contributions to the bias model (Chan et al. 2012;Baldauf et al. 2012). For sufficiently large smoothing scales such non-local contributions may be negligible, which allows for a local description of biasing. Thus, the bias function can be modelled as a Taylor expansion in terms of the matter fluctuations (Fry and Gaztanaga 1993)

δ h = F [δm] ≃ N i=0 bi i! δ i m , (2)
where bi are the bias parameters. However, for smaller smoothing scales, which are typically used for studying the two-and three-point statistics of the galaxy distribution, non-local contributions need to be considered. When fluctuations in the density and velocity potential of the matter distribution are sufficiently small the non-local bias model may be described by its second-order expansion, which is hereafter referred to as non-local quadratic bias model

δg = b1 δm + 1 2 [c2(δ 2 m -δ 2 m ) + g2G2] , (3)
where b1 and c2 ≡ b2/b1 are referred to as linear and quadratic (or second-order) bias parameters respectively. The non-local contribution to the bias function consists in that case of the product of the second-order non-local bias parameter g2 and the smoothed second-order Gallileon G2(r) = β12θv(k1)θv(k2) Ŵ [k12R]e ik 12 •r d 3 k1d 3 k2, (4)

where ki and k12 ≡ k2 -k1 are wave vectors of density oscillations, β12 ≡ 1 -( k1 • k2) 2 represents the modecoupling between density oscillations which describes tidal forces, θv ≡ ∇ 2 Φv is the divergence of the normalised velocity field (v/H/f ) and Ŵ [k12R] is the window function in Fourier space. Note that in the case g2 = 0 equation (3) corresponds to the local quadratic bias model.

The bias parameters are highly relevant for constraining cosmological models via the growth of matter fluctuations, derived from second-order statistics of the observed galaxy distribution. In particular at large scales, the linear bias factor b1 is completely degenerate with the linear growth factor. Hence, growth-independent measurements of b1 can strongly tighten cosmological constraints from galaxy surveys.

Third-order statistics probes the linear and quadratic bias parameters independently of the growth and can be used to break the growth-bias degeneracy in the secondorder statistics. Furthermore, the second-order bias measurements from third-order statistics allow growth measurements from second-order statistics at small scales, where non-linear and non-local terms contribute significantly to the signal. Such a small scale analysis would strongly improve the cosmological constraining power of galaxy surveys.

However, the value of combining second-with thirdorder statistics for constraining cosmological models with high precision arises from a detailed understanding of how exactly the bias parameters enter these statistics at different scales, redshifts and for different samples of tracers.

We therefore investigated in previous works (Hoffmann et al. 2015b;Bel et al. 2015, where the latter is hereafter referred to as paper I) how accurately the linear bias can be measured from different third-order statistics. These studies were based on the large cosmological MICE Grand Challenge (hereafter referred to as MICE-GC) N-body simulation in which haloes were detected as tracers of the cosmic web and associated with galaxies. The fact that the dark matter distribution is accessible in simulations allows for reliable measurements of the linear bias via second-order statistics, which can then be used as a reference for validating linear bias measurements from third-order statistics as well as theory predictions. Note that a reliable reference for validating higher order bias is currently only provided by running separate universe N-body simulations (Wagner et al. 2015;Lazeyras et al. 2016), which is beyond the scope of this article.

Alternatively to growth-independent bias measurements from third-order statistics one can employ bias predictions from the peak-background split model for breaking the growth-bias degeneracy. In Hoffmann et al. (2015a, hereafter referred to paper II) we tested linear bias predictions in the MICE-GC simulations, confirming reports on inaccuracies of these predictions in the literature on a larger mass range thanks to the large volume and resolution of the MICE-GC simulation. As in case of the third-order statistics no reliable measurements for validating the predicted non-linear bias parameters were available. However, an interesting outcome of this analysis was a simple universal relation between the linear and non-linear bias parameters in the peak-background split model, which is independent of redshift and cosmology for halo samples with b1 2.

In the present work we conclude our series of articles, presenting results from an other method for allowing to measure bias parameters in the MICE-GC simulation. This method is based on a direct comparison of the halo and matter density fluctuations and may therefore be seen as the most direct way of measuring bias parameters. We also study a variant of this method, which is based on abundance matching of tracers and matter density contrasts. The interest of this alternative method is that it can be applied in galaxy surveys to estimate to bias function. Both methods deliver linear, as well as non-linear bias measurements. We compare these new measurements with the most reliable results from our previous work (paper I and II, see Table 1), which allows for the validation of how well bias can be measured with each approach.

The different measurements of linear and non-linear bias parameters furthermore allow for a validation of the universal relation between the bias parameters, which we found for peak-background split predictions in paper II. The strength of our comparison emerges from the large halo mass range of roughly 10 12 -10 15 h -1 M⊙, probed by the MICE-GC simulation as well as the fact that we use the same halo mass samples throughout the whole comparison project.

The remainder of this article is organised as follows. In Section 2 we describe the MICE-GC simulation as well as the halo samples on which our analysis is based on. The different methods for obtaining the bias parameters are briefly reviewed in Section 3, while details on bias measurements from the comparison of density contrasts are given in the appendix. We present our results in Section 4 which we summarise and discuss in Section 5.

SIMULATION AND HALO SAMPLES

Our analysis is based on the Grand Challenge run of the Marenostrum Institut de Ciències de l'Espai (MICE) simula- tion suite to which we refer to as MICE-GC in the following.

Starting from small initial density fluctuations at redshift z = 100 the formation of large scale cosmic structure was computed with 4096 3 gravitationally interacting collisionless dark matter particles in a 3072 h -1 Mpc box using the GADGET -2 code (Springel 2005) with a softening length of 50 h -1 kpc. The initial conditions were generated using the Zel'dovich approximation and a CAMB power spectrum with the power law index of ns = 0.95, which was normalised to be σ8 = 0.8 at z = 0. The cosmic expansion is described by the ΛCDM model for a flat universe with a mass density of Ωm = Ω dm + Ω b = 0.25. The density of the baryonic mass is set to Ω b = 0.044 and Ω dm is the dark matter density.

The dimensionless Hubble parameter is set to h = 0.7. More details and validation test on this simulation can be found in Fosalba et al. (2013). Dark matter haloes were identified as Friends-of-Friends groups (Davis et al. 1985) with a linking length of 0.2 in units of the mean particle separation. These halo catalogs and the corresponding validation checks are presented in Crocce et al. (2013). In the present analysis we use the dark matter field as well as haloes in the comoving outputs at redshift z = 0.0 and 0.5. The haloes are thereby divided into the four redshift independent mass samples M0 -M3 (defined in Table 2), which span a mass range from Milky Way like haloes up to haloes of massive galaxy clusters.

BIAS ESTIMATORS

In this section we summarise the different methods for measuring and predicting linear and non-linear bias parameters. An overview of the methods with the abbreviations and references for the corresponding bias is given in Table 3.

δ h -δm relation

In simulations we can measure the density contrasts of haloes and matter (hereafter referred to as δ h and δm respectively), which allows for a direct estimation of the bias parameters bi in equation (2). In this sub-section, we describe the two methods which we employ for these measurements.

The first method has been widely used in the literature and may be seen as the most direct way of measuring bias, provided that the biasing function is deterministic and local. It consists in fitting a polynomial of order n to the δ h (r i )δm(r i ) relation, which can be measured from the smoothed density fields of haloes and matter at different positions r i in a simulation. The scatter between haloes and matter density contrasts in the MICE-GC simulation is shown in Fig. 1 at redshift z = 0.0 for a smoothing scale of R = 60h -1 Mpc.

One can see that the average relation between δm and δ h is different regarding the halo sample: the slope of the most massive one (M3) is higher than the less massive (M0). This is expected because the linear bias which controls the slope of the δ h -δm relation is higher for the sample M3. In general, the best-fit parameters of the polynomial of order n fitted to the scatter plot are then identified as the bias coefficients. Note that the parameter b0 is in principle constrained by requiring that by definition the average of the halo density contrast is null, δ h = 0. In case of a cubic biasing relation it leads to the formal expression

b0 = -b1σ 2 dm c2 2 + c3 6 S 3,dm σ 2 dm , (5)
where S 3,dm and σ 2 dm are respectively the skewness and the variance of the dark matter field smoothed on scale R. Note that c3 is the cubic bias coefficient defined as c3 ≡ b3/b1. We will refer to this method as δ h -δm in the remainder of this article.

In practice, the choice of the order of the polynomial could affect the estimation of the bias coefficients b1, c2 and c3 (see Lazeyras et al. 2016). One can expect that this truncation introduces a dependence of the bias parameters on the smoothing scale, used to estimate the density contrasts. We therefore perform the same convergence test proposed by Lazeyras et al. (2016), which consists in comparing the bias coefficients obtained with a third-, fourth-and fifthorder polynomial. In the case of third-order polynomials we test, in addition, the impact of the choice of b0 on the estimated bias parameters. To do so we let b0 being free and compare the resulting bias parameters to those derived fixing b0 to the expected value from equation (5). Note that, we expect the local deterministic bias to describe the mapping between matter and haloes only on sufficiently large smoothing scales.

The outcome of this analysis is presented is Fig. 2 and3, which show the results of this fitting process for respectively b0, b1 and c2, c3 when applied at various smoothing scale R. For each scale, we estimate the error by performing the fit in 64 Jackknife sub-samples of the full simulation. Examination of the top panel of Fig. 2 shows that b0 changes with scale, while converging to zero. Furthermore its sign depends on the mass of the halo sample. Both effect are expected, as shown by equation ( 5), the sign of b0 can change, depending on the higher order bias coefficients. Moreover, since b0 depends on quantities such as the variance σ 2 and the skewness S3 of matter, we expect a scale dependence. Regarding the biasing coefficients, Figures 2 and3 show that we find a strong scale dependence of our bias measurements for small smoothing radii R, while results converge for scales larger than R = 40 h -1 Mpc (see also Manera and Gaztañaga 2011). In general, we can also conclude that the choice of leaving b0 free or not does not impact the fitted bias parameters. In addition, looking at the two low mass samples we can see that above 40h -1 Mpc, the choice of the order of the fitted polynomial does not impact significantly the measurements. Despite, the visual difficulty of interpreting the results for the two high mass samples we will show later that this choice can have an impact on the fitted third-order bias parameter. We therefore Figure 1. Scatter plot between halo and matter density contrast at z = 0 obtained with a smoothing radius R = 60h -1 Mpc. We represent the 68% (dark colour) and 95% (light colour) contour levels, in red we show the contour levels corresponding to the mass bin M3 while in blue we show the mass bin M0. Black dashed lines show the bias function derived from the δ hδm method, while the green long dashed, red dotted and magenta dotted lines show results from variants of the δP DF method (see text for details).

only use scales between 40 R 80 h -1 Mpc to derive the bias parameters via the δ h -δm method.

The second method for measuring bias parameters from the δ h -δm relation consists in matching the probability distribution functions (PDF) of δ h and δm, it will therefore be referred to as the δPDF method in the following. As the δ h -δm method, it relies on the assumption that the bias function is local and deterministic. However, the δPDF method presents a clear advantage: it can be directly applied to measure galaxy bias from observations. The PDF of matter fluctuations can be, indeed, modelled for a given cosmology and the PDF of galaxies can be directly observed from redshift surveys (see Bel et al. 2016;Di Porto et al. 2014;Marinoni et al. 2005;Bernardeau et al. 2002;Uhlemann et al. 2016). Inferring the bias from the δ h -δm method, on the other hand, requires adequate observations for measurements of δm via weak gravitational lensing and introduces additional uncertainties (Chang et al. 2016;Pujol et al. 2016).

Thus, measuring the bias parameters by matching the 1-point probability distribution of haloes and matter fluctuations (P h (δ h ) and Pm(δm) respectively) represents an alternative way of estimating local bias parameters. It consists in assuming the existence of a local mapping between the matter and haloes fluctuations such that the probability is conserved under a change of variable following this functional relation In both panels, we estimate quantities with respect to the smoothing scale R. We test several methods: in the first we fix b 0 (solid lines) using equation ( 5), then we allow b 0 to vary as a free parameter (short dashed lines), after we fit a fourth (dashed dot-dot-dot lines) and finally a fifth (long dashed lines) order polynomials. In the bottom panels we also represented the results of the fit (magenta) of the linear bias b 1 between 40 and 80h -1 Mpc. We perform this analysis for two comoving outputs at redshift 0.0 (left) and 0.5 (right). In all panels we adopt a colour coding blue, green, orange and red referring respectively to the halo samples M0, M1, M2 and M3. Note that, for clarity, each method applied to a given halo mass is represented with a different shade of the corresponding colour.

C h [δ h ] = Cm[δm], (6)
where Ci is the cumulative distribution function, defined as

Ci[δi] ≡ δ i -1 Pi(δ ′ i )dδ ′ i for either haloes (i = h) or matter (i = m).
The halo density fluctations δ h can thus be expressed as a function of halo density fluctuations δm by inverting equation (6

) (i.e. δ h (δm) = C h -1 [Cm[δm]]
), which is by definition the local bias function. The exponent -1 denotes the reciprocal function. The mapping C -1 h [Cm] is obtained by integrating numerically the probability distribution function of both the halo and matter density contrasts. For technical reasons, we use a cubical smoothing of size 64h -1 Mpc which corresponds in volume to a spherical smoothing of radius R ≃ 40h -1 Mpc (which, we saw, is large enough to consider the bias coefficients as scale independent). Note that when the average number of haloes per cell is smaller than 150 (which is the case for the samples M2 and M3), it is necessary to correct from shot noise effects Figure 3. Top: Second-order bias coefficient c 2 . Bottom: Thirdorder bias coefficient c 3 . We adopt the same colour and line coding described in Fig. 2 (see Di Porto et al. 2014). For reconstructing the halo PDF we use the relation between the discrete probability distribution PN (where N represents the number of haloes inside a cell) and the probability density function of halo density fluctuations P h (δ h ),

PN = ∞ -1 K[N | N (1 + δ h )]P h (δ h )dδ h , ( 7
)
where K is the conditional sampling probability which we assume to be a Poisson distribution (Layzer 1956). In practice the inversion of equation ( 7) requires some subtleties (see Bel et al. 2016). Following Bel et al. (2016), we apply three different reconstruction methods and choose the one which provides the best description of the counting probability distribution PN . Given the high shot noise level in the M 2 and M 3 samples we assume a parametric form for the density probability P h , in practice we tested a skewed-LogNormal (Colombi 1994) and a Gamma expansion (Gaztañaga et al. 2000) and found that the Gamma expansion offers a better matching of the measured counting probability of haloes.

Once we obtained the bias function we devised two ways of estimating the bias parameters. In both cases we fix an effective range on which the Taylor expansion must be valid and on one hand we fit the bias function with a third-order polynomial (referred to as proba. ( 1)) and on the other hand we apply discreet derivatives method (referred to as On one hand, we display the bias coefficients (b 1 , c 2 and c 3 from top to bottom) obtained by fitting the scatter plot with increasing order of the polynomial (black diamonds, green filled circles and orange empty circles) on scales between 40 and 80h -1 Mpc with respect to the average mass of the halo samples. On the other hand, we show the obtained bias coefficients from the two variants of the δPDF method (see text for more details).

proba. ( 2)). We allow the effective range in δ to decrease and stop when the two methods converge to the same values for b1, c2 and c3. An illustration of this technique is displayed in Fig. 1, where we show in green long dashed line the bias function obtained from inversion of equation ( 6) together with the Taylor expansion estimated from the fitting method (orange dotted lines) and from the numerical derivatives method (magenta dotted lines). One can see that the results from the two methods agree well with the obtained bias function from the probability density function. On the other hand, we can see some significant departures from the Taylor expansion obtained from the direct fit to the scatter plot (black lines). In fact, considering Fig. 4 which summarizes the measurements of the bias coefficients from the two presented methods one can remark that despite small discrepancies the overall behaviour of the bias coefficients with respect to the mass is qualitatively in good agreement. In Fig. 4 we also show, with more clarity, the impact of the choice of the order of the fitted polynomial to the scatter plot: it becomes clear that in the halo samples M2 and M3 there is a non-negligible (given the estimated error bars) effect on the estimation of the c3 coefficient while the lowest order bias coefficients are not significantly affected.

two-point cross-correlation

The bias function in equation ( 3) has an important application in the cosmological analysis of galaxy clustering statistics, since it provides a parameterisation of the deviation between the clustering of haloes (or galaxies) and the un-derlying matter distribution. In simulations, where the clustering of matter can be directly measured, one can infer the bias parameters by comparing the clustering statistics of haloes and matter. In this subsection we focus on bias measurements from the two-point correlation ξ as statistical measure for the clustering. It directly probes the linear growth factor, thereby providing strong constraints on cosmological models. The two-point correlation can be defined in configuration space as the mean product of density fluctuations δi at the positions r1 and r2 that are separated by the distance r12 ≡ |r1 -r2|,

ξxy(r12) ≡ δx(r1)δy(r2) = δ 1 x δ 2 y (r12). ( 8
)
The indices x and y refer to the density fields of haloes (or galaxies) and matter, while . . . denotes the average over pairs of any orientation and position. The two-point auto-correlation (hereafter referred to as 2pc) corresponds to the case x = y and will be denoted in the following as ξ h ≡ ξ hh and ξm ≡ ξmm for haloes and matter respectively, whereas the halo-matter cross-correlation (hereafter referred to as 2pcc) will be denoted as ξ × ≡ ξ hm . Inserting the local bias model from equation (3) into equation ( 8) provides, at leading order, a simple relation between ξm and ξ × ,

ξ × (r12) ≃ b1 ξm(r12) + O[ξ 2 m ]. (9)
At large scales (r12 20h -1 Mpc), where the two-point correlation function is small, we expect O[ξ 2 m ] to be negligible (see e.g. paper I), which allows for measurements of the linear bias as

b × ξ (r12) ≡ ξ × (r12) ξm(r12) ≃ b1. (10)
We showed in paper I that b × ξ is a reliable estimate of the linear bias b1, since it agrees with various other estimators from second-order one-and two-point statistics at large scales.

three-point auto-and cross-correlations

Third-order statistics is sensitive to the shapes of largescale density fluctuations. It therefore provides information about the large-scale structure which is not accessible with second-order one or two-point statistics as the latter is defined isotropically in configurations space. Combining second-and third-order statistics therefore allows for breaking the degeneracy between the linear growth factor and the linear bias parameter which otherwise limits the accuracy of cosmological constraints derived from the 2pc at large scales (see e.g. Bernardeau et al. 2002, or paper I). However, this approach relies on the accuracy with which the linear bias can be measured with third-order statistics. The most general third-order statistics in configuration space is the three-point correlation or its reduced counter part (Groth and Peebles 1977). In this subsection we describe our bias measurements from paper I, which are derived from the symmetrised reduced three-point halo-matter-matter crosscorrelation in configuration space (hereafter referred to as 3pcc, see also Pollack et al. 2012). It is defined as

Q × ≡ 1 3 ζ hmm + ζ mhm + ζ mmh ζ hm H , (11)
where

ζ xyz ≡ δx,1δy,2δz,3 (12)
is the three-point halo-matter-matter cross-correlation and

ζ hm H ≡ ξ hm 12 ξ hm 13 + ξ mh 12 ξ hm 23 + ξ mh 13 ξ mh 23 ( 13
)
is the corresponding hierarchical three-point crosscorrelation. The average . . . is made over triangles of any orientation and position. The auto-correlation for halo and matter density fields (Q h and Qm respectively) are defined analogously and will be referred to as 3pc in the following. For measuring the bias accurately using the 3pcc one needs to take into account non-local contributions to the bias function. Neglecting the latter would cause a ∼ 20% error in the estimation of the linear bias parameter (see Manera and Gaztañaga 2011;Chan et al. 2012;Baldauf et al. 2012). We therefore proposed in paper I to measure the linear bias independently of quadratic local and non-local contributions by combining three-point autoand cross-correlations

, b∆Q ≡ -2 Qm Q h -3Q × h . (14)
A similar combination of auto-and cross-correlations allows for the determination of the quadratic local and non-local bias parameters (c2 and γ2) independently of Qm, once the linear bias is determined (for instance via equation ( 14) or the 2pc),

∆Qcg ≡ Q h -Q × = 2 3 1 b1 [c2 + g2Q nloc ] . (15)
We predict the non-local contribution Q nloc from the nonlinear power spectrum, which has been measured in the simulation. Details on this prediction are provided in the appendix of paper I. Combining three-point auto-and crosscorrelation for measuring bias will be referred to as ∆Q method in the following. The results in this work are derived from 3pc measurements which are computed using triangles with fixed legs of 36 and 72h -1 Mpc at 18 opening angles between these legs. The density fields where smoothed with an 8h -1 Mpc Top-hat window function, and the covariance between measurements at different opening angles was estimated from 64 JK samples. Details on these measurements are presented in paper I, where we also study the dependence of 3pc bias measurements on the triangle scale.

third-order cross-moments

Alternatively to the 3pc one can explore the third-order statistical properties of cosmic fields with one-and two-point statistics. Despite the fact that they do not provide access to the full third-order hierarchy probed by the 3pcc (Subsection 3.3) they are useful tools for comparing the statistical properties of the halo distribution with respect to the underlying matter field. These one-and two-point statistics correspond to the 3pc for very specific triangle configurations. The 3pc for triangles which are collapsed into two points corresponds to the correlator C12 (see Bernardeau 1996), while for a further collapse into a single point yields the skewness S3 (see Bernardeau et al. 2002). As in case of the 3pcc, it is possible to define the cross-skewness S × 3 and the cross-correlator C × 12 as

S × 3 ≡ δ h δ 2 m δ h δm 2 (16)
and

C × 12 ≡ δ h,1 δ h,2 δm,2 δ h,1 δ h,2 δ h δm (17)
where δ h and δm are respectively the density contrast of halo and matter density fields which are smoothed with a spherical top-hat window function of radius R. The auto skewness S3 and the auto-correlator C12 of the matter density field are defined analogously. We provide details on how we measure these quantities and how to correct them for shot-noise in paper I. One can express the linear and quadratic bias parameters independently from each other using the combinations

b × τ ≡ S3 -C12 S × 3 -C × 12 ≡ τ τ × (18)
and

c × τ ≡ S3C × 12 -C12S × 3 τ × . ( 19
)
We will refer to this way of measuring bias as the τ × -method in the following. For spherically averaged quantities such as skewness, correlators, cross-skewness and cross-correlators non-local contributions of the second-order described by equation ( 4) lead to an effective second-order local bias

c eff 2 = c2 - 4 3 γ2 b1 , (20)
but do not affect the estimated linear bias parameter b1 unless for highly biased samples (M3, see discussion in paper I). The same effect has been pointed out by Chan et al. (2012) for the δg -δm method for measuring bias, which is described in Section 3.1. Assuming a local Lagrangian bias model, γ2 = -2 7 (b1 -1) we expect the non-local contributions to c eff 2 to be zero for b1 = 0 and 0.3 for b1 = 4. However these expectation will not be accurate as moderate deviations from the Lagrangian bias model have been reported by Chan et al. (2012) and in paper I.

Peak-background split predictions

The peak-background split (hereafter referred to as PBS) model provides predictions for the linear and non-linear bias parameters as a function of halo mass, which are deduced from derivatives of the mass function. Model fits to mass function measurements therefore provide analytical expression of the bias as a function of halo mass. Such PBS bias predictions are essential for various purposes like bias modelling as function of galaxy properties via Halo Occupation Distribution (HOD) modelling, cluster mass calibration or predictions of the cluster count and power spectra covariances (e.g. Cooray and Sheth 2002;Lima andHu 2004, 2005;Lacasa and Rosenfeld 2016). It is therefore important to verify how accurately the PBS model is able to predict bias measurements in simulations.

In the present analysis we validate the accuracy of PBS predictions for the linear, quadratic and third-order bias parameters by comparing them to the various measurements described in the previous subsections. We employ fits to the measured MICE-GC mass function from paper II, based on the models of Tinker et al. (2010) and Warren et al. (2006) as well as a new model presented in paper II, to which we refer to as Tinker, Warren and HBG15 model respectively in the following. The sample M0 has been excluded from the fitting range since we expect halo detection in the corresponding mass range to be very noisy, as we discuss in paper II. However, the predictions can still be made over the full mass range for all the samples M0-M3. For studying the universality of the relation between the different bias parameters we adopt the fitting parameters for the Tinker and Warren model, which are provided in the corresponding articles. Note that these latter fits are based on simulations with different cosmologies than in the MICE-GC simulation. Furthermore Tinker et al. (2010) define haloes as spherical over-densities, not as FoF groups as done by Warren et al. (2006) and in the present study.

RESULTS

bias comparison

We compare in Fig. 5 the linear and non-linear bias parameters from the different measurements and the PBS predictions described in Section 3 using the mass samples M0-M3 (defined in Table 2) at the redshifts z = 0.0 and z = 0.5.

linear bias

In the case of the linear bias (b1, top panel of Fig. 5) we find an overall variation of roughly 5 percent between results from different measurements and PBS predictions at all masses and both redshifts. The measurements from the 2pcc (b × ξ , equation ( 10)) are expected to be the most reliable estimate with percent level accuracy (see e.g. paper I). We therefore use it as a reference for evaluating the accuracy of the other linear bias measurements and predictions in our comparison.

The results in Fig. 5 show that both estimations for the linear bias from third-order statistics, (∆Q and τ × , described in Section 3.3 and 3.4 respectively) are in 1σ agreement with b × ξ . However, the errors of b × τ are roughly 10 percent, while those of b∆Q at the one percent level. We attributed the strong deviation of b × τ from all other results at high masses to an inadequate shot-noise correction in paper I. From the same analysis we expect larger deviations from b × ξ when the linear bias is measured from third-order autoinstead of cross-correlations correlations and when non-local contributions are not taken into account. The two types of linear bias measurements from the δ h -δm relation (described in Section 3.1) are both in percent level agreement with b × ξ which is consistent with results from Manera and Gaztañaga (2011) or Pollack et al. (2012). Results from the δPDF method tend to be overall slightly larger with 7% deviations from b × ξ . A possible reason for the discrepancies can be inaccuracies in the covariance estimation from Jackknife sampling, which might affect the best fit values. Note that we do not expect the discrepancies to result from the truncation of the bias function at third-order based on our convergence tests in Section 3.1.

In addition to the bias measurements we show predictions from the PBS model which are based on different fits to the MICE-GC mass function, presented in paper II. The different predictions show a very good mutual agreement in the mass range spanned by the samples M1-M3 which was used for the fits. For the low mass sample, which was excluded from the fitting range due to noisy halo detection, the different PBS predictions show a stronger variation. For the mass samples M2 and M3 ( 10 13 h -1 M⊙) all predictions lie up to ∼ 7 percent below b × ξ , which is consistent with findings in the literature (e.g. Manera et al. 2010;Manera and Gaztañaga 2011;Pollack et al. 2012;Paranjape et al. 2013a;Lazeyras et al. 2016).

The agreement between predictions from different mass function fits indicates that the inaccuracy of PBS bias is driven by shortcomings of the model, rather than uncertainties in the fits. Such a shortcoming could be the assumption of a constant matter density threshold for the gravitational collapse on which the standard PBS model application to mass functions is based on, as pointed out in recent studies. In particular Paranjape et al. (2013a) showed that the scale-independent PBS bias parameters, reconstructed from scale dependent Lagrangian bias measurements in simulations, are in good agreement with predictions from the excursion set peak (ESP, Paranjape et al. 2013b) model, when a halo mass-dependent scatter in the collapse threshold is included. A percent level agreement between such ESP predictions for PBS bias and direct measurements of the latter from separate universe simulations (which do not rely on a threshold model) was reported by Lazeyras et al. (2016). These direct PBS bias measurements are also in excellent agreement with bias measurements from large-scale Fourier space clustering, confirming results from Baldauf et al. (2015) and Li et al. (2015).

non-linear bias

Measurements and prediction of the second-and third-order bias parameters (c2 ≡ b2/b1 and c3 ≡ b3/b1) are shown in the central and bottom panels of Fig. 5 respectively. Note that for this comparison we do not have a reference for the non-linear bias, such as measurements from separate universe simulations. Furthermore our measurements from third-order statistics provide only bias parameters up to second-order, since they are based on leading order perturbation theory. Measurements of c3 are therefore only obtained from the δ h -δm relation.

Overall we find a stronger variation between the different results for the non-linear bias than for the linear bias. Yet, all the different non-linear bias measurements and predictions show similar mass dependencies at both redshifts. The second-order bias c2 depends only weakly on the halo mass below 10 13 h -1 M⊙ (covered by the mass samples M0 and M1), with a value of ≃ -0.5 at both redshifts. In the mass range ≃ 10 13 -10 14 h -1 M⊙, which corresponds to the mass sample M2, c2 changes from negative to positive values, and increases rapidly with mass to values of up to 2 at z = 0.5 for halo masses of higher than 10 14 h -1 M⊙. The PBS prediction for c2 from all mass function fits show an overall weaker mass dependence than the corresponding measurements, as the former lie above the results from the ∆Q, τ × and δ h -δm methods in the mass range spanned by M0-M2, where c2 is negative. For the high mass sample M3, where c2 is positive at both redshifts we find the opposite. The latter result confirms findings of Manera and Gaztañaga (2011) and Lazeyras et al. (2016) on a larger mass range. An overly weak mass dependence of standard PBS predictions has also been found for Lagrangian bias by Paranjape et al. (2013a).

In the mass range M0-M2 the measurements from the ∆Q, τ × and δ h -δm methods agree mutually at the 1σ level at both redshifts, while results vary strongly for the high mass sample M3. Note that we expect the c2 measurements from the δ h -δm and τ × to be biased by non-local contributions with values between zero for b1 = 0 and 0.3 for b1 = 4, as discussed in Section 3.4, which is comparable with the variation among the different results.

In the bottom panel of Fig. 5 we present a validation of third-order bias prediction from the PBS model with direct measurements from the δ h -δm and δPDF methods. The predictions and measurements are in an overall agreement with each other, as both decrease from positive values of around unity for M0 to negative values of down to -3 for M3 at both redshifts. The zero crossings of c3 occur between 10 12.5 -10 13.5 M⊙/h, which is consistent with the direct PBS bias measurements from separate universe simulations from Lazeyras et al. (2016). Overall the predictions tend to be smaller than the measured values. Note that we used thirdorder polynomials for the bias function in order to measure the bias parameters. From our convergence test in Section 3.1, we expect larger deviations from PBS predictions in the high mass bins M2 and M3 when the third-order bias is measured using higher order polynomials.

Universal relation between bias parameters

Our various bias measurements enable us to validate the universal relation between the linear and non-linear bias parameters b2(b1) and b3(b1), which we found for PBS predictions and linear bias values of b1 2 in paper II. Such a universal behaviour can be expected from a universality of the mass function for different redshifts and cosmologies.

In Fig. 6 we show the linear bias versus the quadratic and third-order bias parameters from the different measurements together with the PBS predictions for redshift z = 0.0 and 0.5.

The PBS predictions labeled as Crocce et al. (2010) are derived from FoF mass function fits of these authors to nested boxes runs from the MICE simulation suite which provides a higher mass resolution in the low mass range. The PBS predictions from the mass functions from Warren et al. (2006), Tinker et al. (2010) and Watson et al. (2013) are based on simulation with cosmologies different from the one of MICE. In addition, Tinker et al. (2010) defined haloes as spherical over-densities instead of FoF groups. These difference can lead to deviations in the mass function (see paper II) and therefore contribute do the differences which we see among the various PBS bias predictions for the bN (b1) relation in Fig. 6. Note that the BPS results from Lazeyras et al. (2016), which also use spherical over-densities as haloes, have not been predicted from the mass function but were directly measured using separate universe simulations. It is interesting to note that the b2(b1) relation from Lazeyras et al. (2016) agrees best with our measurements derived from the ∆Q method, since the latter also delivers very accurate measurements of the linear bias. This suggests that also the c2 measurements and hence the b2 -c2 relations from these two approaches are reliable. Furthermore we include bias measurements from Chan et al. (2012), which are based on the Bispectrum and find a good agreement with our results.

Overall we find consistent results from the different measurements and predictions from simulations with different cosmologies analysed at different redshifts. This indicates a roughly universal behaviour between the linear and higher-order bias parameters. In paper II we derived an analytic expression for bN (b1) from the PBS bias parameters based on the Press and Schechter (1974) (21)

Since this expression is independent from the peak-height, we also do not expect it to depend on how halo masses are defined. The bN (b1) predictions from paper II shown in Fig. 6, have been fitted to PBS predictions for b2(b1) and b3(b1) which are based on the the Tinker mass function model, fitted to the MICE-GC mass function at z = 0.0 in the same article. Fitting the polynomial from equation ( 21) to measurements from ∆Q we find (α0, α1, α2) = (0.77, -2.43, 1).

SUMMARY AND CONCLUSION

This analysis is the last part of a series of articles on the accuracy of bias parameters derived from the MICE Grand Challange (MICE-GC) simulation using clustering statistics and peak-background split (referred to as PBS) predictions (Hoffmann et al. 2015b, paper I, paper II). In the present analysis we studied bias parameters derived from the relation between matter and halo density contrasts (δ h , δm respectively) as an additional method for measuring bias. These measurements are compared to a selection of our most robust previous results using the same four mass samples M0-M3 at redshift 0.0 and 0.5. Thanks to the large volume and resolution of the MICE-GC simulation these samples span a large mass range from between roughly 10 12 and 10 15 h -1 M⊙, which corresponds to Milky Way like haloes and massive galaxy clusters respectively. Our previous have been derived from two-point halomatter cross-correlations (ξ × ), a combination of three-point halo-matter auto-and cross-correlations (∆Q) as well as a combination of the halo-matter cross-skewness and crosscorrelators (τ × ). We thereby employ leading-order modelling of clustering statistics, at which the linear bias parameters from these estimators are not affected by non-local contributions to the bias model. We therefore obtain the linear bias from ξ × and the linear and quadratic bias from ∆Q and τ × . The PBS predictions are based on MICE-GC mass function fits from paper II, using different mass function models, while we study bias parameters up to order three, for which we have corresponding measurements from the δ h -δm relation for validation.

We studied in this work bias measurements from the δ h -δm relation in two ways. The more common method is to fit a polynomial to the δ h -δm relation, measured in the simulation. Alternatively we explore bias measurements obtained from the probability distribution function (referred to as PDF) of δ h and δm. The latter method has the advantage that it can be directly applied to observations, since the PDF of matter density contrasts can be modelled with theory or simulations (Bel et al. 2016;Di Porto et al. 2014;Marinoni et al. 2005;Bernardeau et al. 2002). However, it has not been tested so far how accurately the bias parameters can be determined with this method.

The results of our bias comparison are summarized in Fig. 5. In the case of the linear bias b1 we consider results from ξ × (b × ξ ) as the most reliable since non-linear local and non-local term can safely be neglected. Furthermore, b ξ is highly relevant for cosmology since it is weakening constraints on cosmological parameters inferred from ξ due to its degeneracy with the growth of matter fluctuations.

Linear bias measurements and predictions from all other methods are in a ∼ 5 percent agreement with b × ξ , while ∆Q delivers the most accurate measurements with an overall percent level accuracy. Standard measurements from the δ h -δm relation are in a slightly better overall agreement with b × ξ than those from the δPDF-method. The strongest deviations from b × ξ shows the linear bias predictions from the different PBS model. These predictions are consistently up to 7 percent below b × ξ at high masses as reported in the literature (e.g. Manera and Gaztañaga 2011;Pollack et al. 2012). The fact that we find similar results for different mass function models indicates shortcomings in the standard PBS modelling, i.e. the assumption of a constant matter density threshold for gravitational collapse (Paranjape et al. 2013a,b;Lazeyras et al. 2016). We do not find a clear change of the variation between the different measurements and predictions with mass or redshift.

In the case of the quadratic bias c2 ≡ b2/b1 we find consistent results from the different measurements and predictions as c2 increases from negative values of -0.5 at low halo masses (∼ 10 12 M⊙/h) to positive values of up to ∼ 2 at high masses ( 10 14 M⊙/h). However, the variation between the different results tends to be larger than in the case of b1. In the case of the measurements this effect is presumably caused by the strong assumptions such as the validity of tree-level perturbation theory (∆Q and τ × ), Poissonian shot-noise (τ × ), or a local deterministic bias model (δ h -δm). We therefore do not consider any of these measurements to be sufficiently reliable for being a reference, such as the results from ξ × in the case of the linear bias. However, the fact that the ∆Q method delivers highly accurate measurements of the linear bias suggests that also the quadratic bias is measured reliably with this approach. We find the PBS predictions for c2 to be consistently above (below) all measurements in the low (high) mass range as they show overall weaker mass dependence. The latter finding lines up with reports from Manera and Gaztañaga (2011) and Pollack et al. (2012).

Both measurements of the third-order bias c3 ≡ b3/b1 from the δ h -δm relation agree overall mutually at the 1σ level. These measurements allow for a validation of the corresponding PBS predictions. Results from both methods are similar as c3 is positive below unity in the low mass range and decreases to negative values of down to ∼ -3 in the high masses range. However, deviations are significant, as the PBS predictions tend to be below the measurements at low halo masses while results based on different mass function fits are consistent with each other.

We use our various linear and non-linear bias measurements for validating the universal polynomial relation between linear and non-linear bias (b2(b1) and b3(b1)), which we deduced in paper II from PBS predictions based on the Press and Schechter (1974) mass function. Since this expression is independent of the peak-height we do not expect a strong dependence on halo mass definition. Our measure-ments show an overall agreement with the universal behaviour predicted by the PBS model. Furthermore they agree with results from the literature derived via the Bispectrum in Fourier space or the separate universe approach from simulations with cosmologies different to the one of MICE (i.e. Chan et al. 2012;Lazeyras et al. 2016). We fit a secondorder polynomial to b2(b1) measurements from ∆Q, which we consider as the most reliable c2 estimator as mentioned above.

Such a universal relation between linear and non-linear bias can be useful for reducing errors on the linear bias and the growth from clustering analysis when the latter is affected by c2, for instance in the case of three-point correlations or two-point correlations at small scales.

For applying universal polynomial relations between bias parameters in the analysis of galaxy surveys it would be interesting to show that their universality also holds for halo samples, which are selected by galaxy properties, such as luminosity and colour instead of halo mass. A correlation between the linear and quadratic bias from the 3pc and Bispectra has been reported by for SDSS galaxy samples and mock HOD catalogues by Kayo et al. (2004) and Nishimichi et al. (2007) and compared with PBS predictions. The limited accuracy of their measured b1 -c2 relations could be strongly improved using the methods studied in the present analysis in combination with data from upcoming galaxy surveys.

Figure 2 .

 2 Figure 2. Top: Bias coefficient b 0 . Bottom: Linear bias coefficient b 1 .In both panels, we estimate quantities with respect to the smoothing scale R. We test several methods: in the first we fix b 0 (solid lines) using equation (5), then we allow b 0 to vary as a free parameter (short dashed lines), after we fit a fourth (dashed dot-dot-dot lines) and finally a fifth (long dashed lines) order polynomials. In the bottom panels we also represented the results of the fit (magenta) of the linear bias b 1 between 40 and 80h -1 Mpc. We perform this analysis for two comoving outputs at redshift 0.0 (left) and 0.5 (right). In all panels we adopt a colour coding blue, green, orange and red referring respectively to the halo samples M0, M1, M2 and M3. Note that, for clarity, each method applied to a given halo mass is represented with a different shade of the corresponding colour.

Figure 4 .

 4 Figure4. On one hand, we display the bias coefficients (b 1 , c 2 and c 3 from top to bottom) obtained by fitting the scatter plot with increasing order of the polynomial (black diamonds, green filled circles and orange empty circles) on scales between 40 and 80h -1 Mpc with respect to the average mass of the halo samples. On the other hand, we show the obtained bias coefficients from the two variants of the δPDF method (see text for more details).

  cross-correlation ξ × b 1 Section3.2, paper I combined three-point auto-and cross-correlations ∆Q b 1 , c 2 Section3.3, paper I combined third-order auto-and cross-moments τ × b 1 , c 2 Section3.4, paper I predictions peak-background split (PBS) model b 1 , c 2 , c 3 Section3.5, paper II

  mass function, which has the form bN = N=2 n=0 αnb n .

Figure 5 .

 5 Figure5. Summary of various bias measurements from second-and third-order halo-mater cross-correlations (ξ × , ∆Q and τ × ), δ hδm relations and PBS predictions, described in Section 3. The top, central and bottom panel show the first-second-and third-order bias parameters respectively at z = 0.0 and 0.5 (left and right panel) versus the mean halo mass of each mass samples M0-M3 (slightly shifted along the mass axis for clarity). The subpanel at the bottom of the top panel shows relative deviations from ξ × . The lower and upper limits of the mass samples are marked by vertical grey dashed lines. Error bars denote σ uncertainties.

Figure 6 .

 6 Figure 6. Second-and third-order bias parameters versus the linear bias parameter b 1 (top and bottom panel respectively). Results from ∆Q and τ × ), δ hδm as well as PBS predictions from paper II are the same as those shown in Fig. 5. In addition we show PBS prediction based on mass function fits from the literature. Fits to the b 1 -c 1 and b 1 -c 3 relations from separate universe simulations in combination with the PBS model, given by Lazeyras et al. (2016) are shown as long dashed blue lines.

Table 1 .

 1 Abbreviations for previous articles of this series.

	name in the text	reference
	paper I	Bel et al. (2015)
	paper II	Hoffmann et al. (2015a)

Table 2 .

 2 Halo mass samples used in this study. Np is the number of dark matter particles per halo, N halo is the number of haloes per sample in the comoving output at redshift z = 0.5.

	sample	mass range [10 12 h -1 M ⊙ ]	Np	N halo
	M0	0.58 -2.32	20-80	122300728
	M1	2.32 -9.26	80-316	31765907
	M2	9.26 -100	316-3416	8505326
	M3	100	3416	280837

Table 3 .

 3 Bias estimators, definitions and notations. The term cross refers to halo-matter cross correlations, b and c indicate linear and non-linear bias estimations respectively.

ACKNOWLEDGEMENTS

Funding for this project was partially provided by the Spanish Ministerio de Ciencia e Innovacion (MICINN), project AYA2009-13936, AYA2012-39559 and AYA2015-71825, Consolider-Ingenio CSD2007-00060, European Commission Marie Curie Initial Training Network CosmoComp (PITN-GA-2009-238356) and research project 2014 SGR 1378 from Generalitat de Catalunya. KH is supported by beca FI from Generalitat de Catalunya and MINECO, project ESP2013-48274-C3-1-P. He also acknowledges the Centro de Ciencias de Benasque Pedro Pascual where parts of the analysis were done. The MICE simulations have been developed by the MICE collaboration at the MareNostrum supercomputer (BSC-CNS) thanks to grants AECT-2006-2-0011 through AECT-2010-1-0007. Data products have been stored at the Port d'Informaci Cientfica (PIC).

We thank Martin Crocce, Pablo Fosalba, Francisco Castander and Roman Scoccimarro for interesting and useful comments.