# The universality of the virial halo mass function and models for non-universality of other halo definitions

Giulia Despali, Carlo Giocoli, Raul E. Angulo, Giuseppe Tormen, Ravi K.

Sheth, Giacomo Baso, Lauro Moscardini

► To cite this version:

Giulia Despali, Carlo Giocoli, Raul E. Angulo, Giuseppe Tormen, Ravi K. Sheth, et al.. The universality of the virial halo mass function and models for non-universality of other halo definitions. Monthly Notices of the Royal Astronomical Society, 2016, 456 (3), pp.2486-2504. 10.1093/mnras/stv2842. hal-01440707

HAL Id: hal-01440707

https://hal.science/hal-01440707

Submitted on 9 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The universality of the virial halo mass function and models for non-universality of other halo definitions

Giulia Despali,<fnref n="1" />,<fnref n="2" />★ Carlo Giocoli,<fnref n="1" /> Raul E. Angulo,<fnref n="3" /> Giuseppe Tormen,<fnref n="2" /> Ravi K. Sheth,<fnref n="4" /> Giacomo Baso<fnref n="2" /> and Lauro Moscardini<fnref n="5" />,<fnref n="6" />,<fnref n="7" />

Accepted 2015 December 2. Received 2015 October 31; in original form 2015 July 20

## ABSTRACT

The abundance of galaxy clusters can constrain both the geometry and growth of structure in our Universe. However, this probe could be significantly complicated by recent claims of non-universality-non-trivial dependences with respect to the cosmological model and redshift. In this work, we analyse the dependence of the mass function on the way haloes are identified and establish if this can cause departures from universality. In order to explore this dependence, we use a set of different N-body cosmological simulations (Le SBARBINE simulations), with the latest cosmological parameters from the Planck collaboration; this first suite of simulations is followed by a lower resolution set, carried out with different cosmological parameters. We identify dark matter haloes using a spherical overdensity algorithm with varying overdensity thresholds (virial, 2000, 1000, 500, 200 $\rho_c$ and 200 $\rho_b$) at all redshifts. We notice that, when expressed in terms of the rescaled variable $v$, the mass function for virial haloes is a nearly universal as a function of redshift and cosmology, while this is clearly not the case for the other overdensities we considered. We provide fitting functions for the halo mass function parameters as a function of overdensity, that allow us to predict, to within a few per cent accuracy, the halo mass function for a wide range of halo definitions, redshifts and cosmological models. We then show how the departures from universality associated with other halo definitions can be derived by combining the universality of the virial definition with the expected shape of the density profile of haloes.

**Key words:** methods: numerical – galaxies: haloes – cosmology: theory – dark matter.

## 1 INTRODUCTION

In the cold dark matter (CDM) model, structures – up to protogalactic scales – form through the amplification of small density fluctuations via gravitational instability (Longair 1998; Springel et al. 2005; Mo, van den Bosch & White 2010; Angulo et al. 2012). Dark matter haloes are objects which have been able to break away from the expanding background, and collapse (Press & Schechter 1974; Springel et al. 2001). Small haloes form first, before merging with one another to form ever more massive ones in a hierarchical

process (Lacey & Cole 1993, 1994). As a result of repeated mergers, dark matter haloes grow more massive in time (Tormen, Moscardini & Yoshida 2004). Haloes hosting galaxy clusters represent the most massive and recently formed structures in our Universe. The formation and the merger rates of dark matter haloes are sensitive to the expansion history of the Universe and so can be used to constrain cosmological parameters (Lacey & Cole 1993, 1994; Moreno, Giocoli & Sheth 2008).

In particular, different theoretical studies have shown that the halo mass function and its evolution are important probes of the very early Universe, its expansion history, and the nature of gravity (Press & Schechter 1974; Bond et al. 1991; Lacey & Cole 1993; Sheth, Mo & Tormen 2001). These have shown that, in

<footnote n="1">¹Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388 Marseille, France</footnote>

<footnote n="2">²Dipartimento di Fisica e Astronomia, Università degli Studi di Padova, vicolo dell'Osservatorio 3, I-35122 Padova, Italy</footnote>

<footnote n="3">³Centro de Estudios de Física del Cosmos de Aragón (CEFCA), Plaza San Juan 1, Planta-2, E-44001 Teruel, Spain</footnote>

<footnote n="4">⁴Center for Particle Cosmology, University of Pennsylvania, 209 South 33rd St, Philadelphia, PA 19104, USA</footnote>

<footnote n="5">⁵Dipartimento di Fisica e Astronomia, Alma Mater Studiorum Università di Bologna, Viale Berti Pichat 6/2, I-40127 Bologna, Italy</footnote>

<footnote n="6">⁶INAF, Osservatorio Astronomico di Bologna, via Ranzani 1, I-40127 Bologna, Italy</footnote>

<footnote n="7">⁷INFN, Sezione di Bologna, viale Berti Pichat 6/2, I-40127 Bologna, Italy</footnote>

<footnote>* E-mail: gdespali@gmail.com</footnote>
appropriately scaled units, the mass function can be written in an approximately universal form which is independent of power spectrum and expansion history. Although this universality is only expected to be approximate (Musso & Sheth 2012; Paranjape, Sheth & Desjacques 2013), it vastly simplifies the process of constraining cosmological parameters from observational data sets, so it has served as the basis for fitting functions whose parameters are calibrated using numerical simulations of the dark matter (Del Popolo & Gambera 1998, 1999; Sheth & Tormen 1999, 2002; Jenkins et al. 2001; Warren et al. 2006). As simulated data sets have grown, it has become possible to quantify (small) departures from universality (Tinker et al. 2008; Crocce et al. 2010; Manera, Sheth & Scocci-marro 2010; Wu, Zentner & Wechsler 2010; Courtin et al. 2011; Corasaniti & Achitouv 2011; Murray, Power & Robotham 2013; Watson et al. 2013). Several recent works have also concentrated on the small but significant modifications by neutrinos (Castorina et al. 2014), coupling between dark matter–dark energy (Cui, Baldi & Borgani 2012; Giocoli et al. 2013) and different baryon physics (Cui, Borgani & Murante 2014; Bocquet et al. 2015).

Many cosmological constraints have been obtained from cluster counts (Vikhlinin et al. 2009; Rozo et al. 2010; Planck Collaboration XX 2014) and other galaxy cluster properties (Evrard et al. 2008; Ettori et al. 2009; Giocoli et al. 2012). For an extensive review, see Borgani & Kravtsov (2011). However in the near future, many wide-field surveys are expected to use the cluster mass function to constrain cosmological parameters (Laureijs et al. 2011; Pillepich, Porciani & Reiprich 2012; Boldrin et al. 2015; Sartoris et al. 2015). In the light of these, a mass function calibrated to an accuracy of a few per cent, and flexible enough to account for different halo identification definitions, is of primary importance.

Departures of the halo mass function from universality may depend on how haloes are defined. One of the main goals of this paper is to explore this dependence. Our major result is that, if haloes are defined using the virial density, and the fitting formula includes a parameter which is related to this – as was done by Sheth & Tormen (1999) – then the mass function can be considered universal to within a few per cent. In this respect, our findings confirm those of Courtin et al. (2011): departures from universality result from ignoring the redshift and cosmology dependence of these quantities. Moreover, the departures from universality associated with other halo definitions can be derived from combining the universality of the virial definition with knowledge of the enclosed density profile around haloes.

The paper is organized as follows: In Section 2, we describe the cosmological simulations we use for our study. It also describes our halo finder. We present our reference model for the halo mass function in Section 3. In Sections 4 and 5, we discuss the universality associated with using the virial overdensity to define haloes and show how the parameters of the halo mass function depend on how haloes are defined. The bulk of this analysis is for spherical haloes: Appendix A shows how our results are modified if haloes are allowed to be ellipsoidal, and Appendix B describes how a number of more technical details affect our measurements. In Section 6, we show how the non-virial halo mass functions can be derived from combining knowledge of the dark matter density profile with the virial mass function. In Section 7, we present some comparisons with previous works. We discuss our results and conclude in Section 8. Our analysis suggests that the most massive end of the halo mass function is particularly simple, in the sense that it can be described by a function with fewer free parameters. Appendix C describes how this impacts cluster cosmology. All logarithms where not explicitly stated in the text are in base ten.

# 2 THE NUMERICAL SIMULATIONS

## 2.1 Le SBARBINE simulations

Le SBARBINE simulations are a set of six dark-matter-only cosmological simulations run by the Padova cosmology group. These simulations follow the evolution of $1024^3$ particles, whose motions are assumed to be driven by gravitational instability, using the publicly available code GADGET-2 (Springel 2005). The assumed background cosmology and initial conditions for these runs are consistent with recent Planck results (Planck Collaboration XVI 2014, hereafter Planck13). In particular, we have set: $\Omega_m = 0.307$, $\Omega_\Lambda = 0.693$, $\sigma_8 = 0.829$ and $H = 100\,h\,\text{km s}^{-1} h^{-1}\,\text{Mpc}$ with $h = 0.677$. The initial power spectrum was generated using the CAMB code (Lewis, Challinor & Lasenby 2000), and initial conditions were produced by perturbing a glass distribution with N-GenIC (http://www.mpa-garching.mpg.de/gadget); the realizations have been carefully chosen, in order to follow the initial power spectrum even at large scales and thus to reduce the differences between the linear spectrum and that measured from the simulations.

The parameters of our main simulation set are listed in Table 1. We used a different seed for the random number generator which sets the initial conditions of each simulation, so as to have a sample of independent realizations. Although each box contains the same number of dark matter particles ($1024^3$), the comoving box lengths are different, so the mass resolution in each box is different. The box sizes were chosen so that the set provides good mass resolution down to $10^7 h^{-1}\,\text{Mpc}$.

Le SBARBINE simulations are complemented by a set of lower resolution runs having different cosmological parameters. These all have $512^3$ dark matter particles. In particular, for each set of cosmological parameters we ran two simulations: one with box size $150\,h^{-1}\,\text{Mpc}$ and another with $1000\,h^{-1}\,\text{Mpc}$. These were chosen to ensure good resolution both for intermediate and high-mass haloes. This lower resolution set was produced specifically to test the universality of the halo mass function with respect to the cosmological model (see Section 4.2). The parameters of these other simulations are listed in Table 2.

In addition, we also re-ran three simulations for which the initial conditions were generated using a second order Lagrangian Perturbation Theory (2LPT) algorithm<fnref n="1" /> (Crocce, Pueblas & Scocci-marro 2006). These are copies of Bice ($1024^3$ particles) and of the two $512^3$ simulations with the WMAP7 cosmology (Komatsu et al. 2011) namely (wmap7 and wmap7-big). As we discuss in Appendix B4, there are small differences – not exceeding 5 per cent – between the mass functions in simulations with Zel'dovich versus 2LPT initial conditions. But these differences appear only for high $v$ and high redshifts, which play little role in our calibration of the mass function parameters.

All runs were performed in Padova on ‘Nemo’: a SuperServer Twin 2U Dual Xeon Sandy Bridge composed by four independent node servers each equipped with two Xeon Sandy Bridge 8 Core E5-2670 and 128 GB of RAM, for a total of 64 cores or 128 CPU-threads and 512 GB of RAM.

## 2.2 Halo catalogues

For each stored particle snapshot, we identify haloes using a spherical overdensity (SO) algorithm (e.g. Tormen 1998; Tormen et al.

<footnote n="1">¹ http://cosmo.nyu.edu/roman/2LPT</footnote>
<table><caption>**Table 1.** Features of Le SBARBINE simulations run with Planck13 parameters $\Omega_m = 0.307$, $\Omega_\Lambda = 0.693$, $\sigma_8 = 0.829$ and $h = 0.677$ and containing 1024³ dark matter particles. The last two columns report the total number of haloes identified with the spherical overdensity at redshift $z = 0$ that are resolved with more than 10 and 300 particles, respectively.</caption><thead><tr><th rowspan="2">Name</th><th rowspan="2">Box (h<sup>-1</sup> Mpc)</th><th rowspan="2">z<sub>i</sub></th><th colspan="2">Main set of simulations</th><th rowspan="2">N<sub>h-tot</sub>(z = 0)</th><th rowspan="2">N<sub>h&gt;300</sub>(z = 0)</th><th rowspan="2">Colour</th></tr><tr><th>m<sub>p</sub>(M<sub>☉</sub> h<sup>-1</sup>)</th><th>Soft (kpc h<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Ada</td><td>62.5</td><td>124</td><td>1.94 × 10<sup>7</sup></td><td>1.5</td><td>2264 847</td><td>103 852</td><td>Green</td></tr><tr><td>Bice</td><td>125</td><td>99</td><td>1.55 × 10<sup>8</sup></td><td>3</td><td>2750 411</td><td>129 674</td><td>Cyan</td></tr><tr><td>Cloe</td><td>250</td><td>99</td><td>1.24 × 10<sup>9</sup></td><td>6</td><td>3300 880</td><td>161 580</td><td>Blue</td></tr><tr><td>Dora</td><td>500</td><td>99</td><td>9.92 × 10<sup>9</sup></td><td>12</td><td>3997 898</td><td>191 793</td><td>Magenta</td></tr><tr><td>Emma</td><td>1000</td><td>99</td><td>7.94 × 10<sup>10</sup></td><td>24</td><td>4739 379</td><td>176 633</td><td>Red</td></tr><tr><td>Flora</td><td>2000</td><td>99</td><td>6.35 × 10<sup>11</sup></td><td>48</td><td>5046 663</td><td>75 513</td><td>Orange</td></tr></tbody></table>

<table><caption>**Table 2.** Details of the small set of 10 simulations with different cosmological parameters. Each contains 512³ dark matter particles with initial conditions generated at redshift $z = 99$. For all the models the Hubble parameter is $h = 0.6777$, apart from the WMAP7 cosmology for which $h = 0.704$.</caption><thead><tr><th rowspan="2">Name</th><th colspan="4">Secondary set of simulations</th><th rowspan="2">m<sub>p</sub>(M<sub>☉</sub> h<sup>-1</sup>)</th><th rowspan="2">Colour</th></tr><tr><th>Ω<sub>m</sub></th><th>Ω<sub>Λ</sub></th><th>σ<sub>8</sub></th><th>box (h<sup>-1</sup> Mpc)</th></tr></thead><tbody><tr><td>Tea</td><td>0.2</td><td>0.8</td><td>0.7</td><td>150</td><td>1.396 × 10<sup>9</sup></td><td>Grey-square</td></tr><tr><td>Tea-big</td><td>0.2</td><td>0.8</td><td>0.7</td><td>1000</td><td>4.135 × 10<sup>11</sup></td><td>Grey-square</td></tr><tr><td>Tina</td><td>0.2</td><td>0.8</td><td>0.9</td><td>150</td><td>1.396 × 10<sup>9</sup></td><td>Grey-triangle</td></tr><tr><td>Tina-big</td><td>0.2</td><td>0.8</td><td>0.9</td><td>1000</td><td>4.135 × 10<sup>11</sup></td><td>Grey-triangle</td></tr><tr><td>Vera</td><td>0.4</td><td>0.6</td><td>0.7</td><td>150</td><td>2.791 × 10<sup>9</sup></td><td>Brown-square</td></tr><tr><td>Vera-big</td><td>0.4</td><td>0.6</td><td>0.7</td><td>1000</td><td>8.271 × 10<sup>11</sup></td><td>Brown-square</td></tr><tr><td>Viola</td><td>0.4</td><td>0.6</td><td>0.9</td><td>150</td><td>2.791 × 10<sup>9</sup></td><td>Brown-triangle</td></tr><tr><td>Viola-big</td><td>0.4</td><td>0.6</td><td>0.9</td><td>1000</td><td>8.271 × 10<sup>11</sup></td><td>Brown-triangle</td></tr><tr><td>Wanda (wmap7)</td><td>0.272</td><td>0.728</td><td>0.81</td><td>150</td><td>1.898 × 10<sup>9</sup></td><td>Blue-circle</td></tr><tr><td>Wanda-big (wmap7)</td><td>0.272</td><td>0.728</td><td>0.81</td><td>1000</td><td>5.624 × 10<sup>11</sup></td><td>Blue-circle</td></tr></tbody></table>

2004; Giocoli, Pieri & Tormen 2008). We chose this rather than the Friends-of-Friends (FoF) method of Davis et al. (1985), because we believe it to be slightly closer to physical models of halo formation, and because it is quite similar to how mass is defined in observational data.

For each particle distribution, we estimate the local dark matter density at the position of each particle by calculating the distance $d_{i,10}$ to the tenth nearest neighbour. In this way, we assign to each particle a local density $\rho_i \propto d_{i,10}^{-3}$. We then sort particles in density and choose as centre of the first halo the position of the densest particle. We grow a sphere around this centre, and stop when the mean density within the sphere falls below a desired critical value. At this point, we assign all particles within the sphere to the newly identified halo, and remove them from the global list of particles. We then choose the densest particle of the ones remaining, and repeat (i.e. grow a sphere around it until the mean enclosed density falls below threshold, etc.). We continue in this manner until none of the remaining particles has a local density large enough to be the centre of a 10 particle halo (as we discuss shortly, we apply a more stringent cut when we fit for the mass function); particles not assigned to any halo are called 'field' or 'dust' particles.

For the critical overdensity, we adopt six different definitions: 2000, 1000, 500 and 200 $\rho_c(z)$, 200 $\rho_b$ and the virial value. We chose these values of overdensity since they are (or they are very close to) the commonly used ones: 200 $\rho_b$ is motivated by the spherical collapse model in an Einstein-de-Sitter universe and – together with 200 $\rho_c$ – is a popular choice (Tinker et al. 2008); moreover, 200 $\rho_c$ is often used to define galaxy cluster masses; 500 $\rho_c$ (and the higher overdensities) are used in X-ray analyses, and in general in observations that are able to resolve only the inner parts of haloes. The virial overdensity depends on redshift and cosmological model

(e.g. Bryan & Norman 1998); we use the numerical solutions of Eke, Cole & Frenk (1996).

The comoving density of the background is

$$ \rho_{\text{com}} \equiv \rho_b = \rho_c(0)\Omega_m(0) = \rho_c(z)\Omega_m(z), \quad (1) $$

where $\rho_c(0) \equiv 3H_0^2/8\pi G \simeq 2.775 \times 10^{11} h^{-1} M_\odot h^3 \text{Mpc}^{-3}$ is the critical density. For the Planck cosmology adopted in the main set of simulations, $\Delta_{\text{vir}}(z=0) \simeq 319 \rho_b$ is greater than 200 $\rho_b$ (corresponding to $\simeq 98 \rho_b$) and lower than all the other thresholds we consider. At high redshifts both 200 $\rho_b$ and 200 $\rho_c$ converge to the virial definition, the first from below and the second from above.

Fig. 1 shows a schematic representation of the haloes identified with our SO finder. In contrast to previous work, in which an FoF catalogue is used as the basis for subsequent SO identifications, we run our halo finder code from scratch for each threshold. Hence, although a halo in one catalogue may be present in another, this is not necessarily true. In particular, while the halo centred at 1 is common to all catalogues, the ones at 2 and 3 belong only to 1000 $\rho_c$ and to the virial catalogue, respectively. This happens, as can be noticed from their density profiles presented in the right part of the figure, because while halo 1 is dense enough to go from overdensity 2000 $\rho_c$ to 200 $\rho_b$ – so is in common to all catalogues – haloes 2 and 3 reach only 1000 $\rho_c$ and the virial overdensity, respectively, and are present only in those catalogues.

Halo identification in each simulation snapshot was done with the aim of studying the evolution the halo mass function. For each snapshot–density threshold combination, we saved a catalogue containing all the information about the identified haloes. So as not to be biased by the mass and force resolution, we only consider systems resolved with at least 300 particles (Macciò et al. 2007; Macciò, Dutton & van den Bosch 2008; Velliscig et al. 2015). Therefore, while the $z=0$ catalogues of each simulation contain many haloes
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic representation of the halo identifications in the particle density distribution at $z=0$. Different colours represent the various overdensities, from 2000 $\rho_c$ down to 200 $\rho_b$. Since our SO halo finder starts to grow the sphere starting from the densest particle, the halo centred at 1 is common to all the catalogues, whereas the ones at 2 and 3 belong only to 1000 $\rho_c$ and to the virial catalogue, respectively. This is clearer in the right-hand panel of the figure which shows the density profiles of the three systems.</figcaption></figure>

over a wide mass range, only the high-mass end of the higher redshift outputs is reliably measured.

We also ran an ellipsoidal overdensity (EO) finder (Despali, Tormen & Sheth 2013; Despali, Giocoli & Tormen 2014; Bonamigo et al. 2015), which we used for estimating the triaxial properties of the collapsed systems. In what follows we concentrate on the results for SO haloes; a brief summary of the corresponding results for ellipsoidal haloes can be found in Appendix A. These are broadly similar to spherical haloes, although the best-fitting mass function parameters differ slightly from those for the SO case.

### 3 MODEL FOR THE HALO MASS FUNCTION

Let $dn/dln M$ denote the comoving number density of haloes in a logarithmic bin $dln M$ around mass $M$. Then, the mass fraction in such haloes is

$$f(M)dM = \frac{M}{\rho_b} \frac{dn}{d\ln M} d\ln M. \qquad (2)$$

It is usual to define

$$\sigma^2(M, z) = \int \frac{dk}{k} \frac{k^3 P_{\text{lin}}(k, z)}{2\pi^2} W^2[kR(M)], \qquad (3)$$

where $P_{\text{lin}}(k, z)$ is the initial power spectrum extrapolated to redshift $z$ using linear theory, $W(x) \equiv 3j_1(x)/x$, and $R(M)$ is given by requiring $M/(4\pi R^3/3) = \rho_b$ (recall that $\rho_b$ is comoving, so it is independent of redshift). In principle, $\sigma(M, z)$ depends only on $P_{\text{lin}}(k, z)$ and the smoothing filter $W$. In practice, the number of Fourier modes that can be effectively sampled in the initial conditions depends on some computational limits (e.g. box size, number of particles, etc.). As a result, there are differences between the actual power spectra in a box and the theoretical mean value. Appendix B illustrates the impact on $\sigma^2(M)$.

We compute $\sigma^2(M)$ for each box using the actual realization of $P(k)$ in it. This reduces bias and scatter in the mass function, especially in the high-mass tail of each simulation, allowing us to reach higher precision (down from ten to a few per cent) particularly in the smaller boxes within which cosmic variance would otherwise contribute substantially. In particular, the difference between the theoretical input power spectrum and the one measured from each

simulation can be seen in Fig. B1, while the associated difference in the mass function is shown in the bottom panel of Fig. B2. Appendix B also discusses the impact of box size on the mass functions, and compares a number of different prescriptions that account for finite box-size effects on $\sigma^2(M)$.

Although a number of workers so far have parametrized the mass function in terms of $\sigma$ alone (i.e. they write $f(M)dM = f(\sigma)d\sigma$ and work with $f(\sigma)$), Sheth & Tormen (1999) were careful to parametrize in terms of

$$\nu \equiv \delta_c^2(z)/\sigma^2(M), \qquad (4)$$

where $\delta_c(z)$ is the critical linear theory overdensity $\delta_{\text{lin}}(z)$ required for spherical collapse divided by the growth factor (Carroll, Press & Turner 1992). This quantity, which depends weakly on $\Omega$ for the $\Lambda$CDM family of models, is well approximated by

$$\delta_{\text{lin}}(z) \approx \frac{3}{20} (12\pi)^{2/3} [1 + 0.0123 \lg \Omega(z)] \qquad (5)$$

(Kitayama & Suto 1996). The rationale for including it came from the fact that they used the virial density to find haloes, and the same model which predicts this virial density also predicts $\delta_c$. Courtin et al. (2011) highlight the fact that this factor becomes increasingly important at high masses; failure to include it can masquerade as non-universality.

For this reason, our reference parametrization for the halo mass function is that of Sheth & Tormen (1999):

$$f(M)dM = f(v)dv \qquad (6)$$

where

$$v f(v) = A \left(1 + \frac{1}{v'^p}\right) \left(\frac{v'}{2\pi}\right)^{1/2} e^{-v'^2/2} \qquad (7)$$

with $v' = av$. The parameters $(a, p, A_0)$ define the high-mass cut-off, the shape at lower masses, and the normalization of the curve, respectively. In addition, because $P_{\text{lin}}(k)$ determines the value of the mass variance $S(M) = \sigma^2(M)$ which enters in the definition of $v$, the initial power spectrum plays an important role. Since it enters in the denominator of $v$, one might say that the mass function is 'non-perturbative' in $P_{\text{lin}}(k)$.

### 4 THE UNIQUENESS OF THE VIRIAL OVERDENSITY

This section studies the universality of the halo mass function for haloes identified using the virial overdensity. We do so by finding the set of parameters $(a, p, A_0)$ that best fit the $z=0$ data from the Planck13 simulations. We then use measurements at other redshifts and cosmologies to test for universality. We start with the virial overdensity because we believe it to be the most physically motivated choice for identifying haloes. Later in the paper, we study haloes identified using the other SO catalogues.

Our model for the mass function, equation (7), has three free parameters $(a, p, A_0)$, whose values we adjust so as to minimize the Chi square with respect to the measured binned mass function:<fnref n="2" />

$$\chi^2(a, p, A_0) = \sum_i \frac{(\log(v f(v))_i - \log(v f(v))_{\text{fit}})^2}{e^2 \log(v f(v))_i}, \qquad (8)$$

<footnote n="2">² Note that Sheth & Tormen (1999) only varied $a$ and $p$; they fixed the value of $A$ by requiring that the integral over all masses give $\rho_b$.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Halo mass function at four redshifts, for all SO haloes with more than 300 dark matter particles in the SO catalogues, identified using the redshift-dependent virial density. Green, cyan, blue, magenta red and orange symbols show results for Ada, Bice, Cloe, Dora, Emma and Flora. Black dashed line, same in all panels, shows the result of fitting the $z = 0$ points to equation (7). The middle panels show the residuals (in log space) from this best fit. The dashed vertical lines show the maximum $v$ for which the bins contain at least 100 objects. The two solid vertical lines show the minimum $v$ at which we may expect differences of at least 5 per cent in the halo mass function accounting for 2LPT initial condition, adopting the rescaling of equation (11) by Reed et al. (2013), for the simulations starting at $z = 99$ and 124 (right and left lines, respectively). Red solid line shows the best-fitting mass function obtained using the data from all the snapshots with $z \le 1.25$. The bottom panels show residuals from the fit to the scaled counts from all snapshots up to $z = 1.25$. The departures at high redshift are reduced, so universality is even more pronounced.</figcaption></figure>

where the sum is over binned counts. The bins are equally spaced in $\log_{10}(v)$, with $\Delta\log_{10}(v) = 0.05$ and we neglect covariances between the binned counts when fitting.<fnref n="3" /> We discard bins with fewer than 30 objects (typically the high-$v$ bins), to set a limit on the Poissonian error $\epsilon_{\log(v/f(v))}$ in a bin. There were no significant differences in the best-fitting parameters when we repeated the analysis using only bins with at least 100 objects.

## 4.1 The halo mass function at different redshifts and overdensities

### 4.1.1 Virial haloes

The best-fitting values for the $z = 0$ virial overdensity halo counts are

$$ (a, p, A_0) = (0.794 \pm 0.005, 0.247 \pm 0.009, 0.333 \pm 0.001). \quad (9) $$

The dashed black curves, same in all panels of Fig. 2, show equation (7) with these parameters. The symbols in the different panels show the measured virial halo mass function at four different redshifts as labelled. Green, cyan, blue, magenta, red and orange points show results from the six simulations of the main set: Ada, Bice, Cloe, Dora, Emma and Flora. Note the excellent agreement among the simulations despite the very different volume that each samples. This is because the way we calculate $v(\delta_e, \sigma)$ largely eliminates the impact of cosmic variance.

The middle panels of Fig. 2 show logarithmic residuals from the best fit to the $z = 0$ counts. These indicate the goodness of our fit,

which was calibrated in log space. The residuals are close to zero at small and intermediate values of $v$. They are larger at high $v$, where the best fit predicts more haloes than we find in the simulations at $z = 0$. Some of this is due to the Poisson uncertainty in the small number counts; at higher redshifts the effect is reduced since the high $v$ tail becomes populated by lower mass haloes. The smallness of the residuals in the other panels indicates that the $z = 0$ model is a good description of the virial halo mass function at higher redshifts as well.

The vertical lines in the bottom panels show the regime of influence of two numerical effects: First, the dashed grey vertical line shows the maximum $v$ for which there are at least 100 (rather than 30) haloes in the bin. This matters only at high redshifts and in the high-mass regime and, as mentioned before, the best-fitting parameters are not significantly different; Secondly, we addressed the effect of not using 2LPT initial conditions for our simulations: we rescaled our points using equation (11) of Reed et al. (2013) (calibrated only for $1 \lesssim v \lesssim 5$), which gives an estimate of the bias in the halo mass function between ZA (Zel'dovich approximation) and 2LPT methods; the two solid vertical lines show the minimum $v$ at which we may expect differences of at least 5 per cent, for the simulations starting at $z = 99$ (right) and $z = 124$ (left). Such discrepancies would mainly affect the high-$z$ and high-$v$ data points; they do not affect the counts in the $z \le 1.25$ range we will adopt to calibrate the parameters of equation (7). The rescaling equation by Reed et al. (2013) and the ZA/2LPT difference are discussed in detail in Appendix B4.

From the figure and the discussion above, we conclude that there is no significant systematic deviation from universality within 8 per cent for $v \lesssim 10$ at redshifts $z \le 5$. Therefore, we can combine measurements at many different redshifts to increase the precision in the estimate of $(a, p, A_0)$. In addition to increasing the statistics,

<footnote n="3">³ We do not include a bin that counts the mass fraction which is not assigned to haloes. See Manera et al. (2010) for an algorithm which does not require binned counts, and does account for the unassigned mass fraction.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Change in the best-fitting mass function when one parameter at a time is changed by up to ±30 per cent (in increments of 10 per cent). Blue and red curves show cases where the parameter values are larger or smaller than the reference one.</figcaption></figure>

adding the high-redshift data fills-in the high-ν tail, allowing it to play a greater role in determining the best-fitting parameters. We determined the best-fitting parameters using the virial halo data from all the snapshots between z = 0 and 1.25 (a total of 15 snapshots). At z ≥ 1.25, the data no longer samples the whole mass function, preventing a reliable fit. The resulting best-fitting parameters are:

$$ 
\begin{aligned}
a &= 0.7663 \pm 0.0013 \\
p &= 0.2579 \pm 0.0026 \\
A_0 &= 0.3298 \pm 0.0003.
\end{aligned} 
\qquad (10) $$

The differences between this all-z fit and the previous z = 0 fit are of the order of $\Delta a \simeq 3$ per cent, $\Delta p \simeq 4$ per cent and

$\Delta A_0 \simeq 1$ per cent. This leads to per cent level differences in the mass function which reduce the residuals at high redshift. The red solid line in each panel of Fig. 2 shows this all-z best fit: it can hardly be distinguished from that for z = 0 only (black dashed curve) but, as shown by the residuals in the lower panels, it traces the high-ν part of the mass function more closely.

Fig. 3 shows how the virial mass function changes as each of the parameters (a, p, A₀) is varied by ±30 per cent, while the other two are kept fixed. While the mass function is sensitive to both A₀ and a – which modify the normalization and the high-mass cut-off, respectively – this is less true for p. Even if p changes by ±30 per cent, the mass function changes by 5 per cent at most, for $ν \lesssim 0.1$ and $ν > 10$. In what follows, we use the values in equation (10) as our reference model.

### 4.1.2 Other overdensity thresholds

We now study the mass functions associated with the other five density thresholds. It is tempting to identify these thresholds with an ‘effective formation redshift’, and hence with an effective value of $\delta_c(z)$ when defining $ν$. Since the role of a in equation (7) is simply to rescale $ν$, one might wonder if these other halo definitions lead to universal mass functions which differ from those for the virial overdensity only in the value of a.

Table 3 shows the values of the best-fitting parameters at z = 0 for the different threshold densities. Clearly the other parameters, p and A₀, also depend strongly on how the haloes were identified. Moreover, as we will show below, universality clearly does not hold for any of these other definitions: in all cases, the best-fitting parameters for the higher redshift counts depart significantly from those at z = 0.

To illustrate this, Fig. 4 shows the halo mass functions at z = 0, 1, 2 and 5 for haloes identified using 200 $\rho_b$ as density threshold. The dashed curve, which is the same in all four top panels, shows the best fit to the z = 0 data points. The middle panels, which show the residuals with respect to this fit, show that it overestimates the counts at all higher redshifts. These departures from universality are in agreement with previous work on 200 $\rho_b$ haloes (e.g. Tinker et al. 2008). The solid curves in the upper panels show the result of rescaling our universal virial counts as described in the next

<table><caption>**Table 3.** Dependence of best-fitting parameters on the overdensity used to identify SO haloes. In the top part, we show the three parameters calculated at z = 0; then we report those obtained by fitting the virial halo counts from all snapshots up to z = 1.25 (i) of the Planck simulations and (ii) of all the simulations with different cosmologies together. The bottom row shows the parameters which best fit the $M_{vir} \ge 3 \times 10^{13} h^{-1} M_{\odot}$ counts from all cosmologies and all redshifts.</caption><thead><tr><th>ρ (SO)</th><th>a</th><th>p</th><th>A</th></tr></thead><tbody><tr><td></td><td colspan="3" style="text-align: center;">z = 0</td></tr><tr><td>200 ρ<sub>b</sub></td><td>0.739 ± 0.005</td><td>0.206 ± 0.008</td><td>0.360 ± 0.001</td></tr><tr><td>Δ<sub>vir</sub></td><td>0.794 ± 0.005</td><td>0.247 ± 0.009</td><td>0.333 ± 0.001</td></tr><tr><td>200 ρ<sub>c</sub></td><td>0.903 ± 0.006</td><td>0.322 ± 0.009</td><td>0.287 ± 0.001</td></tr><tr><td>500 ρ<sub>c</sub></td><td>1.166 ± 0.009</td><td>0.344 ± 0.012</td><td>0.236 ± 0.001</td></tr><tr><td>1000 ρ<sub>c</sub></td><td>1.462 ± 0.012</td><td>0.349 ± 0.015</td><td>0.197 ± 0.001</td></tr><tr><td>2000 ρ<sub>c</sub></td><td>1.821 ± 0.017</td><td>0.413 ± 0.017</td><td>0.158 ± 0.001</td></tr><tr><td></td><td colspan="3" style="text-align: center;">All z – Planck cosmology</td></tr><tr><td>Δ<sub>vir</sub></td><td>0.7663 ± 0.0013</td><td>0.2579 ± 0.0026</td><td>0.3298 ± 0.0003</td></tr><tr><td></td><td colspan="3" style="text-align: center;">All z and cosmologies</td></tr><tr><td>Δ<sub>vir</sub></td><td>0.7689 ± 0.0011</td><td>0.2536 ± 0.0026</td><td>0.3295 ± 0.0003</td></tr><tr><td></td><td colspan="3" style="text-align: center;">All z and cosmologies – cluster counts: M<sub>vir</sub> &gt; 3 × 10<sup>13</sup> M<sub>☉</sub> h<sup>-1</sup></td></tr><tr><td>Δ<sub>vir</sub></td><td>0.8199 ± 0.0010</td><td>0</td><td>0.3141 ± 0.0006</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Same as Fig. 2, but for haloes identified using 200 $\rho_b$ instead of $\Delta_{vir}$. Middle panels show residuals with respect to the best fit at $z=0$ (dashed curve, same in all top panels); these show that the mass function is not universal across all redshifts. Lower panels show residuals with respect to our rescaled model [equation (12); solid curves, different in each top panel].</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Same as Fig. 2, but for haloes identified using 200 $\rho_c$; middle panels show that the universality is broken in the opposite sense to when the threshold was 200 $\rho_b$, and bottom panels show that our rescaled model accounts for this quite well.</figcaption></figure>

section; the bottom panels show the residuals with respect to it, which provides a much better fit.

Fig. 5 shows a similar analysis of haloes identified using 200 $\rho_c(z)$. In this case, the best-fitting relation at $z=0$ underestimates the counts at higher $z$. Finally, Fig. 6 shows haloes identified using threshold values of 500, 1000 and 2000 $\rho_c(z)$, at redshifts $z=0$ and 1. The trends are qualitatively similar to those for 200 $\rho_c(z)$, with the $z=0$ fit underestimating the counts at higher $z$.

## 4.2 The halo mass function for different cosmologies

Having established that the virial mass function is universal with respect to redshift, we now test if it is universal across other back-

ground cosmological models. To do so, we use the halo catalogues in our secondary set of 512³ simulations, whose properties are summarized in Table 2.

First, for virial haloes, we used the best-fitting parameters of equation (10) which were calibrated using the 1024³ particle simulations of a Planck13 cosmology. The four panels in Fig. 7 show the logarithmic difference from this best fit at redshifts $z=0, 1, 2, 5$. Following the colour code of Table 2, grey points are for $\Omega_m = 0.2$, brown for $\Omega_m = 0.4$ and light blue for $\Omega_m = 0.272$ (the WMAP7 cosmology); squares and triangles represent $\sigma_8 = 0.7$ and 0.9, respectively. For halo mass definitions that are close to the virial value the residuals are relatively small, indicating that our relation – calculated exclusively from the Planck13 cosmology
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 6. Same as Fig. 2, but for haloes identified using 500, 1000 and 2000 $\rho_c$ instead of $\Delta_{vir}$, at $z=0$ and 1.</figcaption></figure>

data – remains valid for the other cosmological models as well. In addition, although we do not show this explicitly, for these cosmological models also, universality is broken when considering other overdensities. This demonstrates universality of the virial relation

with respect to other cosmological models. However, departures from this universality may arise if one considers more extreme departures from Planck13. For example, (Courtin et al. 2011) explore models in which the growth factor, growth history and virial density differ more radically from that of the Planck cosmology, and find correspondingly larger departures from universality.

Fig. 8 shows the results from all our simulations and redshifts together. Different colours represent the different simulations (with the same colour code as in previous figures), and different symbols distinguish the four redshifts (and not the different simulations). The best fit to all the virial halo catalogues out to $z = 1.25$ has

$$ 
\begin{aligned}
a &= 0.7689 \pm 0.0011 \\
p &= 0.2536 \pm 0.0026 \\
A_0 &= 0.3295 \pm 0.0003,
\end{aligned} 
\qquad (11) $$

and is shown by the solid black line. Lower panel shows the logarithmic residuals. These best-fitting parameter values equal – at the per mil level – the values obtained in the last section from our first set of simulations, which assumed a Planck13 cosmology (equation 10), thus further confirming the universality as a function of the background cosmological model.

Table 3 compares the three sets of best-fitting values for the virial haloes. The bottom row shows the parameters which best fit the counts of cluster-mass haloes: $M_{vir} \ge 3 \times 10^{13} h^{-1} M_{\odot}$. Note that at these high masses $p = 0$, so that the mass function is like that of Press & Schechter (1974). We discuss this further in Appendix C in the context of cosmological constraints in the $\Omega_{m}-\sigma_8$ plane from cluster counts.

Covariances between the best-fitting parameters are shown in Fig. 9. Contours show $1\sigma$, $2\sigma$ and $3\sigma$ levels for each pair of parameters. While $A_0$ and $a$ are not strongly degenerate, $p$ correlates with both $a$ and $A_0$. Manera et al. (2010) show that such correlations can be understood as resulting from requiring the model to reproduce the total measured mass fraction in haloes (the mass fraction is measured with much greater precision than is the detailed shape of the mass function; also see Fig. 3).

## 5 RESCALING THE MASS FUNCTION: A UNIVERSAL PARAMETRIZATION

We now describe a simple method that allows one to derive the halo abundances associated with other overdensities at any redshift by a straightforward rescaling of the virial halo mass function.

We calculated the best-fitting parameters for all the six halo catalogues at all outputs between $z=0$ and 1.25. For the overdensities $\ge 500 \rho_c(z)$, we only fitted the data points out to $z \sim 0.4$ so as to be sure of having good statistics for both the shape and the cut-off of the mass function, i.e. $p$ and $a$, respectively. At higher redshifts the small number of haloes which probe a smaller range of $\nu$ cannot break degeneracies in the determination of the three parameters. Nevertheless, even after restricting to $z \sim 0.4$, we obtain very good results. As the reference virial halo mass function we used the fit from equation (11). Fig. 10 shows how the best-fitting parameters vary as a function of $x = \log (\Delta(z)/\Delta_{vir}(z))$.

As noted in Tinker et al. (2008, albeit for a different functional form), the best-fitting parameters ($a$, $p$, $A_0$) are smooth functions of the critical density threshold. As can be seen in the top panel of Fig. 10, the normalization $A_0$ decreases linearly with the overdensity – a natural consequence of the decrease in halo mass – while $a$ and $p$ both increase with threshold overdensity. The trends
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Cosmology independence of the virial halo mass function. We show the logarithmic difference from the best fit calculated using the Planck cosmology for redshifts z = 0, 1, 2 and 5. Following the colour code of Table 2, grey points indicate the results of the simulations with Ω_m = 0.2, brown for Ω_m = 0.4 and light blue for Ω_m = 0.272; squares and triangles represent σ_8 = 0.7 and 0.9, respectively.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Universality of the virial halo mass function. Symbols show measurements from all our simulations – both from the main and secondary set – at four redshifts: z = 0, 1, 2, 5. Although we show measurements at these four redshifts, the fit, which was calibrated on all simulations, used only the z ≤ 1.25 snapshots.</figcaption></figure>

can be described using linear or quadratic polynomials:

$$ 
\begin{align} 
a &= 0.4332x^2 + 0.2263x + 0.7665, \\ 
p &= -0.1151x^2 + 0.2554x + 0.2488, \\ 
A &= -0.1362x + 0.3292, 
\end{align} 
\tag{12} $$

shown by the blue solid lines in Fig. 10. While a and A behave very regularly, the determination of p is less certain. This does not have a big influence on the rescaling method, since the mass function is less sensitive to p than to the other parameters (see Fig. 3). The red dashed curves in Fig. 10 show the trends obtained by fitting to the z = 0 counts only:

$$ 
\begin{align} 
a &= 0.3881x_0^2 + 0.2776x_0 + 0.7837, \\ 
p &= -0.07459x_0^2 + 0.2016x_0 + 0.2518, \\ 
A &= -0.1337x_0 + 0.3315, 
\end{align} 
\tag{13} $$

where $x_0 \equiv \log(\Delta(z_0)/\Delta_{vir}(z_0 = 0))$. These trends are very similar to the previous ones.

The above relations show that both a and p increase with increasing threshold, qualitatively consistent with Manera et al. (2010) who studied FoF rather than SO haloes. (They found that a and p increase as the FoF linking length is decreased, and it is well known that shorter linking lengths return denser haloes.) In addition, upon noting that a multiplies δ_c(z) in equation (7), the increase of a with density threshold is qualitatively consistent with the notion that the denser inner parts of a halo virialized at higher redshift, so mass functions for higher density thresholds are qualitatively like those for higher redshifts. Unfortunately, the agreement is not quantitative. To see this, consider the thresholds 2000 and 200 ρ_c which differ by a factor of 10. If the mass associated with the denser threshold virialized at (1 + z_2000)³ ~ 10, we would expect the associated δ_c to be larger by a factor of (1 + z_2000) ~ 10¹/³. If this increase is to be provided by increasing a, then a must be larger by

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Covariance between fitted parameters of the virial halo mass function (the 1024³ and 512³ runs, and all redshifts to z ≤ 1.25): contours show 1σ, 2σ and 3σ reference levels for each pair of parameters.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Dependence of best-fitting parameters *a*, *p*, $A_0$ on $x = \log(\Delta(z)/\Delta_{vir}(z))$ for $z = 0$ to 1.25. Different coloured symbols represent different overdensities; smooth curve shows equation (12).</figcaption></figure>

$10^{2/3}$. This is considerably larger than the ratio of *a* values for 2000 and 200 $\rho_c$ in Table 3.

Even though we do not have a quantitative physical understanding of the scaling in equations (13) and (12) – a question we return to in Section 6 – they can still be used to predict the mass functions at overdensities of current interest that are not considered in our work. These precise fitting formulae are able to follow the evolution of the non-virial mass functions at all redshifts, at the level of a few

per cent – which is comparable to the intrinsic uncertainty in the mass function.

To show this, we return to the solid lines and the bottom panels in Figs 4–6. The black solid curves in Figs 4–6 show equation (12), recall that the dashed curves show the best $z = 0$ fits) and the bottom panels show the logarithmic difference from our model. The agreement is indeed very good: the residuals are generally smaller than a few per cent; this is quite acceptable given the intrinsic scatter of the mass function is of this order. The middle and bottom panels show that our model curves are able to account both for the change in the normalization and the tilt of the mass function as $z$ changes.

To conclude: our results provide an efficient way for predicting the halo mass function for any redshift and overdensity in a Planck ΛCDM-like model. Once $\log(\Delta(z)/\Delta_{vir}(z))$ has been calculated, it is straightforward to obtain the other parameters.

Returning to the other cosmological models, Fig. 11 shows the logarithmic differences from our reference rescaling models for all the mass functions of the secondary set, at redshifts $z = 0, 1$. At both redshifts, the residuals are very small showing that our fitting relations also work for any cosmological model once mass and redshift are appropriately rescaled to $v$.

## 6 MATCHED HALOES

So far, we used the halo catalogues derived independently at each overdensity threshold from the particle density distribution saved in the snapshot files. As already noted, the number of haloes identified using different overdensities at a certain redshift is not the same. This means that we are not looking at a ‘rescaled version’ of the virial population, but at different halo samples (see Fig. 1). For example, virial haloes may contain many smaller higher density peaks, each of which could correspond to a massive substructure within the virial radius. On the other hand, the 200 $\rho_b$ threshold is lower than the virial one and so – at low redshifts – the identified haloes will be larger and more massive than the virial ones. In general, in this latter case we may have fewer haloes than in the virial catalogues, since the particles corresponding to more than one smaller virial halo may be included in a single 200 $\rho_b$ halo (cf. the virial halo 3 in Fig. 1).

In this section, we analyse the halo mass function of ‘matched haloes’ which we create as follows. For each virial halo, we select from the different halo overdensity catalogues the object whose centre of mass is closest to that of the virial halo. For the haloes identified with a threshold that is denser than virial, we only keep the object which corresponds most closely to the virial one. To ensure good resolution (and hence a match) in the inner regions, we only considered virial haloes with more than $10^3$ particles. In this way, we define a series of approximately concentric spheres of ever higher density within the virial radius: a density profile.

If we average together all the profiles for a narrow bin in virial mass, then we can use this mean profile, and the scatter around it, to model the ‘matched halo’ mass function associated with different overdensity thresholds $\Delta \ge \Delta_{vir}$, because

$$ \frac{dn_{\Delta}}{dM_{\Delta}} = \int dM_{\text{vir}} \frac{dn_{\text{vir}}}{dM_{\text{vir}}} p(M_{\Delta}|M_{\text{vir}}). \quad (14) $$

If $p(M_{\Delta}|M_{\text{vir}})$ is sharply peaked around $M_{\Delta} = g_{\Delta}(M_{\text{vir}}) M_{\text{vir}}$ say, then

$$ \frac{dn_{\Delta}}{dM_{\Delta}} = \frac{dn_{\text{vir}}}{dM_{\text{vir}}} \frac{dM_{\text{vir}}}{dM_{\Delta}}, \quad (15) $$
<figure><img src="image_15.png" /><figcaption>Figure 11. Mass functions for other cosmological models: logarithmic difference from our best-fitting models in the secondary set of 512³ simulations, at two different redshifts ($z = 0, 1$) and for all the overdensities. Grey, brown, light blue show results for $\Omega_m = 0.2, 0.4$ and $0.272$. Squares and triangles show results for $\sigma_8 = 0.7$ and $0.9$. For the virial haloes, we use the fit obtained from stacking all $z$ outputs of the Planck cosmology. For the other cases, the residuals are with respect to equation (12).</figcaption></figure>

so that

$$[\nu f(\nu)]_{\Delta} = g_{\Delta}(M_{\text{vir}}) [\nu f(\nu)]_{\text{vir}} \quad (16)$$

as in Bocquet et al. (e.g. 2015). But if the scatter around the mean profile is not negligible, then the full convolution of equation (14) must be performed.

In principle, $p(M_{\Delta}|M_{\text{vir}})$ could be determined directly by fitting an NFW profile Navarro, Frenk & White (1996) to each virial halo in the stack. This functional form has just one free parameter, the concentration $c$, so $p(M_{\Delta}|M_{\text{vir}})$ is simply related to the mean $c$ at fixed $M_{\text{vir}}$, and the scatter around this mean, which is known to be lognormal with width $\sigma_{\ln c|m} = 0.25$.

In practice, we have not measured the profile shapes. Rather, we use the fact that the concentration of a halo is related to its mass accretion history (Zhao et al. 2009), which we also do not measure. Rather, we estimate it, and hence $c$, using the model of Giocoli et al. (2012). We then use the mean $c$ along with the assumption that the scatter is negligible, to compute the value of $g_{\Delta}$ to insert in equation (16) for the desired $\Delta$.

The symbols in the left-hand panels of Fig. 12 show the mass function of the haloes in the different catalogues matching the virial systems, at redshift $z = 0$ and 1 for the catalogues constructed considering 200, 500, 1000 and 2000 $\rho_c$ as thresholds. The solid curves, following the corresponding data points, show the virial mass function rescaled to the different overdensities.

The panels on the right of Fig. 12 show the logarithmic residuals of the data points with respect to the corresponding rescaled models. They all remain well below 10 per cent – apart from the high-$\nu$ tail – indicating that the precision of the rescaled mass functions is of the same order as of the ‘true’ virial mass function. Some of the discrepancy may be caused by the fact that not all haloes follow NFW profiles all the way down to the very centre (Einasto 1965; Retana-Montenegro et al. 2012; Ludlow et al. 2013; Dutton & Macciò 2014).

The case of 200 $\rho_b$ must be treated separately, since the more massive haloes identified using 200 $\rho_b$ may include more than one smaller virial overdensity haloes in their outskirts which we exclude (e.g. object 3 in Fig. 1). Nevertheless, in this case also, the rescaled relation captures the behaviour of the mass function measured in the simulations with a precision comparable to the other overdensities. Recall that the rescaled mass function at 200 $\rho_b$ at $z = 1$ is almost the same as the virial mass function, since $\Delta_{\text{vir}}$ at high redshift is nearly equal to 200 $\rho_b$.

Fig. 13 shows how different these rescaled mass functions are from the original fits of Table 3. The upper panels show the two mass functions – the rescaled one in dashed lines and the original ones in solid lines – for all the density thresholds; the lower panel shows the differences in halo counts between the two cases. Note that the matched mass functions ‘lose’ haloes: because some of the smaller denser systems are halo substructures. This is true for all the overdensities except for 200 $\rho_b$: as we said, in this case we are outside the virial haloes and so we may in fact have more haloes than those found at the virial overdensity.

## 7 COMPARISON WITH PREVIOUS WORK

Fig. 14 shows the residuals between our virial mass function $\nu f(\nu)$ – the one calibrated from all $z \le 1.25$ outputs of the 1024³ and the 512³ simulations – and the mass functions derived in previous work (as indicated). The shaded region in each panel shows the range in $\nu$, at $z = 0$, over which the various authors have calibrated their mass functions. The two vertical dot-dashed lines show the range in $\nu$ over which we have performed our calibration.

To facilitate comparison, the main features of these works are listed in Table 4. Most of these authors used an FoF algorithm to identify the haloes in their simulations: while $b = 0.2$ is the most common choice for the linking length, Manera et al. (2010) use three different values ($b = 0.15, 0.168, 0.2$); moreover, Warren
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Top left: $z=0$ halo mass functions of the matched haloes. Different colours show different overdensities, while the black solid curves show the corresponding rescaled mass functions, obtained from equation (16). Top right: logarithmic residuals from the rescaled mass functions of the upper panel. The scatter with respect to the rescaled mass function is of the same order of the intrinsic scatter of the mass function (shown by the virial case). Bottom: same as top, but for $z=1$.</figcaption></figure>

et al. (2006) discuss how to correct the FoF masses, claiming that the identification is influenced by how well an halo is resolved and thus the same linking length may identify structures with different enclosed overdensities. The same argument is discussed also in Courtin et al. (2011) and More et al. (2011) and they predict a linking length equivalent to our virial overdensity (at $z=0$) of 0.193 and 0.206, respectively. Sheth & Tormen (1999) and Tinker et al. (2008) use SO algorithms, while Watson et al. (2013) compares the results of the two methods (in Fig. 14, we show only their universal FoF fit). At intermediate masses all the analytical mass functions agree quantitatively. Larger differences arise for more massive systems: in this range, the precision of the fit is strongly affected by the

resolution of the simulation and consequently by the number of high-mass haloes.

The main purpose of our work is to analyse the impact of the density threshold chosen to identify the haloes on the universality of the mass function: the majority of the identification algorithms use a threshold $\Delta$ of 178 or 200 for all redshifts (e.g. Tinker et al. 2008 and Watson et al. 2013). Since we have argued that ignoring the redshift evolution of $\Delta$ causes most of the observed non-universality one must carefully match density thresholds when comparing different works. This is shown more clearly in the lower panels of Fig. 14, where we compare our results with those of Tinker et al. (2008): we tried to match their SO thresholds as closely as possible to ours,
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Difference between the mass functions of Table 3 (solid lines) and the rescaled ones from the matched haloes as in Fig. 12 (dashed lines); the lower panel shows the residuals between the two relations on the halo counts.</figcaption></figure>

to test their compatibility. The figure shows – from top to bottom – the residuals between: (i) their $\Delta = 200$ fits with our virial fits; (ii) their $\Delta = 300$ fits with our virial fits; (iii) their $\Delta = 200$ fits with our $200 \rho_b$ fits. The agreement for cases (ii) and (iii) is much better than for (i). Since our virial overdensity at $z=0$ is closer to $\Delta = 300$ than 200, this shows that different analyses converge provided one compares compatible halo identification schemes (Knebe et al. 2013). Moreover, not only are our $200 \rho_b$ fits close to the main result of Tinker et al. (2008), but the non-universality they report is compatible with the redshift evolution of our $200 \rho_b$ mass function (Fig. 4): the systematic deviation from universality is primarily due to the choice of the overdensity threshold used in the halo identification method. The agreement, once overdensities have been correctly matched, is even more reassuring when one notes that Tinker et al. (2008) fit their counts to a functional form that differs from our equation (7).

Additional reassurance comes from the fact that, while one may expect trends arising from differences in halo identification methods, our SO results lie in between the $l = 0.168$ and 0.15 FoF results of Manera et al. (2010), just as they should.

As a final remark we note that, even when the computational schemes are theoretically similar, the mass function for structures identified with different codes (in the same simulation!) can differ by up to 10 per cent (Knebe et al. 2011). This might explain some of the discrepancies between our best-fitting model and those of other works. Finally, we underline that there could be other effects depending on the assumed cosmological parameters as explored by Murray et al. (2013).

<figure><img src="image_18.png" /><figcaption>**Figure 14.** Comparison with some previous work. In all panels, the dotted vertical grey lines show our range in mass (for virial haloes at $z < 1.25$) and the grey region show an estimate of the range in mass of the other works; we calculated this last converting their mass ranges in $\nu$, so it may not be perfectly the same but it gives an estimate of the overlapping regions. Top: we show the relative residuals of our virial mass function $\nu f(\nu)$ with respect to the models of other works. From top to bottom, we compare our model with Sheth & Tormen (1999), Jenkins et al. (2001), Warren et al. (2006), Watson et al. (2013) and Manera et al. (2010). For this last, we show the results obtained from FoF catalogues with linking lengths of 0.2, 0.168 and 0.15, which they fit to the same functional form we do. Our results lie between their fits to the smaller linking catalogues, as they should, over the whole range of masses. Bottom: we compare our results with those of Tinker et al. (2008) at different overdensities. First, we show the relative residuals of their $\Delta = 200$ reference model with our best fit calculated at the virial density; secondly, we compare their model with our best fit at $200 \rho_b$. Finally, we test the opposite situation, comparing our virial model with their mass function at $\Delta = 300$, an overdensity closer to the value of $\Delta_{vir}$ adopted in this work.</figcaption></figure>

# 8 DISCUSSION AND CONCLUSIONS

In this paper, we analysed a set of cosmological simulations in order to study the halo mass function and its dependence on redshift, cosmology and the halo identification method. In what follows we summarize our main results.
<table><caption>**Table 4.** Schematic (i.e. incomplete) summary of the main features of the other works to which we compare our results in Fig. 14 (Sheth & Tormen 1999; Jenkins et al. 2001; Warren et al. 2006; Tinker et al. 2008; Manera et al. 2010; Watson et al. 2013). We list the resolution and scale of their simulations and the main algorithm with the associated threshold used to identify the haloes.</caption><thead><tr><th></th><th>N<sub>part</sub></th><th>Box (h<sup>-1</sup> Mpc)</th><th>Cosmological models</th><th>Algorithm</th><th>Threshold</th></tr></thead><tbody><tr><td>Sheth & Tormen 1999</td><td>(3 sim) 256<sup>3</sup></td><td>85, 141</td><td>S/O/ΔCDM</td><td>SO</td><td>Δ = Δ<sub>vir</sub></td></tr><tr><td>Jenkins et al. 2001</td><td>256<sup>3</sup>, 512<sup>3</sup>, 10<sup>9</sup></td><td>84–3000</td><td>τ/ΔCDM</td><td>FoF</td><td>b=0.2</td></tr><tr><td>Warren et al. 2006</td><td>1024<sup>3</sup></td><td>96–3072</td><td>ΔCDM</td><td>FoF</td><td>b corrected</td></tr><tr><td>Watson et al. 2013</td><td>3072<sup>3</sup>–6000<sup>3</sup></td><td>11–6000</td><td>ΔCDM (WMAP5)</td><td>FoF (SO)</td><td>b=0.2 (Δ=178)</td></tr><tr><td>Manera et al. 2010</td><td>(49 sim) 640<sup>3</sup></td><td>1280</td><td>ΔCDM</td><td>FoF</td><td>l=0.15, 0.168, 0.2</td></tr><tr><td>Tinker et al. 2008</td><td>512<sup>3</sup>–1024<sup>3</sup></td><td>80–1280</td><td>ΔCDM (WMAP1–3)</td><td>SO</td><td>Δ=n x Δ<sub>bg</sub></td></tr></tbody></table>

(i) We demonstrated and confirmed the universality of the virial halo mass function, by comparing the measured mass function at many different redshifts and in simulations with different cosmological models. As stated in Sheth & Tormen (1999) and Courtin et al. (2011), the halo mass function at any time, and for any cosmological model, can be well described by a single functional form once mass and redshift are appropriately parametrized in terms of $v = \delta_c^2(z)/\sigma^2(M)$.

(ii) We showed that only the virial overdensity leads to a universal halo mass function: most of the non-universality seen in other works arises from not using the virial value when identifying haloes. Commonly used values of 178 or 200× the critical or background density induce non-universal trends.

(iii) We derived three different sets of the best-fitting parameters for the virial halo mass function:

(a) $a = 0.794$, $p = 0.247$ and $A = 0.333$ – from the $z = 0$ virial mass function in the Planck cosmology;

(b) $a = 0.7663$, $p = 0.2579$ and $A = 0.3298$ – from the virial mass function in the Planck cosmology, using all the data points up to $z = 1.25$ (15 snapshots × 6 runs);

(c) $a = 0.7689$, $p = 0.2536$ and $A = 0.3295$ – from the mass functions in all cosmologies, using all the data points up to $z = 1.25$ (15 snapshots × 16 runs).

The last two are in remarkable agreement – at the per mil level for the best-fitting parameters, and sub per cent for the mass function), illustrating the level of the universality as a function of cosmology. In general, the normalization is the most stable parameter, while $a$ and $p$ change between the first and the other two cases: when using only the $z = 0$ points, the high-$v$ tail is less well resolved, and so the fit is less precise in the determination of the full shape of the mass function.

(iv) We presented a simple rescaling method which allows one to estimate the three parameters of the fitting function through first- or second-order scaling relations. The three parameters ($a$, $p$, $A_0$) are smooth functions of the overdensity, for any redshifts. Equations (12) are able to describe the change in slope and normalization of the mass function (as a function of redshift and overdensity) with good precision (Figs 4, 5, 6 and 11).

(v) Finally, we studied the mass function of ‘matched haloes’: the counterparts of virial haloes at different overdensities. We provide an efficient rescaling method, equation (16), which uses knowledge of the mass density profile and the concentration–virial mass relation, to estimate their mass function with a high precision (Fig. 12).

These rescaling methods are useful for comparing analyses which use different definitions of halo mass. In particular, they can be used to translate our universal virial halo mass functions to the non-universal form associated with halo definitions which are closer to those commonly used in observational studies – ranging from the X-ray and Sunyaev-Zel'dovich (SZ) to the visible and near-infrared.

We conclude that – over the range of redshifts and cosmological models in our simulation set – the virial halo mass function is, to within 5–8 per cent for a wide range of masses, a universal function of redshift and cosmology. Non-universal behaviour can be an artefact induced by the halo identification method and by the choice of the overdensity threshold. Other, true departures from universality may be sought in the other components of the universe or in more extreme cosmological models. Future extremely well-resolved simulations should allow a per cent level estimate of the universality of the mass function.

## ACKNOWLEDGEMENTS

GD, CG and RA thank Marceau Limousin and LAM (Laboratoire d'Astrophysique de Marseille), for organizing and allowing the meeting in Marseille where this work started, and Y. Rasera and P.S. Corasaniti for comments on an early version. CG thanks CNES for financial support. GD has been partially financed by the ‘Ing. A. Gini’ Fellowship. This work has been carried out thanks also to the support of the OCEVU Labex (ANR-11-LABX-0060) and the A*MIDEX project (ANR-11-IDEX-0001-02) funded by the ‘Investissements d'Avenir' French government programme managed by the ANR. The Cosmology Group in Padova thanks Vincenzo Mezzalira for installing and managing the ‘Nemo’ cluster. We thank the anonymous referee for his/her useful comments that considerably improve the presentation of our results. Great gratitude goes to the cosmology group at UPENN – Philadelphia – from GD and CG for the hospitality.

## REFERENCES

* Angulo R. E., Springel V., White S. D. M., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046

* Bocquet S., Saro A., Dolag K., Mohr J. J., 2015, preprint (arXiv:1502.07357)

* Boldrin M., Giocoli C., Meneghetti M., Moscardini L., Tormen G., Biviano A., 2015, preprint (arXiv:1505.03515)

* Bonamigo M., Despali G., Limousin M., Angulo R., Giocoli C., Soucail G., 2015, MNRAS, 449, 3171

* Bond J. R., Cole S., Efstathiou G., Kaiser N., 1991, ApJ, 379, 440

* Borgani S., Kravtsov A., 2011, Adv. Sci. Lett., 4, 204

* Bryan G. L., Norman M. L., 1998, ApJ, 495, 80

* Carroll S. M., Press W. H., Turner E. L., 1992, ARA&A, 30, 499

* Castorina E., Sefusatti E., Sheth R. K., Villaescusa-Navarro F., Viel M., 2014, J. Cosmol. Astropart. Phys., 2, 49

* Corasaniti P. S., Achitouv I., 2011, Phys. Rev. D, 84, 023009

* Courtin J., Rasera Y., Alimi J.-M., Corasaniti P.-S., Boucher V., Füzfa A., 2011, MNRAS, 410, 1911

* Crocce M., Pueblas S., Scoccimarro R., 2006, MNRAS, 373, 369

* Crocce M., Fosalba P., Castander F. J., Gaztañaga E., 2010, MNRAS, 403, 1353

* Cui W., Baldi M., Borgani S., 2012, MNRAS, 424, 993

* Cui W., Borgani S., Murante G., 2014, MNRAS, 441, 1769

* Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
Del Popolo A., Gambera M., 1998, A&A, 337, 96
Del Popolo A., Gambera M., 1999, A&A, 344, 17
Despali G., Tormen G., Sheth R. K., 2013, MNRAS, 431, 1143
Despali G., Giocoli C., Tormen G., 2014, MNRAS, 443, 3208
Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359
Einasto J., 1965, Tr. Astrofiz. Inst. Alma-Ata, 5, 87
Eke V. R., Cole S., Frenk C. S., 1996, MNRAS, 282, 263
Ettori S., Morandi A., Tozzi P., Balestra I., Borgani S., Rosati P., Lovisari L., Terenziani F., 2009, A&A, 501, 61
Evrard A. E. et al., 2008, ApJ, 672, 122
Giocoli C., Pieri L., Tormen G., 2008, MNRAS, 387, 689
Giocoli C., Meneghetti M., Ettori S., Moscardini L., 2012, MNRAS, 426, 1558
Giocoli C., Marulli F., Baldi M., Moscardini L., Metcalf R. B., 2013, MNRAS, 434, 2982
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
Kitayama T., Suto Y., 1996, MNRAS, 280, 638
Knebe A. et al., 2011, MNRAS, 415, 2293
Knebe A. et al., 2013, MNRAS, 435, 1618
Komatsu E. et al., 2011, ApJS, 192, 18
Lacey C., Cole S., 1993, MNRAS, 262, 627
Lacey C., Cole S., 1994, MNRAS, 271, 676
Laureijs R. et al., 2011, preprint (arXiv:1110.3193)
Lewis A., Challinor A., Lasenby A., 2000, Astrophys. J., 538, 473
Longair M. S., ed., 1998, Galaxy Formation. Astrophysics Library, Springer
Ludlow A. D. et al., 2013, MNRAS, 432, 1103
Macciò A. V., Dutton A. A., van den Bosch F. C., Moore B., Potter D., Stadel J., 2007, MNRAS, 378, 55
Macciò A. V., Dutton A. A., van den Bosch F. C., 2008, MNRAS, 391, 1940
Manera M., Sheth R. K., Scoccimarro R., 2010, MNRAS, 402, 589
Mo H., van den Bosch F. C., White S., 2010, Galaxy Formation and Evolution. Cambridge Univ. Press, Cambridge
More S., Kravtsov A. V., Dalal N., Gottlöber S., 2011, ApJS, 195, 4
Moreno J., Giocoli C., Sheth R. K., 2008, MNRAS, 391, 1729
Murray S. G., Power C., Robotham A. S. G., 2013, MNRAS, 434, L61
Musso M., Sheth R. K., 2012, MNRAS, 423, L102
Musso M., Sheth R. K., 2014a, MNRAS, 438, 2683
Musso M., Sheth R. K., 2014b, MNRAS, 439, 3051
Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
Paranjape A., Sheth R. K., Desjacques V., 2013, MNRAS, 431, 1503
Pillepich A., Porciani C., Reiprich T. H., 2012, MNRAS, 422, 44
Planck Collaboration XX, 2014, A&A, 571, A20
Planck Collaboration XVI, 2014, A&A, 571, A16 (Planck13)
Press W. H., Schechter P., 1974, ApJ, 187, 425
Reed D. S., Smith R. E., Potter D., Schneider A., Stadel J., Moore B., 2013, MNRAS, 431, 1866

Retana-Montenegro E., van Hese E., Gentile G., Baes M., Frutos-Alfaro F., 2012, A&A, 540, A70
Rozo E. et al., 2010, ApJ, 708, 645
Sartoris B. et al., 2015, preprint (arXiv:1505.02165)
Sheth R. K., Tormen G., 1999, MNRAS, 308, 119
Sheth R. K., Tormen G., 2002, MNRAS, 329, 61
Sheth R. K., Mo H. J., Tormen G., 2001, MNRAS, 323, 1
Springel V., 2005, MNRAS, 364, 1105
Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
Springel V. et al., 2005, Nature, 435, 629
Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709
Tormen G., 1998, MNRAS, 297, 648
Tormen G., Moscardini L., Yoshida N., 2004, MNRAS, 350, 1397
Velliscig M. et al., 2015 MNRAS, 454, 3328
Vikhlinin A. et al., 2009, ApJ, 692, 1060
Warren M. S., Abazajian K., Holz D. E., Teodoro L., 2006, ApJ, 646, 881
Watson W. A., Iliev I. T., D'Aloisio A., Knebe A., Shapiro P. R., Yepes G., 2013, MNRAS, 433, 1230
Wu H.-Y., Zentner A. R., Wechsler R. H., 2010, ApJ, 713, 856
Zhao D. H., Jing Y. P., Mo H. J., Bnörner G., 2009, ApJ, 707, 354

## APPENDIX A: RESULTS FOR ELLIPSOIDAL HALOES

We identified haloes using the ellipsoidal halo finder of (Despali et al. 2013) and repeated the analysis described in the main text. Although the best-fitting parameters are systematically different, the results are otherwise consistent with those found for SO haloes: the virial overdensity yields universality, and others do not.

Fig. A1 shows the measured mass function for EO virial haloes and the residuals with respect to the best-fitting relation (calculated at $z=0$), similarly to what was done in Fig. 2. Table A1

<table><caption>**Table A1.** Parameters of the best-fitting mass function for EO haloes at $z=0$.</caption><thead><tr><th>ρ (EO)</th><th>a</th><th>p</th><th>A</th></tr></thead><tbody><tr><td>200 ρ<sub>b</sub></td><td>0.6730 ± 0.004</td><td>0.1783 ± 0.007</td><td>0.4237 ± 0.001</td></tr><tr><td>Δ<sub>vir</sub></td><td>0.7369 ± 0.004</td><td>0.2089 ± 0.007</td><td>0.3894 ± 0.001</td></tr><tr><td>200 ρ<sub>c</sub></td><td>0.8286 ± 0.005</td><td>0.2776 ± 0.009</td><td>0.3335 ± 0.001</td></tr><tr><td>500 ρ<sub>c</sub></td><td>1.0223 ± 0.007</td><td>0.3417 ± 0.009</td><td>0.2672 ± 0.001</td></tr><tr><td>1000 ρ<sub>c</sub></td><td>1.2576 ± 0.009</td><td>0.3803 ± 0.012</td><td>0.2181 ± 0.001</td></tr><tr><td>2000 ρ<sub>c</sub></td><td>1.6088 ± 0.015</td><td>0.3824 ± 0.018</td><td>0.1755 ± 0.001</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>**Figure A1.** Same as Fig. 2, but for the haloes in the EO catalogue.</figcaption></figure>
summarizes how the best-fitting $z = 0$ parameters depend on overdensity threshold. They behave regularly, just as for the SO haloes shown in Fig. 10, and are well modelled by equations (A1) (equivalent to the blue curves in Fig. 10).

For EO haloes, we have obtained analogous relations to those of equation (13):

$$ 
\begin{aligned}
a &= 0.7057 + 0.2125x + 0.3268x^2, \\
p &= 0.2206 + 0.0.1937x - 0.04570x^2, \\
A_0 &= 0.3953 - 0.1768x.
\end{aligned} 
\quad (A1) $$

## APPENDIX B: ISSUES FROM THE INITIAL CONDITIONS

The initial power spectrum $P(k)$ plays an important role when translating from halo mass $M$ to the scaled variable $\nu = \delta_c^2/\sigma^2(M)$, since it determines the value of the mass variance $S(M) = \sigma^2(M)$.

In principle, $\sigma(M)$ should be the same for any simulation run with the same cosmological parameters. In practice, the number of Fourier modes that can be effectively sampled in the initial conditions depends on some computational limits: (i) the box size determines the minimum mode in the power spectrum that can be sampled by a specific simulation; (ii) each of our simulations started from a different random realization of the displacement field. If not properly accounted for, these effects may increase the scatter – and possibly bias – the halo mass function measured in the simulation with respect to theoretical model predictions.

For example, using the theoretical linear power spectrum when computing $\sigma$ yields a precision which cannot be reduced below about 10 per cent. However, future wide field missions (e.g. Euclid; Laureijs et al. 2011) require per cent level precision. As we describe below, to achieve this, we calculate $\sigma(M)$ from the actual realization of the initial power spectrum in each simulation box – scaled using linear theory to $z = 0$ – and not from the theoretical linear power spectrum.

### B1 Power spectrum of each realization

Fig. B1 shows some detailed results on the computational effects introduced in the initial $P(k)$ regarding (i) the random number seed choice and (ii) the box size. The top panel shows $S = \sigma^2(M)$ as a function of mass $M$, for all the simulations of the main set; the black line shows the theoretical power spectrum from CAMB. To compute $S(M)$ for each simulation, we integrated the initial power spectrum from the minimum mode resolvable in each box (see equation 3). The bottom panel shows the relative difference between the measured $S(M)$ and the one calculated using the theoretical $P(k)$. For comparison, the lower panel shows the initial power spectrum for each simulation of our main set and the relative differences with respect to the theoretical one. The use of the actual power spectrum for each simulation allows us to achieve a more precise estimate of $\nu$ for the box, and hence greater precision on the measured mass function.

### B2 The random seed for the initial displacement field

The six simulations from the main set are independent, meaning that they have all been generated from different realizations of the displacement field. However, these simulations also differ in box size, etc. To isolate the effect of the random seed, we re-ran two of the simulations of the secondary set ('uno' and 'wmap7') using two different random seeds. We found few per cent-level differences

<figure><img src="image_20.png" /><figcaption>Figure B1. Top: measured mass variance and residuals with respect to the theoretical calculation. Bottom: initial power spectrum measured from each simulation and relative difference with respect to the theoretical case.</figcaption></figure>

in the high-mass tail, with one of the two seeds generating greater departures from the universality. See e.g. the last blue triangle in the $z=0$ panel of Fig. 2 or the grey and brown ones in the corresponding panel of Fig. 7. This is despite the fact that we use the actual realization of the power spectrum when computing $\sigma(M)$. Thus, even when all other parameters are held fixed, some of the scatter in the measured halo mass function is introduced by the initial seed.

### B3 Mass variance definition

The top part of Fig. B2 compares different definitions for the mass variance $S(M)$. The relative differences are presented with respect to the theoretical prediction – i.e. integrating the theoretical power spectrum(black solid curve). The green line in the top panel (orange
<figure><img src="image_21.png" /><figcaption>**Figure B2.** Top: comparison of the variance as a function of mass using the different prescriptions described in the text to account for the finite box size, for our smallest and largest simulation boxes. Bottom: differences between the measured mass functions when the prescriptions shown in the top panel are implemented, and when the box size is effectively infinite.</figcaption></figure>

line in the middle panel) shows the prediction for $S(M)$ in our smallest (largest) simulation Ada (Flora) where we sum over the modes in the initial conditions, starting from the k mode corresponding to the box size. The blue short dashed line represents one way of accounting for the limited box size of the simulation: we subtract from the theoretical one the mass variance computed using a top-hat filter with scale

$$R_{\text{Box}} = \left( \frac{3 M_{\text{Box}}}{4\pi \rho_{\text{b}}} \right)^{1/3} = L_{\text{Box}} \left( \frac{3}{4\pi} \right)^{1/3}. \qquad (\text{B1})$$

The red long dashed curve shows the case in which we compute $S(M)$ accounting for the fact that the periodic boundary conditions of the simulation mean that the Fourier modes within the box are

<figure><img src="image_22.png" /><figcaption>**Figure B3.** Difference in the measured mass function for different IC methods. For four redshifts (in different colours), we show the logarithmic residuals between the results obtained using the Zel'dovich Approximation (ZA – obtained with N-GenIC) and second-order lagrangian perturbation theory (2LPT – obtained with 2LPTic) to generate the initial conditions. The solid curves are calculated from Reed et al. (2013), where the authors present a fit to the expected ratio between the two methods, which is in good agreement with our measurements.</figcaption></figure>

constrained to give the background density on the scale of the box. This constraint modifies the expression for the variance to

$$s \rightarrow s \left( 1 - \frac{S_x^2}{s S_{\text{Box}}} \right) = s \left( 1 - \frac{S_x^2}{S_{\text{Box}}^2} \frac{S_{\text{Box}}}{s} \right),$$

where

$$S_x(M) = \int \frac{dk}{k} \frac{k^3 P_{\text{lin}}(k, z)}{2\pi^2} W[kR(M)]W[kR_{\text{box}}] \qquad (\text{B2})$$

represents the cross-correlation variance between the lagrangian scale of the halo and that of the simulation box (e.g. Musso & Sheth 2014a,b). For large boxes $S_{\text{Box}} \ll s$ so $S_x/S_{\text{box}} \approx 1$ and the correction to s becomes vanishingly small.

The bottom panels show the $z=0$ mass function residuals in Ada and Flora which result from using different prescriptions for accounting for the finite box size (with respect to the case in which the box size is infinite). For the larger box, Flora, the effect is negligible, but for the smaller box, Ada, it is not. If the mass variance does not properly account for the limited box size, the rescaling to $\nu$ results in non-negligible bias at small masses and large scatter at large masses.

## B4 ZA or 2LPT?

For all our simulations, we used N-GenIC to generate the initial conditions (IC) from a glass distribution; N-GenIC uses first-order perturbation theory, calculated through the ZA. Not including second-order perturbations in the IC has a strong impact on the results obtained at outputs 'close' to the initial conditions Crocce et al. (2006): the differences are larger for simulations with small boxes and/or late starting redshifts. We chose $z=99$ as a starting redshift for all our simulations with a box size larger than $100 h^{-1}$ Mpc; for Ada, which has a box of only $62.5 h^{-1}$ Mpc, we generated the ICs at $z=124$. To test the impact of ICs on our results, we ran an identical copy of Bice and of the two WMAP7 simulations, using exactly the same parameters and seeds, but generating the IC with 2LPTic instead of N-GenIC. The symbols in Fig. B3 show the differences in the measured mass functions; different coloured solid curves show the relation of Reed et al. (2013):

$$\frac{dn_{\text{ZA}}}{dn_{\text{2LPT}}} = e^{-0.12 \frac{a_{\text{f}}}{\sigma_{\text{f}}} v_{\text{1}}^{2.5}}, \qquad (\text{B3})$$
even though it was only calibrated for $1 \lesssim \nu_1 \lesssim 5$. Note that their $\nu$ is defined differently than ours: their $\nu_1 \equiv \sqrt{\nu} = \delta_c(z)/\sigma(M, z)$. At the four redshifts shown, the differences are on average very small – less than 5 per cent at small and intermediate masses and less than 10 per cent in the very high-$\nu$ tail – which lies within the intrinsic scatter (e.g. that due to the random initial seed) seen in the mass function. A small difference between the two definitions can be seen at high redshift (i.e. at $z=5$ in the figure), but this also lies within the intrinsic scatter of the mass function. Since our best fits are calibrated using data at $z \le 1.25$, they are not affected by second-order effects in the initial condition density field. For example, the best-fitting parameters obtained using Reed-rescaled points (all $z$, all cosmologies) are $a = 0.7603$, $p = 0.2549$ and $A_0 = 0.330$ – similar at per mil level to the ones calculated with the original data.

## APPENDIX C: CLUSTER MASS FUNCTION

One of the most important application that necessitate of an accurate and well-calibrated mass function is the study of the observed cluster counts (Vikhlinin et al. 2009; Rozo et al. 2010; Planck Collaboration XX 2014; Sartoris et al. 2015). Degeneracies between fitted parameters mean that the best-fitting parameters to fits restricted to cluster mass haloes ($M \ge 10^{13} h^{-1} M_{\odot}$) may differ from those returned from fitting a larger range of masses (we refer to these as the CMF and HMF, for cluster and halo mass functions, respectively). This appendix discusses the expected level of systematic bias this may induce on cosmological constraints derived from the observed number of clusters per square degree.

The black circles in the top panel of Fig. C1 show the halo mass function extracted from all cosmologies and redshifts from our simulation suite. The main text shows that a good fit to the bins that contain at least 30 haloes (red triangles) is given by a Sheth & Tormen (1999) mass function with the following parameters: $A_0 = 0.3295$, $p = 0.2536$ and $a = 0.7689$. However it is worth mentioning that according to the points density distribution in the considered area, this curve may be less accurate in describing the shape of the most massive haloes: the cluster mass function. To test the difference, and to better describe the objects more massive than $M_{vir} \ge 3 \times 10^{13} h^{-1} M_{\odot}$ that are typically associated with groups and clusters of galaxies, we have performed a fit to the orange crosses only. The best-fitting parameters are $A_0 = 0.8199$, $a = 0.3141$, and $p = 0$, as reported in Table 3. The value $p = 0$ is due to the fact that it is mainly the small masses which determine $p$ (Sheth & Tormen 1999). In the first bottom subpanel, we show the relative difference – in log space – of the two fits for the halo mass function tail. For comparison, the other two subpanels show the relative residuals of the data points with no error bars to the two corresponding fits.

A relative difference in the number of clusters per square degree may appear using the fitting function computed using only the clusters (hereafter cluster mass function) or all haloes over a broader range of masses (halo mass function). To quantify this difference in Fig. C2, we present the relative difference in cluster counts – $M_{vir} > 3 \times 10^{13} M_{\odot} h^{-1}$ – per square degree in the $\Omega_m-\sigma_8$ plane between the CMF and the HMF fits. To compute the count $N$, we have integrated the HMF and CMF fits over $M_{vir} > 3 \times 10^{13} M_{\odot} h^{-1}$ and comoving volume. The figure shows that in our reference Planck13 cosmology the relative difference in the cluster counts between the two fits is ~4 per cent. It increases to ~10 per cent for small values of both $\Omega_m$ and $\sigma_8$. At high $\Omega_m$ and $\sigma_8$ the difference between the HMF and CMF-derived counts is smaller.

<figure><img src="image_23.png" /><figcaption>**Figure C1.** Halo mass function for all considered cosmologies and redshifts extracted from our simulation suite (black circles). Red triangles show bins with at least 30 counts, while the orange crosses also require that the haloes be more massive than $3 \times 10^{13} M_{\odot} h^{-1}$. The solid and dashed curves represent the best fit to the red triangles and orange crosses – accounting for the associated Poisson error bars. The first subpanel shows the relative difference – in log space – of the two fits, while the other two present the relative difference of the data points with respect to their corresponding best-fitting curve.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure C2.** Relative difference in the cluster counts – systems with $M_{vir} > 3 \times 10^{13} M_{\odot} h^{-1}$ per square degree – in the $\Omega_m-\sigma_8$ plane, resulting from the best-fitting parameters obtained for clusters (CMF) and all haloes (HMF) as presented in Fig. 8.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure C3.** Relative cosmological constraints between the CMF and the HMF in the Ω_m-σ_8 plane for systems with M_{vir} > 3 × 10^{13} M_⊙ h^{-1}.</figcaption></figure>

In Fig. C3, we present the relative cosmological constraints for the cluster counts in the Ω_m-σ_8 plane between the HMF and the CMF. The relative difference is well below 10 per cent in the whole range and the diagonal shape shows the degeneration region in the parameter space.

This suggests that, notwithstanding the small difference between the two fits, caution may be necessary when using fits to the halo mass function for precision cosmology. For cluster counts, the CMF reported in the last row of Table 3 may be more appropriate than the HMF, and may yield better than 4 per cent accuracy.

This paper has been typeset from a TeX/LATeX file prepared by the author.