# Comparing halo bias from abundance and clustering

###### Abstract

We model the abundance of haloes in the volume of the MICE Grand Challenge simulation by fitting the universal mass function with an improved Jack-Knife error covariance estimator that matches theory predictions. We present unifying relations between different fitting models and new predictions for linear () and non-linear ( and ) halo clustering bias. Different mass function fits show strong variations in their performance when including the low mass range () in the analysis. Together with fits from the literature we find an overall variation in the amplitudes of around % in the low mass and up to % in the high mass (galaxy cluster) range (). These variations propagate into a % change in predictions and a % change in or . Despite these strong variations we find universal relations between and or for which we provide simple fits. Excluding low mass haloes, different models fitted with reasonable goodness in this analysis, show percent level agreement in their predictions, but are systematically % lower than the bias directly measured with two-point halo-mass clustering. This result confirms previous findings derived from smaller volumes (and smaller masses). Inaccuracies in the bias predictions lead to % errors in growth measurements. They also affect any HOD fitting or (cluster) mass calibration from clustering measurements.

###### keywords:

galaxies: haloes – galaxies: abundances – methods: analytical – methods: statistical – dark matter – large-scale structure of Universe## 1 Introduction

Observations of structures in the large-scale distribution of galaxies are a powerful tool for constraining cosmological models. However, such constraints require a model which connects the galaxy distribution to the matter density field. Various observations support the gravitational instability paradigm in which galaxies form in potential wells, generated by the gravitational collapse of dark matter into haloes. The relation between the halo and full matter density fields ( and respectively) is therefore a crucial ingredient for a precise large-scale structure analysis. In fact, the uncertainties in this relation strongly increase the errors in the Dark Energy equation of state or gravitational growth index from future galaxy surveys (e.g. Eriksen and Gaztañaga, 2015).

A formal approach for describing the halo-matter density relation was suggested by Fry and Gaztanaga (1993) as the Taylor expansion around the matter density contrast, , at the same position, known as bias function

(1) |

where , is the mean density and denotes the spatial position. For the construction of the density field is commonly smoothed with a top-hat filter of size . The coefficients are the so called bias parameters, while we will investigate non-linear bias in terms of the ratios and . The relation in equation (1) corresponds to a local bias model in which the density of galaxies is fully determined by the matter density at the same position, while environmental effects are not considered. Recent studies demonstrated that the local model is inadequate as tidal forces of the surrounding large-scale structure generate non-local contributions to the bias function (e.g. Chan et al., 2012; Baldauf et al., 2012). However, the two-point correlation, commonly used to study galaxy clustering, is primarily sensitive to the linear bias parameter at scales between Mpc. Due to the level of precision achieved in our analysis we will not take non-linear and non-local bias contributions to the two-point correlation into account. Note that, especially at smaller scales, such a negligence would not be appropriate (e.g. Saito et al., 2014; Biagetti et al., 2014).

Besides the clustering also the abundance of haloes as a function of halo mass (known as the mass function) is related to the bias function. This relation can be understood with the peak-background split model (hereafter referred to as PBS model, Bardeen et al., 1986; Cole and Kaiser, 1989; Mo and White, 1996). In this model, large-scale density fluctuations are superposed with fluctuations at small scales. These large-scale density fluctuations modulate the background cosmology (i.e. the mean density and the Hubble rate) around small-scale fluctuations (e.g. Martino and Sheth, 2009). The critical density contrast for gravitational collapse therefore depends on the environment. In regions with large-scale overdensities more small-scale fluctuations collapse to haloes than in underdense regions. This effect modifies the abundance of haloes and also their spatial distribution as they follow the pattern of the peaks of large-scale fluctuations. Haloes are therefore biased tracers of large-scale fluctuations in the full matter density field. For a given matter power spectrum the halo bias parameters can be predicted from the mass function via the PBS model.

The PBS bias predictions can be used to determine the dark matter clustering from observed galaxy distributions if the halo masses of a given tracer sample are known (or the other way round). Such an analysis requires that the bias parameters, predicted from the mass function, are equivalent with the bias which affects the clustering. Studies of this equivalence have revealed that the PBS predictions for the linear bias are around below measurements from two-point clustering statistics. Such deviations might result from assumptions of the PBS model, such as spherical collapse, or a local bias relation (e.g. Mo et al., 1997; Desjacques et al., 2010; Paranjape et al., 2013; Schmidt et al., 2013). Further numerical effects, like the definition of haloes in N-body simulations, or systematic effects such as the parametrisation and fitting procedure of the mass function might contribute to the discrepancy between the bias from PBS and clustering (e.g. Hu and Kravtsov, 2003; Manera et al., 2010). Predictions of the PBS for the relation between halo mass and bias are also employed in Halo Occupation Distribution models to predict the bias as a function of galaxy properties, such as luminosity or color (e.g. Cooray and Sheth, 2002; More et al., 2011; Coupon et al., 2012; Carretero et al., 2015). Inaccuracies of the PBS can affect such halo model predictions for galaxy bias or the average number of galaxies per halo. Moreover, haloes of equal mass could have different galaxy occupation, depending on their environment (e.g. Pujol and Gaztañaga, 2014). Besides clustering analysis the PBS can be employed for estimating the lower mass threshold (or mass-observable relation) of observed galaxy samples. This so-called self calibration method (Lima and Hu, 2004, 2005) uses the fact that both, the clustering and the abundance of haloes, depend on halo mass. Inaccuracies of the PBS model can change the estimation of halo mass thresholds and therefore change the cosmological parameters inferred from such an analysis (e.g. Wu et al., 2010; Manera and Gaztañaga, 2011).

The broad application of the PBS model in large-scale structure analysis and the precision of abundance and clustering measurements from incoming observational data calls for a detailed validation of the PBS bias predictions. The purpose of this analysis is to pursue the study of deviations between halo bias measurements from clustering and PBS predictions using the wide mass range of the MICE Grand Challenge (hereafter referred to as MICE-GC) simulation (Fosalba et al., 2015a, b; Crocce et al., 2013; Carretero et al., 2015; Hoffmann et al., 2015). We thereby focus on the effect of mass function parametrisation and fitting on PBS bias predictions. The mass function fits are affected by the error estimations. Our analysis therefore includes a detailed study of the mass function error and covariance which leads us to an improvement of the standard Jack-Knife estimator. The study of PBS bias predictions includes non-linear bias parameters which are important for an analysis of higher-order correlations of the large-scale halo distribution and two-point correlations at small scales. We further compare the mass function fits and bias predictions with results from the literature based on different simulations, to verify a universal behaviour of these quantities.

This paper is organised as follows. In Section 2 we present the MICE-GC simulation and mass function fits. In Section 3 we present new galaxy bias predictions which we compare with the literature and find a universal relation between bias parameters. In Section 4 we compare these predictions with the bias directly measured from the two-point halo-matter cross-correlations of the MICE-GC simulation. The comparison with higher-order clustering will be presented in a separate paper (Bel, Hoffmann & Gaztañaga, in preparation). A summary is given together with our conclusions in Section 5. In Appendix B we present a new method to improve the Jack-Knife covariance matrix estimation. This method can also be easily generalised to other statistics, such as the two-point correlation function (Hoffmann et al., in preparation).

## 2 Simulation & Halo Mass Function

Our analysis is based on dark matter haloes, identified in the comoving outputs of the MICE-GC simulation at the redshifts and . Starting from small initial density fluctuations at redshift the formation of large-scale cosmic structure was computed with gravitationally interacting collisionless particles in a Mpc box using the GADGET - 2 code (Springel, 2005) with a softening length of kpc. The initial conditions were generated using the Zel’dovich approximation and a CAMB power spectrum with the power law index of , which was normalised to be at . The cosmic expansion is described by the CDM model for a flat universe with a mass density of = . The density of the baryonic mass is set to and is the dark matter density. The dimensionless Hubble parameter is set to . More details and validation tests on this simulation can be found in Fosalba et al. (2015a).

Dark matter haloes were identified as Friends-of-Friends groups (hereafter referred to as FoF groups, Davis et al., 1985) with a redshift independent linking length of in units of the mean particle separation. These halo catalogues and the corresponding validation checks are presented in Crocce et al. (2013). To study the galaxy bias as a function of halo mass we divide the haloes into the four redshift independent mass samples M0, M1, M2 and M3, shown in Table 1. These samples span a mass range from Milky Way like haloes up to massive galaxy clusters. In our analysis we consider mass function fits over different mass ranges which we label as M0123, M123, M23 or M012, following the notation in Table 3.

sample | mass range [] | ||
---|---|---|---|

M0 | |||

M1 | |||

M2 | |||

M3 |

### 2.1 Mass function definition and measurement

The unconditional mass function, , is defined as the comoving number density of haloes with masses between and . The mass function can be written in a form which is nearly independent of redshift, cosmology and initial power spectrum (Press and Schechter, 1974; Bond et al., 1991; Sheth and Tormen, 1999) as

(2) |

where is the mean comoving mass density. The height of density peaks is defined as

(3) |

where is the critical density for spherical collapse (which is the exact solution value for the spherical collapse in an Einstein-de Sitter universe). The variance of matter density fluctuations, , smoothed with a spherical top-hat window with radius , can be calculated as

(4) |

where is the spherical top-hat window in Fourier space and is the linear power spectrum. Note that refers to the matter density field when it appears as lower index and to the mass, enclosed by , when used as a variable. We measure the mass function in the MICE-GC simulation at redshift and convert it to , to predict the halo bias parameters , and via the PBS theory (Bardeen et al., 1986; Cole and Kaiser, 1989; Mo and White, 1996). We do not apply the halo mass correction suggested by Warren et al. (2006) for low mass resolution, since we analyse haloes down to particles, while this correction was only proposed for larger numbers of particles per halo. Furthermore, it is not clear that the FoF mass, corrected in such a way is closer to the halo mass on which the PBS model is based on. But note that our results do not depend on this correction as illustrated in Fig. 4. More details about how we measure the mass function are given in Appendix A.

Our measurements of at and are shown as symbols in Fig. 1. As expected, they agree visually with the idea of a weak redshift dependence for FoF mass functions when using a redshift independent linking lengths (e.g. Press and Schechter, 1974; Jenkins et al., 2001; More et al., 2011). Errors and covariances of the measurements were derived with a new estimator which combines the JK approach with predictions for sampling variance from the power spectrum (see Appendix B). We also show in Fig. 1 fits to the measurements, based on the mass function parametrisation of Tinker et al. (2010, equation (5)). The model, fitted over the mass range M123 (that is, excluding the low mass sample M0, see Table 3) is in reasonable agreement with the measurements. Including the sample M0 (haloes with less than particles) to the fitting range leads to poor fits of the model. The fits at both redshifts differ by less than for , confirming low redshift dependence from the measurements. The redshift dependence is stronger when lower masses are included in the fitting range, possibly because of redshift dependent noise in the low mass FoF detection. At larger masses (ln() we find up to deviations, which are comparable with the mass function errors. We have verified that our conclusions also hold for fits over the higher mass range, M23, and different mass function binnings. A detailed analysis of the mass function fits, including fits of other mass function models over different mass ranges and different binnings as well as a comparison with fits compiled from the literature, can be found in the next section.

### 2.2 Mass function fits

In order to predict the halo bias from the mass function via the PBS approach we fit different mass function models to the measurements. Several systematic effects, such as the choice of the mass function model or the mass range over which the model is fitted can limit the accuracy of the PBS bias predictions (e.g. Manera et al., 2010; Manera and Gaztañaga, 2011). The objective of the subsequent analysis is to find out how strongly these effects impact the predicted linear, quadratic and third-order bias coefficients. In particular we aim to verify if the disagreement between PBS predictions for the linear bias and the corresponding measurements from two-point correlations, presented in Section 4, is driven by possible shortcomings of the mass function fits. We therefore study in this subsection the fitting performances of different mass function models.

The latest model in our analysis with the highest number of free parameters is the expression given by Tinker et al. (2010) (hereafter referred to as Tinker model). It can be written as

(5) |

where are the free parameters. We have redefined the parameters so that fixing certain parameters delivers expressions which correspond to the mass function models suggested by Press and Schechter (1974), Sheth and Tormen (1999) and Warren et al. (2006) (hereafter referred to as PS, ST and Warren model respectively). The corresponding parameter constraints are summarised in Table 2 together with the abbreviations for the reference of each model. This unification of notation allows a more direct comparison between models. In Table 2 we also propose a constrain, which constitutes a new mass function fit. Its advantage is that it has as many free parameters as the Warren model, but matches the mass function better when we fit over the whole mass range, as we show later. In our analysis we will focus on the models of ST, Warren, Tinker and our proposal.

model | reference | constraints |
---|---|---|

Tinker | Tinker et al. (2010) | free |

Warren | Warren et al. (2006) | |

ST | Sheth and Tormen (1999) | , |

PS | Press and Schechter (1974) | , , |

proposal | this work |

We determine the best fitting parameters for each mass function model by minimising

(6) |

with and . and are derived from our new JK estimator, introduced in Appendix B. For searching the best fitting parameters we implemented a Monte Carlo Markov Chain algorithm to explore the parameter space.

In Fig. 2 we show the significance of the deviations between the mass function measurements and the best fits by the different models. Results are shown at the redshifts and , while at each redshift we fit the mass function over the different mass ranges, which are shown in Table 3. The first includes all halo mass samples (M0123), the second and the third exclude the low mass samples (M123, M23) and the fourth mass range excludes the highest mass sample (M012). For each fitting range we show fits based on seven different mass function binnings, dividing the mass range into logarithmic bins. We find that the deviations between fit and measurement can vary with the binning. However, we also see trends which are independent of this systematic effect.

All mass function models show a clear dependence of the best fit on the chosen mass range, while this dependence is weaker for the ST model. The Tinker parameterization is the model which best fits the measurements at both redshifts and all mass ranges. This can be attributed to the fact that it contains the highest number of free parameters. The best fit parameters for the Tinker model are given in Table 4. For fits over the whole mass range (M0123) our proposed mass function model seems to match the measurements almost as good as the Tinker model, while having one free parameter less. It also has the advantage of producing stable values for the parameters regardless of the range used for the fit. When the fits are performed only at the highest mass range (M23) the Tinker and the Warren mass functions fit the data equally well, while the proposed model is a slightly worse fit. The ST model delivers the poorest fits in all cases. At we find strong deviations between fits and measurements when the fitting range includes the low mass sample M0. This indicates that the FoF detection of low mass haloes can be strongly affected by shot-noise, while this effect is stronger at higher redshift.

mass range | halo masses [] | |
---|---|---|

M0123 | ||

M123 | ||

M23 | ||

M012 |

mass range | |||||||
---|---|---|---|---|---|---|---|

M0123 | |||||||

M0123 | |||||||

M123 | |||||||

M123 | |||||||

M23 | |||||||

M23 |

For studying the goodness of the best fits for the different mass function models we present their best fit parameters and the corresponding values per degree of freedom () in Fig. 3, where the refer to the number of mass function bins used for the fit. Results are shown for fits over the mass ranges M0123, M123 and M23, which correspond to the different minimum peak heights given by the x-axis. For clarity we show here only results at redshift , while we find similar results at . For each fit we show mean results with standard deviations from the seven mass binnings mentioned previously. In addition to the results derived by taking the covariance between different mass function bins into account in the fitting procedure we show results which were computed neglecting the covariance. We find that neglecting the covariance can lead to different best fit parameters, especially when the low mass range, where the off-diagonal elements of the covariance have the highest amplitudes, is included in the analysis (see Appendix B). However, the bias predictions are only weakly affected by the negligence of the covariance (see Fig. 15). The conclusions of this article about the comparison between bias predictions and measurements does not dependent strongly on the covariance use.

The results, shown in the bottom panel of Fig. 3, are very high when the mass functions are fitted over the whole range (lowest minimum peak height). This poor performance, which is even apparent for the Tinker model with its five parameters, is probably related to the fact that our mass function measurements are not reliable in the low mass range. In fact the M0123 sample includes haloes with down to particles. For such low numbers of particles per halo we expect strong systematic effect in the halo mass estimation and therefore on the mass function (e.g. Warren et al., 2006; More et al., 2011). Furthermore, the halo samples might be contaminated with spuriously linked FoF groups. If the analysis is performed using only the high mass sample M23 (highest minimum peak height), the values for the best fit models drop down to values between unity and four. If we perform the fits ignoring off-diagonal elements in the covariance matrix we obtain substantially lower values, especially when the fits are performed over the whole mass range. This demonstrates that the covariance cannot be neglected in the fit for the evaluation of the fitting performance of a mass function model. This statement is even true in the high mass range where the covariance is dominated by shot-noise. This is important as the goodness of the fit is the way to validate the predictions.

We also see in Fig. 3 that the can change for different mass function binnings, which can already be seen in Fig. 2. This dependence on the binning is also apparent when the off-diagonal elements of the covariance matrix are neglected in the fit. However, the best fit values of each model and the corresponding bias estimations are only weakly affected by this systematic effect.

Interestingly the best fit parameters of the Tinker model have the same values as the ones from the Warren model when the fit is performed on the higher mass M23 sample. Consequently the minimum are the same in both cases. This indicates that the parameter , which is set to zero in the Warren model is not required for fitting the high mass range, but becomes necessary, when the low mass range is included in the fit. The values of our proposed model are smaller than those for the Warren model for minimum peak heights of (M123). This agrees with the visual impression, gained from Fig. 2, that our proposed model delivers better mass function fits than the model of Warren, unless the analysis is restricted to the highest mass range M23. We come to the same conclusion when analysing the mass function at .

### 2.3 Mass function universality

In Fig. 1 we demonstrated that the mass function, when expressed in terms of the peak-high , depends only weakly on redshift. To verify if this universality also holds for other cosmologies we compare our mass function fits to the Tinker and Warren model with fits to the same models, compiled from Warren et al. (2006, Table 8), Tinker et al. (2010, Table 4, ), Crocce et al. (2010, Table 2) and Watson et al. (2013, Table 2, FoF Uni.). Crocce et al. (2010) and Watson et al. (2013) fit mass functions to the Warren model. Note that Crocce et al. (2010) also used simulations from the MICE simulation suite, with the same cosmology as MICE-GC, but rely on the nested boxes approach to cover a similar mass range, while having a higher resolution in the low mass end than MICE-GC. Tinker et al. (2010) used spherical overdensities to define haloes. A universal behaviour would not only be useful for PBS bias predictions, but also for constraining with galaxy luminosity functions, statistics of the initial density field and various other application (see e.g. White, 2002).

We compare our mass function fits with those from the literature in Fig. 4. We find that the different mass function fits agree at the level in the low mass end, but differ by up to at high masses with a significance of about in terms of error in the measurement. Departures from universality are expected for different cosmologies but can also result from systematic effects, such as the halo mass definition (e.g. Lacey and Cole, 1994; Sheth et al., 2001; Jenkins et al., 2001; White, 2002; Reed et al., 2007; Lukić et al., 2007; Tinker et al., 2008; Crocce et al., 2010; Courtin et al., 2011; More et al., 2011; Bhattacharya et al., 2011; Castorina et al., 2014). Furthermore, the fitting procedure affects the presented comparison as well.

The comparison between the Warren fit from Crocce et al. (2010) and from this work reveals the strong impact of the latter systematic effects on the fit. These two fits agree well in the high mass end, when the fit is performed over the whole mass range M0123. Interestingly we find that these fits differ more strongly from the measurements in the high mass end than fits from other simulations. Excluding lower masses from the fit (M23) leads to a better agreement between our fit and the measurement in the high mass end and therefore to a stronger difference between the results from Crocce et al. (2010) and ours. The lower amplitude of the Crocce et al. (2010) fit at low masses indicates that the low halo mass MICE-GC halo samples include more spuriously linked FoF groups, which can be expected from the low resolution as we concluded before in this section. Furthermore, a lower mass resolution leads to an overestimation of halo masses. Correcting this effect as suggested by Warren et al. (2006) and done by Crocce et al. (2010) results in a decrease of the amplitude, which is shown as grey symbols in Fig. 4. The fact that our Warren corrected mass function is lower than all mass function fits in the low mass range () indicates that the Warren correction leads to an underestimation of halo masses when it is applied on FoF groups with order of particles. For intermediate masses () our Warren corrected measurements are in better agreement with the results from Crocce et al. (2010) than those without Warren correction. A comparison between the Warren corrected MICE-GC mass function at and the Crocce et al. (2010) fit, presented by Crocce et al. (2013), also shows higher amplitude of the prediction compared to the measurement in the highest mass bin at and an opposite trend for lower masses.

## 3 PBS bias predictions

The bias parameters , introduced in equation (1), can be obtained from derivatives of the halo mass function via the PBS approach (Bardeen et al., 1986; Cole and Kaiser, 1989; Mo and White, 1996). Following Scoccimarro et al. (2001) we derive the first-, second- and third-order bias parameters from the mass function fits as

(7) |

(8) |

(9) |

where the parameters and are given by the spherical collapse model. and are computed from the fitted parameters in the mass function models as shown in Table 5. Note that the non-linear bias parameters (equations (8) and (9)), derived from the expressions in Table 5, are here presented for the first time for the Tinker model. Applying the parameter constraints from Table 2 delivers the equivalent expression for the PS, ST and the Warren models, as well as for our proposed model.

Predictions for , and , derived from the Tinker mass function fits at , are shown as a function of FoF halo mass in Fig. 5. The results are based on mass function fits over the whole mass range (M0123) and fits over the higher mass ranges (M123 and M23). The predictions for the different fitting ranges agree in the high mass end where the fitting ranges overlap and the mass function fits agree as well (see Fig. 2). In the low mass end we find a clear, but relatively weak dependence of the linear bias prediction on the fitting range. In the case of and this dependence is stronger and reaches to higher halo masses. This indicates that second- and third-order derivatives of the mass function, used to derive and , cannot be measured as reliable as first-order derivatives, used to derive . We see the same trends when employing the ST and Warren mass function models as well as for our proposed model, while in these cases the dependence on the fitting range is weaker (see Appendix C). We also find a similar behaviour of the bias predictions at .

The absolute deviations between bias prediction from the Tinker mass function, fitted over the range M123 and other predictions are shown in Fig. 6. These other predictions are based on Tinker and Warren fits over different mass ranges and fits for the same models compiled from the literature. We do not show relative deviations to avoid singularities at the zero crossings of and . For the linear bias we find absolute deviations between the different predictions of , which roughly corresponds to relative deviations of around . The relative deviations for and are around , but can go up to more than . Mass function fits over the high mass range M23 to the Tinker and Warren models deliver almost identical bias predictions, which can be expected since also the fitted parameters are very similar (see Fig. 3). In the high mass end these two bias predictions agree with prediction from the fit to the Warren mass function given by Watson et al. (2013). Comparing our results to those of Crocce et al. (2010) we find a reasonable agreement for bias predictions based on the Warren model fitted over the whole mass range M0123.

### 3.1 Universal relation between bias parameters

A universal behaviour of the mass function, as studied in Section 2.2, would suggest that the bias parameters, derived from the mass function are universal as well, when they are expressed as a function of peak height . Our comparison with the literature shows that both, the mass function from different simulations and the bias parameters derived from these mass functions (especially and ) can differ significantly from each other. These disagreements might not only arise from different cosmologies, but also systematic effects, as discussed previously.

We now aim to verify the universality of the relation between the bias parameters. Such a universal behaviour would be useful for reducing uncertainties in linear bias measurements from third-order statistics (e.g. Manera and Gaztañaga, 2011; Hoffmann et al., 2015). In Fig. 7 we show the PBS prediction of the second- and third-order bias parameters, and , as a function of the prediction for the linear bias . We find a agreement for the relations for large values of the linear bias (). These relations appear to be well described by second- and third-order polynomials in the case of and respectively, with

(10) |

as we demonstrate in the same figure. This finding can be expected from expressing and as functions of with the PS model (Table 2). For this model the parameters can be directly predicted as for and for . However, we find smaller rms values with respect to the Tinker and Warren predictions for the and relations, when we leave as free parameters. We show values for from fits to the Tinker predictions in Fig. 7.

For predictions based on our fits over the whole mass range, M0123, we find deviation from this universal behaviour, while these results involve the low mass samples which we found to be unreliable previously, possibly due to low mass resolution and noise in the halo detection. For lower values the different predictions differ more strongly from each other. However, a weakly universal relation, especially between and , might already help to improve constraints from third-order statistics as these two parameters are usually treated as independent. A comparison of the relation, predicted by the PBS with measurements from combined second- and third-order clustering was presented by Saito et al. (2014), who also find that this relation is consistent with redshift independence. We will pursue the study of this relation with different measurements of and in a future analysis (Bel, Hoffmann & Gaztañaga in preparation).

## 4 Bias prediction versus measurements from clustering

In the previous section we found that the PBS bias predictions depend on the employed mass function model and the mass range over which the models are fitted. We now aim to verify how the predictions for the linear bias, , in these different cases compare to linear bias measurements from the two-point halo-matter cross-correlation. A comparison of second- and third-order bias parameter predictions with other measurements will be presented in a future analysis (Bel, Hoffmann & Gaztañaga in preparation).

The two-point cross-correlation between halo- and matter density fields, , can be measured as the mean product of smoothed fluctuations of each density field, , at the positions and as a function of the scale ,

(11) |

The measurements for the four halo mass samples M0-M3 at the redshifts and are shown in the top panels of Fig. 8. The amplitude increases with halo mass as expected from the PBS predictions. The growth of matter fluctuations further contributes to a change with redshift. At around Mpc shows a local maximum which results from baryonic acoustic oscillations in the initial power spectrum of the simulation.

A relation between the two-point halo-matter cross-correlation and the two-point matter auto-correlation, , via the halo bias can be obtained by inserting the local bias model from equation (1) into equation (11),

(12) |

At large scales (Mpc) we expect the higher-order contributions to be negligible, which allows for measurements of the linear bias as

(13) |

The measurements of are shown in the bottom panel of Fig. 8. We fit between Mpc, where the scale-independence is a good approximation. Non-linear terms impact at smaller scales, but also around the scale of baryonic acoustic oscillations. Comparing these bias measurements from the cross-correlation to those from the auto-correlation, shown in Hoffmann et al. (2015), we find that non-linearities have a stronger effect on the autocorrelation. However, differences in the bias from auto- and cross-correlations are small compared to differences between these measurements and the PBS predictions. We present a detailed analysis on the impact of non-linearities on bias from second-order statistics in Bel et al. (2015).

To compare the PBS predictions for the linear bias with the measurements from the two-point correlation we calculate the average bias prediction in each of the mass samples M0-M3, weighted with the halo number density ,

(14) |

and are the lower and upper limits of each mass sample M, given in Table 1. PBS predictions, based on fits to the Tinker model over the mass range M123, are compared with the measurements in the bottom panel of Fig. 8. For the high mass samples M2 and M3 we find clear deviations between measurements and predictions as the latter are significantly too low, on all scales.

The dependence of these deviations on the mass function model and the mass range in which the models are fitted is shown in Fig. 9. Fitting the mass function over the whole mass range, M0123, delivers predictions which tend to be , below the measurements, except for the low mass samples M0 and M1 at , for which we find up to deviations in the opposite direction. We find the strongest variations between bias predictions from different models when, i) the low mass range at is included in the mass function fitting range, or ii) when the bias is predicted for mass samples which are not within the fitting range (e.g. bias predictions for the mass sample M1, based on fits over the mass range M23). The first case i) might be explained by noise, contaminating the FoF halo detection, which results in the poor mass function fits shown in Fig. 2 (see discussion in Section 2.2). In this latter figure we also see that the mass function fits outside the fitting range can strongly differ for different models. This could cause the strong differences in the bias predictions, described above as case ii). We do not see that the deviations between PBS bias predictions and measurements decrease when the analysis is restricted to the higher mass range. This is true for both redshifts ( and ) and consistent with results from Manera and Gaztañaga (2011).

However, restricting the fitting range to the higher mass range M23 we find a good agreement between the linear bias predictions from different mass functions models at and . The fact that the fitting performance strongly differs for the different models (see Fig. 3), while all models predict a linear bias with similar deviations to the measurements, suggests that the goodness of the mass function fit is not the only reason for these deviations, as mentioned in the introduction to this article. These results line up with reports of Manera et al. (2010), who also find the linear PBS bias prediction to lie below measurements from the power spectrum and two-point correlations, especially at high halo masses. As in our case their result is independent of the employed mass function model and the way it is fitted to the measurements.

Furthermore, these authors investigate if differences between the predictions and measurements are related to the mass definition of haloes. They therefore perform their analysis using FoF groups with different linking lengths, as well as spherical overdensities to define halo masses. In both cases they find that the PBS model underpredicts the linear bias measurements. In fact one could expect that the halo mass should be higher than those of FoF groups in order to match the PBS predictions (since shifting the measurements in Fig. 9 to higher masses would decrease the deviations between measurements and predictions). However, halo masses defined by spherical overdensities tend to be below those of FoF groups (Tinker et al., 2008). This should lead to higher measurements of the linear bias for spherical overdensities within a given mass range than corresponding measurements for FoF groups, as found by Tinker et al. (2010). The underprediction of linear bias measurements by the PBS model, which we see for high mass haloes in Fig. 9, is therefore probably a lower bound. The consideration above also suggests that applying the Warren correction on the FoF masses could increase the differences between the PBS bias predictions and the measurements. Hence, these differences might not only be related to the halo mass definition, but also to assumptions of the PBS model, such as spherical collapse or a local bias relation (e.g. Schmidt et al., 2013; Paranjape et al., 2013). The conclusion, that bias predictions are only weakly dependent on the employed mass function model does not hold for the higher-order bias predictions and (see Fig. 5, 6 and 15).

### 4.1 Bias ratios

The degeneracy between the growth of matter fluctuations and the bias of observed galaxy samples is one of the largest uncertainties for constraints of cosmological models derived from large-scale structure observations. With estimations of the typical host halo masses of such galaxy samples the PBS model can be employed to predict the bias of these samples to break the growth-bias degeneracy. Besides the galaxies host halo mass estimation, the inaccuracy of the PBS bias prediction constitutes an additional source of error in this approach. Here we aim to quantify the impact of such inaccuracies on measurement of the linear growth factor. The considered growth measurements are based on the ratio of the correlation functions of galaxy samples at two different redshifts, and , multiplied with the inverse ratio of the bias of these samples (see e.g. Hoffmann et al., 2015). The bias ratio needs to be estimated or predicted, while its uncertainties propagate linearly into the growth measurements.

In Fig. 10 we show the PBS bias ratio predictions for the redshifts and and all combinations of the four halo mass samples M0-M1 at these two redshifts. The predictions are based on fits of the Tinker model to the mass function of the mass range M123, which we found to be reliable at both redshifts previously. We find an overall variation of for the higher mass range M123, while deviations are stronger when the low mass sample M0 at redshifts is taken into account. This variation is stronger than uncertainties expected from the measurements. The strong deviations for the low mass range are expected due to the poor mass function fit including M0 at (see Section 2.2). The error in the bias ratio will propagate into error of the growth factor measurement. This uncertainty is lower than the uncertainties found for growth measurements based on bias ratio estimations from the three-point correlation (see Hoffmann et al., 2015). However, the estimation of the galaxies host-halo mass will introduce additional limitation in breaking the growth-bias degeneracy. Furthermore, the precision of any HOD fitting or mass interpretation from clustering measurements will be affected at similar level.

## 5 Summary and Conclusion

We investigated bias predictions from the PBS model, derived via fits to MICE-GC FoF mass functions. The accuracy of this model was tested by comparing its predictions for the linear bias to direct measurements from two-point correlations. In order to verify how the bias predictions are affected by the goodness of the mass function fit, we study the performance of four mass function models, fitted over different mass ranges at the redshifts and . These fits are based on a new mass function error and covariance estimator.

We show that the models of Press and Schechter (1974), Sheth and Tormen (1999) and Warren et al. (2006) are special cases of the mass function expression suggested by Tinker et al. (2010), as they correspond to certain values of free parameters in the Tinker model (see Table 2). This finding motivated us to propose a new model by fixing a different free parameter. The fitting performance of the mass function models, quantified by the minimum values per number of analysed mass function bins (), shows strong variations among different models and fitting ranges (see Fig. 3). All models match the measurements better when the low mass range is excluded from the analysis. This indicates resolution effects, given that we analyse FoF groups with down to particles. We find that the model of Tinker et al. (2010) shows the best overall performance, which can be expected since it contains the highest number of free parameters. Our proposed model delivers results similar to those from the Tinker model when the whole mass range is analysed, while it has one free parameter less. These two models outperform the model of Warren et al. (2006) for fits over the whole mass range. A restriction to the high mass range ( ) leads to very similar fitting performance of the Warren et al. (2006) and Tinker et al. (2010) models with minimum values close to unity, while our proposal is slightly worse. Fits to the model of Sheth and Tormen (1999) show the most significant deviations to the measurements in all cases. These findings are independent of the mass function binning. We find that the inclusion of the covariance into the analysis substantially increases the minimum values of the best fits and also has an impact on the best fit parameters. However, our PBS bias predictions are only very weakly affected by the mass function covariance, especially when the higher mass range is analysed, where errors are shot-noise dominated.

The results described above can be affected by the way the mass function errors and the covariance between different mass function bins are estimated. We therefore conducted a detailed study of these quantities which is presented in Appendix B. Given the one MICE-GC realisation, we rely on the internal JK error estimator which we compared to theory predictions. The comparison reveals that the JK method is in good agreement with the predicted mass function error only in the shot-noise dominated high mass range (), but overestimates the predictions by up to in the lower mass range, where the errors are dominated by sampling variance. We show that this difference arises because the standard JK estimator assumes a wrong scaling relation between sampling variance and sample volume. By introducing an improved scaling relation, predicted from the linear matter power spectrum, we are able to propose a new mass function error estimator. Deviations between errors of our new estimator and the predictions are less than (see Fig. 11). The advantage of the new estimator with respect to predictions is that it does not rely on a model for halo bias and does not depend on the power spectrum normalisation. This approach to JK error estimations can also be applied to other statistics, such as two-point correlation functions (Hoffmann et al. in preparation).

The presence of non-zero off-diagonal elements in the mass function covariance suggests that a similar covariance can be found in the luminosity function or the stellar mass function, as reported in the literature (e.g. Smith, 2012; Benson, 2014). The latter work demonstrated that the negligence of covariance in the stellar mass function significantly affects parameter constraints in semi-analytic models of galaxy formation. In this case a correct estimation of the error and covariance might be important for a correct interpretation of observations within such models.

Our FoF mass function measurements show no significant () change between the redshifts and for haloes with more than particles (corresponding to a lower halo mass limit of , see Fig. 1). When including lower masses, the redshift dependence is stronger, possibly because of redshift dependent noise in the low mass FoF detection. In order to investigate a dependence on cosmology we compare our results with mass function fits from the literature. We find variations between in the low mass end and up to in the high mass end (see Fig. 4). This finding is in agreement with other studies on departures from the mass function universality (see Section 2.3 for references). The advantage of the MICE-GC simulation is the large volume, which leads to smaller errors in the high mass range and therefore allows for a more significant assessment of the aforementioned variations in the mass function amplitude.

Our analysis demonstrated that numerical and systematic effects, such as halo mass definition, resolution effects or the fitting procedure can contribute to variations in the mass function amplitude with a similar impact as differences in the cosmology (see Section 2.2). This result indicates that understanding such numerical and systematic effects is important for discriminating different cosmologies using the mass function. In fact it was shown in the literature that uncertainties in the mass function in the order of magnitude that we find in our analysis, can strongly affect constraints on the matter density, the dark-energy equation of state, the power spectrum amplitude or neutrino properties from surveys such as DES or Euclid (see e.g. Crocce et al., 2010; Wu et al., 2010; Costanzi Alunno Cerbolini et al., 2013; Weinberg et al., 2013; Appleby et al., 2013; Basse et al., 2014; Bocquet et al., 2015).

After comparing fits from different mass function models to MICE-GC measurements, we study the bias prediction, derived from the best performing model of Tinker et al. (2010). Note that the non-linear bias parameter expressions for this model are presented in this work for the first time. We find that the bias prediction depends on the mass range over which the model was fitted as the amplitude of the linear bias predictions varies by around for different fitting ranges. For the second- and third-order bias parameters the amplitude can vary by more than . These dependences of the bias predictions on the fitting range are comparable with variations obtained when employing fits to other mass function models. Furthermore we find deviations with similar amplitudes in a comparison with bias prediction from mass function fits to other simulations, compiled from the literature (see Fig. 6).

A universal behaviour of the mass function would suggest that the bias parameters, derived from the mass function are universal as well. Despite the strong variation among different bias predictions we find a tight universal relation between and or for across different simulations and mass function models. For smaller values, these relations are more dependent on the mass function fit, but still quite tight. Using the PS mass function model we derive that the second- and third-order bias parameters and can be expressed as second- and third-order polynomials of the linear bias (see Fig. 7). These findings suggests that the linear bias can, at least, constrain the non-linear bias parameters. This could be used to improve the linear bias measurements from third-order statistics.

A common application of the PBS model is to predict the linear bias from clustering. We measured the latter directly from the two-point halo-matter cross-correlation at large scales in the MICE-GC and compare it to the PBS predictions. The comparison was conducted using four different mass samples at the redshifts and . Excluding the low mass sample M1 with less than particles per halo from the analysis, which we expect to be affected by noise, we find that the linear bias, predicted from the PBS, model lies below results from the two-point correlation (see Fig. 9). This effect is similar at the redshifts and and independent of the employed mass function model and the way it is fitted to the measurements, confirming previous findings (Manera et al., 2010; Manera and Gaztañaga, 2011). Including the low mass sample delivers similar results, but with a larger scatter among the models. From the analysis in the higher mass ranges we conclude that shortcomings in the fitting performance of the mass function model are not the main reason for the discrepancy between PBS predictions for the linear bias and the corresponding measurements from clustering. An alternative reason for such discrepancies might be given by the overestimation of halo masses by the FoF algorithm, as those tend to be larger than halo masses of spherical overdensities (Tinker et al., 2008). However, from our results in Fig. 9 we conclude that shifting the linear bias measurements of FoF halo samples to lower masses would increase the deviations between measurements and PBS predictions. Hence, if FoF masses are overestimations of the halo masses described by the PBS model then the differences between linear bias predictions and measurements, found in this analysis, constitute lower bounds for the inaccuracy of the predictions. This indicates that simple assumptions of the PBS model, such as a local bias model or spherical collapse might limit the accuracy of the linear bias predictions. We will present a comparison between the non-linear bias parameters from predictions and measurements in Bel, Hoffmann & Gaztañaga (in preparation).

The deviations between linear bias predictions and measurements will affect at similar level the precision of any HOD fitting or mass interpretation from clustering measurements. We demonstrate the impact of these deviations on growth measurements from two-point correlations. Such measurements are based on the ratio of the linear bias at two different redshifts. Ignoring the unreliable low mass range we find deviations between PBS predictions for the bias ratios and measurements from the two-point correlation. This inaccuracy would propagate linear into measurements of the linear growth factor, based on PBS bias predictions.

## Acknowledgements

Funding for this project was partially provided by the Spanish Ministerio de Ciencia e Innovacion (MICINN), project AYA2009-13936, Consolider-Ingenio CSD2007- 00060, European Commission Marie Curie Initial Training Network CosmoComp (PITN-GA-2009-238356) and research project 2009- SGR-1398 from Generalitat de Catalunya. JB acknowledges useful discussions with Emiliano Sefusatti and support of the European Research Council through the Darklight ERC Advanced Research Grant (#291521). KH is supported by beca FI from Generalitat de Catalunya. 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Ã³ CientÃfica (PIC).

We thank Martin Crocce, Pablo Fosalba, Francisco Castander, Fabien Lacasa and Shun Saito for interesting and useful comments.

## References

- Appleby et al. (2013) Appleby S.A., Linder E.V., Weller J., 2013, Phys.Rev.D, 88, 043526
- Baldauf et al. (2012) Baldauf T., Seljak U., Desjacques V., McDonald P., 2012, Phys.Rev.D, 86, 083540
- Bardeen et al. (1986) Bardeen J.M., Bond J.R., Kaiser N., Szalay A.S., 1986, ApJ, 304, 15
- Basse et al. (2014) Basse T., Eggers Bjælde O., Hamann J., Hannestad S., Wong Y.Y.Y., 2014, Journal of Cosmology and Astroparticle Physics, 5, 021
- Bel et al. (2015) Bel J., Hoffmann K., Gaztañaga E., 2015, ArXiv e-prints
- Benson (2014) Benson A.J., 2014, MNRAS, 444, 2599
- Bhattacharya et al. (2011) Bhattacharya S., Heitmann K., White M., Lukić Z., Wagner C., Habib S., 2011, ApJ, 732, 122
- Biagetti et al. (2014) Biagetti M., Desjacques V., Kehagias A., Riotto A., 2014, Phys.Rev.D, 90, 045022
- Bocquet et al. (2015) Bocquet S., et al., 2015, ApJ, 799, 214
- Bond et al. (1991) Bond J.R., Cole S., Efstathiou G., Kaiser N., 1991, ApJ, 379, 440
- Carretero et al. (2015) Carretero J., Castander F.J., Gaztañaga E., Crocce M., Fosalba P., 2015, MNRAS, 447, 646
- Castorina et al. (2014) Castorina E., Sefusatti E., Sheth R.K., Villaescusa-Navarro F., Viel M., 2014, Journal of Cosmology and Astroparticle Physics, 2, 049
- Chan et al. (2012) Chan K.C., Scoccimarro R., Sheth R.K., 2012, Phys.Rev.D, 85, 083509
- Cole and Kaiser (1989) Cole S., Kaiser N., 1989, MNRAS, 237, 1127
- Cooray and Sheth (2002) Cooray A., Sheth R., 2002, Phys. Rept., 372, 1
- Costanzi Alunno Cerbolini et al. (2013) Costanzi Alunno Cerbolini M., Sartoris B., Xia J.Q., Biviano A., Borgani S., Viel M., 2013, Journal of Cosmology and Astroparticle Physics, 6, 020
- Coupon et al. (2012) Coupon J., et al., 2012, A&A, 542, A5
- Courtin et al. (2011) Courtin J., Rasera Y., Alimi J.M., Corasaniti P.S., Boucher V., Füzfa A., 2011, MNRAS, 410, 1911
- Crocce et al. (2010) Crocce M., Fosalba P., Castander F.J., Gaztañaga E., 2010, MNRAS, 403, 1353
- Crocce et al. (2013) Crocce M., Castander F.J., Gaztanaga E., Fosalba P., Carretero J., 2013, ArXiv e-prints
- Davis et al. (1985) Davis M., Efstathiou G., Frenk C.S., White S.D.M., 1985, ApJ, 292, 371
- Desjacques et al. (2010) Desjacques V., Crocce M., Scoccimarro R., Sheth R.K., 2010, Phys.Rev.D, 82, 103529
- Eriksen and Gaztañaga (2015) Eriksen M., Gaztañaga E., 2015, ArXiv e-prints
- Fosalba et al. (2015a) Fosalba P., Crocce M., Gaztañaga E., Castander F.J., 2015a, MNRAS, 448, 2987
- Fosalba et al. (2015b) Fosalba P., Gaztañaga E., Castander F.J., Crocce M., 2015b, MNRAS, 447, 1319
- Fry and Gaztanaga (1993) Fry J.N., Gaztanaga E., 1993, ApJ, 413, 447
- Hoffmann et al. (2015) Hoffmann K., Bel J., Gaztañaga E., Crocce M., Fosalba P., Castander F.J., 2015, MNRAS, 447, 1724
- Hu and Kravtsov (2003) Hu W., Kravtsov A.V., 2003, ApJ, 584, 702
- Jenkins et al. (2001) Jenkins A., Frenk C.S., White S.D.M., Colberg J.M., Cole S., Evrard A.E., Couchman H.M.P., Yoshida N., 2001, MNRAS, 321, 372
- Lacey and Cole (1994) Lacey C., Cole S., 1994, MNRAS, 271, 676
- Lima and Hu (2004) Lima M., Hu W., 2004, Phys.Rev.D, 70, 043504
- Lima and Hu (2005) Lima M., Hu W., 2005, Phys.Rev.D, 72, 043006
- Lukić et al. (2007) Lukić Z., Heitmann K., Habib S., Bashinsky S., Ricker P.M., 2007, ApJ, 671, 1160
- Manera and Gaztañaga (2011) Manera M., Gaztañaga E., 2011, MNRAS, 415, 383
- Manera et al. (2010) Manera M., Sheth R.K., Scoccimarro R., 2010, MNRAS, 402, 589
- Martino and Sheth (2009) Martino M.C., Sheth R.K., 2009, MNRAS, 394, 2109
- Mo and White (1996) Mo H.J., White S.D.M., 1996, MNRAS, 282, 347
- Mo et al. (1997) Mo H.J., Jing Y.P., White S.D.M., 1997, MNRAS, 284, 189
- More et al. (2011) More S., Kravtsov A.V., Dalal N., Gottlöber S., 2011, ApJS, 195, 4
- Norberg et al. (2009) Norberg P., Baugh C.M., Gaztañaga E., Croton D.J., 2009, MNRAS, 396, 19
- Paranjape et al. (2013) Paranjape A., Sheth R.K., Desjacques V., 2013, MNRAS, 431, 1503
- Press and Schechter (1974) Press W.H., Schechter P., 1974, ApJ, 187, 425
- Pujol and Gaztañaga (2014) Pujol A., Gaztañaga E., 2014, MNRAS, 442, 1930
- Reed et al. (2007) Reed D.S., Bower R., Frenk C.S., Jenkins A., Theuns T., 2007, MNRAS, 374, 2
- Robertson (2010) Robertson B.E., 2010, ApJ, 713, 1266
- Saito et al. (2014) Saito S., Baldauf T., Vlah Z., Seljak U., Okumura T., McDonald P., 2014, Phys.Rev.D, 90, 123522
- Schmidt et al. (2013) Schmidt F., Jeong D., Desjacques V., 2013, Phys.Rev.D, 88, 023515
- Scoccimarro et al. (2001) Scoccimarro R., Sheth R.K., Hui L., Jain B., 2001, ApJ, 546, 20
- Sheth and Tormen (1999) Sheth R.K., Tormen G., 1999, MNRAS, 308, 119
- Sheth et al. (2001) Sheth R.K., Mo H.J., Tormen G., 2001, MNRAS, 323, 1
- Smith (2012) Smith R.E., 2012, MNRAS, 426, 531
- Smith and Marian (2011) Smith R.E., Marian L., 2011, MNRAS, 418, 729
- Springel (2005) Springel V., 2005, MNRAS, 364, 1105
- Tinker et al. (2008) Tinker J., Kravtsov A.V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D.E., 2008, ApJ, 688, 709
- Tinker et al. (2010) Tinker J.L., Robertson B.E., Kravtsov A.V., Klypin A., Warren M.S., Yepes G., Gottlöber S., 2010, ApJ, 724, 878
- Valageas et al. (2011) Valageas P., Clerc N., Pacaud F., Pierre M., 2011, A&A, 536, A95
- Warren et al. (2006) Warren M.S., Abazajian K., Holz D.E., Teodoro L., 2006, ApJ, 646, 881
- Watson et al. (2013) Watson W.A., Iliev I.T., D’Aloisio A., Knebe A., Shapiro P.R., Yepes G., 2013, MNRAS, 433, 1230
- Weinberg et al. (2013) Weinberg D.H., Mortonson M.J., Eisenstein D.J., Hirata C., Riess A.G., Rozo E., 2013, Phys. Rept., 530, 87
- White (2002) White M., 2002, ApJS, 143, 241
- Wu et al. (2010) Wu H.Y., Zentner A.R., Wechsler R.H., 2010, ApJ, 713, 856

## Appendix A Mass function measurements

The mass function measurements are based on a rewritten form of equation 2:

(15) |

where is the mean halo mass in each logarithmic mass bin. If the mass bins are chosen to be exactly equal in logarithmic space, the mass function amplitude slightly oscillates in the low mass range due to mass resolution effects. Since the errors are smallest in the low mass range this artefact can significantly affect the fits, causing a strong dependence of the fits on the number of mass function bins. Aiming to minimise this mass discreteness effect we determine the minimum and maximum number of particles per halo in each logarithmic mass bin, and respectively. The width of the mass bin in then recalculated as , where is the particle mass. The value of of each bin is calculated from the mean mass of haloes in the bin. The term in equation 15 is derived directly from equation(4).

## Appendix B Covariance

In order to fit the mass function we estimate the errors and the covariance between different mass bins. A direct measurement of these quantities would require a large set of independent realisations of the simulation. Since just one realisation of the MICE-GC simulation was run we estimate the errors and the covariance using the Jack-Knife (hereafter referred to as JK) sampling technique. To validate these estimations we compare the results to theoretical predictions, which we will describe first.

### b.1 Covariance prediction

Following Crocce et al. (2010) we derive the covariance prediction for the comoving halo number density from the linear bias relation at large scales,

(16) |

where is the number density of haloes with mass in a volume (in our case the simulation volume) around position , is the matter density contrast in the same volume and is the linear halo bias factor (as before refers to the matter density field when it appears as lower index and to the halo mass when it is used as variable). The last term, , corresponds to noise. We will assume to be Poisson shot-noise and therefore independent of . The predictions for the unconditional mass function can be related to those for the halo number density via equation (2). For the sake of simplicity the following considerations are based on the latter. The covariance matrix for number densities of haloes in the mass bins and is defined as

(17) |

where and denotes the average over statistically independent volumes (note that the introduced here is not related to the used in equation (6) for calculating the values of mass function fits). Inserting the expression for the number density of haloes in mass bin from equation (16) leads to

(18) |

The variance of matter fluctuations can be derived from the power spectrum via equation (4), while is the total mass within the volume in which the mass function is measured (in our case the total mass in the simulation). Since this mass corresponds to a very large smoothing radius we can compute from the linear power spectrum. The sampling variance contribution to the covariance is therefore given by

(19) |

If the noise term is Poissonian it averages out when taking the mean over many independent volumes. The contribution of shot-noise to the covariance is then given by

(20) |

while here is the Kronecker delta. Based on these considerations we can write the total covariance as

(21) |

A more formal derivation for this relation is given by Smith and Marian (2011), see also Robertson (2010); Valageas et al. (2011); Smith (2012). The diagonal elements of the covariance matrix correspond to the predictions for the mass function variance,

(22) |

as given by Crocce et al. (2010). For fitting the mass function we work with the normalised covariance ) and differences normalised to . (note that here refers to the variance of the mass function in the mass bin and not to the variance of the matter field, ).

### b.2 Jack-Knife estimation of covariance

For mass function fits in observations the covariance prediction is of limited use since it requires knowledge about the bias and the power spectrum in advance. This problem might be solved with an iterative approach for the fit, starting from an initial guess for the power spectrum and the linear bias factor. Another possibility to obtain the covariance without knowledge of the bias and the power spectrum is to estimate it with the JK sampling technique. Testing this approach we construct JK samples by subtracting cubical sub-volumes (hereafter referred to as JK cells) with the size from the total simulation volume . The basic assumption of the JK approach is that the error scales with the size of the subtracted volume (e.g. Norberg et al., 2009). We follow the common approach by rescaling the covariance with the factor , which leads to

(23) |

Again , but now is the average over the different JK samples , is the comoving number density of haloes in the mass bin in each JK sample and is the corresponding halo number density in the whole simulation volume. Note that the rescaling factor, , is only weakly justified and can be improved, as we show in Section B.4. As in the case of the predictions the diagonal elements of are the JK estimation for the variance () and we normalise . Note that we can use the JK approach for studying the covariance between low and high mass bins because of the large mass range of the MICE-GC simulation. This analysis would not be possible using nested boxes, where the different mass ranges are covered by different realisations with different box sizes (e.g. Warren et al., 2006; Crocce et al., 2010; Tinker et al., 2010).

### b.3 Covariance prediction versus Jack-Knife estimation

A comparison between the error prediction from equation (21) and the corresponding JK estimation from equation equation (23) (with ) is shown in Fig. 11 for the redshift . The error predictions are based on linear bias predictions from mass function fits to the Tinker model over the whole mass range, for which we expect uncertainties of around (see Section 3). From the prediction we expect the error to be dominated by sampling variance in the low mass end and by shot-noise in the high mass end. At halo masses of both sources are predicted to contribute equally to the total error. The JK error estimation is in good agreement with the predictions in the high mass end (). This indicates that the JK method is working well for different JK cell volumes when the error is dominated by shot-noise. Furthermore, the shot-noise is well described by a Poisson distribution. At halo masses lower than we find the JK error to be up to higher than the prediction.

This overestimation is consistent with results reported by Crocce et al. (2010) using the same simulation box size as the MICE-GC, while for smaller simulation boxes they find the JK error to be lower than the prediction. The fact that the overestimation of the JK error in the low mass end is larger for smaller JK cells indicates that the JK assumption of a linear relation between errors and volume is inadequate when sampling variance is the dominating source for error. However, increasing the size of the JK cells results in a smaller number of samples and therefore a stronger noise on the estimated error. In Fig. 11 we also show a new JK error estimation, which is in good agreement with the predictions at all masses. This new JK error is based on an improved scaling between the sampling variance in a JK cell and in the whole simulation box using the linear power spectrum, as explained in Section B.4.

In Fig. 12 we compare the normalised covariance of the mass function between the mass bins and , predicted via equation (21) with the JK estimation from equation (23) using JK samples at . The shape of the covariance is in good agreement with results from Smith and Marian (2011). The low mass bins are highly correlated because of sampling variance, while high mass bins are uncorrelated as their errors are dominated by shot-noise. For the comparison of the variances we find a reasonable agreement between the prediction and the JK estimations, especially in the high mass end. In the low mass end the covariance seems to be overestimated by the JK approach, while the new JK method reproduces the prediction well.

We show a more detailed comparison of the covariance amplitudes in Fig. 13, fixing one mass bin and varying the second mass bin . For JK cells we find the normalised JK covariance amplitudes to be higher than the predictions with differences of up to . Using larger JK cells this overestimation slightly decreases, while results become more noisy. Again the improved estimation is in better agreement with the prediction. We have verified that our conclusions also hold for redshift .

### b.4 Improved JK estimator

We now aim to understand the disagreement between the predicted mass function error and the corresponding JK estimation in order to improve the latter. The JK samples are constructed by subtracting haloes in JK cells of the size from the total halo distribution. The number of haloes in the remaining JK sample is then given by . The volume of a JK sample is given by . From the definition of the number density, (), and the deviation from the mean over the total volume one can derive . Note that this relation also holds for the number density contrast, . Hence, subtracting an overdense JK cell from the total volume generates a slightly underdense JK sample. This result leads to a relation between the variances of the number density, , in the JK cells, and the corresponding variance for the JK samples,

(24) |

As for equation (17) denotes the ensemble average. The variance of the JK samples is therefore simply related to the variance at the scales of the JK cells. Note that is not the variance at the scale of the JK sample volume, , since the JK samples are not independent from each other.

From the linear bias model we assume that the variance of the halo number density results from shot-noise () and sampling variance (), as explained in Section B.2. The latter contribution to the total variance of the JK cells, measured via equation (24), is therefore given by

(25) |

Since , the shot-noise for JK cells is related to the shot-noise of the whole box as . To obtain the sampling variance at the scale of the simulation box, , we multiply with a rescaling factor

(26) |

which can be predicted from the linear matter power spectrum. This prediction is based on the assumption that, at large smoothing scales, the sampling variance of the halo number density is related to the dark matter variance by the linear bias factor, . Since is constant at large scales (see Fig. 8), it cancels out in the rescaling factor, hence . The prediction is then based on , computed from the linear matter power spectrum via equation (4). We can now write the expression for the sampling variance of the simulation box, based on equation (25) in the general case of the covariance

(27) |

Note that we have now dropped the index in for simplicity and to be consistent with equation (23) for the standard JK estimator. As in the latter equation the lower indices refer to the mass bins and . With the Poisson shot-noise term from equation (20), the total covariance is then given by equation (21) as . The resulting expression constitutes a new error estimation for the mass function which combines direct measurements of sampling variance via the JK sampling with predictions for the rescaling factor and the shot-noise. This new estimator can be written more explicitly as

(28) |

As before the diagonal elements correspond to the variance, . Note that for Poisson shot-noise dominated errors the sampling variance can be approximated as and the new estimator reduces to the shot-noise term, . In this case we derive from equation (27) that . For large numbers of JK samples () this expression is equivalent to . The left hand side of this relation is the standard JK estimator. This consideration explains the good agreement between standard JK estimator with the improved JK estimator and the predictions at high masses, where the errors and the covariance are shot-noise dominated (Fig. 11, 12 and 13). In the low mass range our new method is in much better agreement with the predictions than the standard JK error estimator (23). This can be understood with the following consideration. For large numbers of JK samples () the new estimator corresponds to the standard JK estimator if . Since , this condition is equivalent to . The JK approach can therefore be described as the assumption that for large . We show , computed from the linear power spectrum via equation (4) in Fig. 14. The JK assumption is in a clear disagreement with the prediction which causes a too high and therefore an overestimation of sampling variances at the scale of the simulation box for the standard JK assumption.

The advantage of the new JK estimation with respect to the prediction is that it does not require knowledge of the halo bias. Furthermore, this new approach is independent of the normalisation of the power spectrum as it cancels out in the rescaling factor (equation (26)). However, the large scale power spectrum still needs to be known for accurate rescaling of the sampling variance via . For simulations the linear power spectrum is given. In this case the new method can be used instead of running several realisations for deriving mass function errors and covariances. In observations the large scale power spectrum can only be assumed. However, with such an assumption the accuracy of the error estimation can still be improved with respect to the standard JK method, which also implies the strong assumption of . The advantage of the new JK estimation with respect to using independent subvolumes for the error estimation is that the JK samples cover larger volumes with larger average numbers of massive haloes. The covariance between the low- and high mass end of the mass function is therefore better sampled by JK samples than subvolumes. We employ our new method for the error and covariance estimation using samples.

## Appendix C PBS bias predictions

We show in Fig. 15 the PBS bias predictions based on the mass function models studied in this analysis. The different predictions are based on fits over the four mass ranges M0123, M123, M23 and M012, defined in Table 3. The figure is analogous to Figure 5. We find that the linear bias parameter is less sensitive to the mass function model and the mass function fitting range than the non-linear bias parameters á and . I addition to Figure 5 we show that the bias predictions become unstable when the low mass sample, M1, is included in the analysis. Furthermore we show bias predictions, based on mass function fits over the range M123, which were derived neglecting the covariance between different mass function bins. We find that the for this example the mass function covariance has a smaller impact on the bias prediction than the choice of the mass function model, or the mass function fitting range. We expect the impact of the mass function covariance to increase, when low mass samples are included in the fit. However, the low mass range is hard to access for analysis of halo abundance in the MICE-GC, due to resolution effects (see Section 2.2).