# Euclid preparation. XXIV. Calibration of the halo mass function in $\Lambda(\nu)$CDM cosmologies

T. Castro, A. Fumagalli, R. E. Angulo, S. Bocquet, S. Borgani, C. Carbone, J.

Dakin, K. Dolag, C. Giocoli, P. Monaco, et al.

► To cite this version:

T. Castro, A. Fumagalli, R. E. Angulo, S. Bocquet, S. Borgani, et al., Euclid preparation. XXIV. Calibration of the halo mass function in $\Lambda(\nu)$CDM cosmologies. Astronomy and Astrophysics - A&A, 2023, 671, pp.A100. 10.1051/0004-6361/202244674. obspm-04079759

HAL Id: obspm-04079759

https://hal-obspm.ccsd.cnrs.fr/obspm-04079759

Submitted on 24 Apr 2023

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
# Euclid preparation. XXIV. Calibration of the halo mass function in $\Lambda(\nu)$CDM cosmologies

Euclid Collaboration: T. Castro1,2,3*, A. Fumagalli4,1,2,3, R. E. Angulo5,6, S. Bocquet7, S. Borgani4,3,1,2, C. Carbone8,
J. Dakin9, K. Dolag10,11, C. Giocoli12,13, P. Monaco4,3,1,2, A. Ragagnin14,3,1, A. Saro4,1,3,2, E. Sefusatti3,1,2,
M. Costanzi4,3,1, A. M. C. Le Brun15, P.-S. Corasaniti15, A. Amara16, L. Amendola17, M. Baldi14,18,19, R. Bender20,10,
C. Bodendorf11, E. Branchini21,22, M. Brescia23,24, S. Camera25,26,27, V. Capobianco27, J. Carretero28,29,
M. Castellano30, S. Cavuoti24,31, A. Cimatti32,33, R. Cledassou34,35, G. Congedo36, L. Conversi37,38, Y. Copin39,
L. Corcione27, F. Courbin40, A. Da Silva41,42, H. Degaudenzi43, M. Douspis44, F. Dubath43, C.A.J. Duncan45,46,
X. Dupac37, S. Farrens47, S. Ferriol39, P. Fosalba48,49, M. Frailis1, E. Franceschi18, S. Galeotta1, B. Garilli8, B. Gillis36,
A. Grazian50, F. Grupp20,10, S.V.H. Haugan51, F. Hormuth52, A. Hornstrup53, P. Hudelot54, K. Jahnke55, S. Kermiche56,
T. Kitching57, M. Kunz58, H. Kurki-Suonio59, P.B. Lilje51, I. Lloro60, O. Mansutti1, O. Marggraf61, F. Marulli62,18,19,
M. Meneghetti18,19, E. Merlin30, G. Meylan40, M. Moresco62,18, L. Moscardini62,18,19, E. Munari1, S.M. Niemi63,
C. Padilla28, S. Paltani43, F. Pasian1, K. Pedersen9, V. Pettorino47, S. Pires64, G. Polenta65, M. Poncet34, L. Popa66,
L. Pozzetti18, F. Raison20, R. Rebolo67,68, A. Renzi69,70, J. Rhodes71, G. Riccio24, E. Romelli1, R. Saglia10,20,
D. Sapone72, B. Sartoris10,1, P. Schneider61, G. Seidel55, G. Sirri19, L. Stanco70, P. Tallada Crespi73,29, A.N. Taylor36,
R. Toledo-Moreo74, F. Torradeflot73,29, I. Tutusaus58, E.A. Valentijn75, L. Valenziano18,19, T. Vassallo1, Y. Wang76,
J. Weller10,20, A. Zacchei1,3, G. Zamorani12, S. Andreon77, S. Bardelli18, E. Bozzo43, C. Colodro-Conde67, D. Di
Ferdinando19, M. Farina78, J. Graciá-Carpio20, V. Lindholm59, C.Neissner28, V. Scottez54,79, M. Tenti19, E. Zucca18,
C. Baccigalupi80,3,1,2, A. Balaguera-Antolínez67,68, M. Ballardini81,82,18, F. Bernardeau83, A. Biviano1,3,
A. Blanchard84, A. S. Borlaff85, C. Burigana81,86,87, R. Cabanac84, A. Cappi88,18, C.S. Carvalho89, S. Casas90,
G. Castignani62,18, A. Cooray91, J. Coupon43, H.M. Courtois92, S. Davini93, G. De Lucia1, G. Desprez43, H. Dole44,
J.A. Escartin20, S. Escoffier56, F. Finelli18,87, K. Ganga94, J. Garcia-Bellido95, K. George7, G. Gozalias196,
H. Hildebrandt97, I. Hook98, S. Ilic99,34,84, V. Kansal64, E. Keihanen59, C.C. Kirkpatrick59, A. Loureiro36,100,101,
J. Macias-Perez102, M. Magliocchetti78, R. Maoli103,30, S. Marcin104, M. Martinelli30, N. Martinet105, S. Matthew36,
M. Maturi17,106, R.B. Metcalf62,18, G. Morgante18, S. Nadathur16, A.A. Nucita107,108,109, L. Patrizii19, A. Peel110,
V. Popa66, C. Porciani61, D. Potter111, A. Pourtsidou112,36, M. Pöntinen96, A.G. Sánchez20, Z. Sakr84,17,113,
M. Schirmer55, M. Sereno18,19, A. Spurio Mancini57, R. Teyssier114, J. Valiviita115, A. Veropalumbo116, M. Viel3,2,1,80

*(Affiliations can be found after the references)*

March 17, 2023

## ABSTRACT

Euclid's photometric galaxy cluster survey has the potential to be a very competitive cosmological probe. The main cosmological probe with observations of clusters is their number count, within which the halo mass function (HMF) is a key theoretical quantity. We present a new calibration of the analytic HMF, at the level of accuracy and precision required for the uncertainty in this quantity to be subdominant with respect to other sources of uncertainty in recovering cosmological parameters from Euclid cluster counts. Our model is calibrated against a suite of N-body simulations using a Bayesian approach taking into account systematic errors arising from numerical effects in the simulation. First, we test the convergence of HMF predictions from different N-body codes, by using initial conditions generated with different orders of Lagrangian Perturbation theory, and adopting different simulation box sizes and mass resolution. Then, we quantify the effect of using different halo finder algorithms, and how the resulting differences propagate to the cosmological constraints. In order to trace the violation of universality in the HMF, we also analyse simulations based on initial conditions characterised by scale-free power spectra with different spectral indexes, assuming both Einstein-de Sitter and standard $\Lambda$CDM expansion histories. Based on these results, we construct a fitting function for the HMF that we demonstrate to be sub-percent accurate in reproducing results from 9 different variants of the $\Lambda$CDM model including massive neutrinos cosmologies. The calibration systematic uncertainty is largely sub-dominant with respect to the expected precision of future mass-observation relations; with the only notable exception of the effect due to the halo finder, that could lead to biased cosmological inference.

**Key words.** Galaxy Clusters, Cosmology, Large-Scale Structure

## 1. Introduction

Structure formation in the Universe follows a hierarchical process, with larger structures forming from the collapse of smaller

<footnote>* e-mail: tiago.batalha@inaf.it</footnote>
ones. Galaxy clusters sit at the top of this hierarchy as the most massive virialized objects in the Universe. Their cosmological abundance and evolution make them competitive cosmological probes of both the geometry of our Universe and the growth of density perturbations (Allen et al. 2011; Kravtsov & Borgani 2012, for reviews).

Number-counts experiments represent the main cosmological probe from cluster surveys (Holder et al. 2001; Rozo et al. 2010; Hasselfield et al. 2013; Planck Collaboration XX. 2014; Bocquet et al. 2015; Mantz et al. 2015; Planck Collaboration XXIV. 2016; Bocquet et al. 2019; Abbott et al. 2020; Costanzi et al. 2021; Lesci et al. 2022). The idea behind number counts is remarkable in its simplicity; assuming the halo mass function (HMF) is precisely predicted in a given cosmological model, one can confront the number of observed clusters in a survey with this theoretical prediction and, from this, constrain cosmological parameters.

Press & Schechter (1974) presented the first theoretical model for the HMF based on mapping the peaks of a primordial, Gaussian density field to the number of collapsed structures at low redshift. Later, Bond et al. (1991) introduced the excursion-set formalism that embedded the Press–Schechter derivation in a more rigorous and rich mathematical framework. Despite the success of these models in providing an insightful, qualitative description of the abundance of halos, their agreement with increasingly accurate *N*-body simulations is quantitatively poor (see, e.g. Sheth & Tormen 1999). Later, Maggiore & Riotto (2010a,b,c) extended the excursion-set solution to more realist collapse barriers leading to better agreement with *N*-body simulations (see also Corasaniti & Achitouv 2011a,b).

Several works presented semi-analytical extensions of the Press–Schechter formalism calibrated to reproduce the results of numerical simulations (see, e.g. Sheth & Tormen 1999; Sheth et al. 2001; Jenkins et al. 2001; Warren et al. 2006; Tinker et al. 2008; Watson et al. 2013; Despali et al. 2016). These extensions rely on the feature of the HMF being predominantly universal, that is, being independent of the underlying cosmological model, if expressed as a function of the Gaussian density field statistics.

The departure of the HMF from universality depends on several aspects, and modelling it is a daunting task. Recent works have started using emulators to circumvent the complex analytical modelling to predict the HMF (McClintock et al. 2019; Nishimichi et al. 2019; Bocquet et al. 2020). Although practical and straightforward, this strategy has limitations as emulators are not rigorously supported by a robust underlying model and are known to perform poorly outside the regime in which they have been built. Furthermore, building an emulator does not necessarily lead to a better understanding of nature, lacking a theoretical legacy. Still, it is essential to note that emulators have pushed the precision and accuracy of the HMF predictions to a few percent, outperforming fitting functions, which range in accuracy from 10 to 30 percent.

*N*-body simulations are central to both approaches as they are the only theoretical tool to with which to stringently assess the non-linear regime where galaxy clusters are deep-seated. The ever-increasing requirements for larger and high-resolution simulations pushed the development of fast and efficient gravity solvers, each of them relying on different algorithms with different approximations. Keeping the validity, accuracy, and precision of those approximations under control is utterly necessary in order to guarantee the robustness of the final HMF model (see e.g. Angulo & Hahn 2022, for a detailed discussion of the methodology).

The luminous matter composition of our Universe, although subdominant, is known to impact the formation of structure in a significant way (Cui et al. 2014; Velliscig et al. 2014; Bocquet et al. 2016; Castro et al. 2020). The modelling of the baryonic feedback in hydrodynamical simulations is a controversial subject in vogue. However, at the scales of galaxy clusters, it is well established that baryonic feedback cannot disrupt structures; instead, its net effect is to rearrange the halo composition, changing its mass with respect to the same object simulated using a collisionless scheme. As hydrodynamical simulations are thousands of times more expensive than purely gravitational *N*-body simulations, the standard approach is to characterise the HMF using the latter and model how baryonic physics changes the mass of halos in post-processing (see, e.g. Schneider & Teyssier 2015; Aricò et al. 2021). This method will not be addressed in this paper.

The impact of systematic effects in determining the HMF has already been addressed. For instance, Knebe et al. (2011) and Garcia & Rozo (2019) showed that differences in the halo finder algorithm result in changes in the HMF that can be as large as several percent. Salvati et al. (2020) and Artis et al. (2021) claimed that such discrepancies on the HMF model in future surveys could lead to severe biases in the cosmological constraints, raising awareness of the degree of accuracy required by the quality of next generation data.

In the present work, we aim to calibrate an HMF fitting function with a target accuracy of 1 percent for objects more massive than a few times $10^{13} h^{-1} M_{\odot}$, as demanded by the large number of galaxy clusters expected to be observed by the ESA *Euclid* satellite (see Sartoris et al. 2016). In particular, the wide, *Euclid* survey will cover 15 000 deg$^2$ of the extra-galactic sky (Laureijs et al. 2011) in optical and near-infrared (NIR) bands. By relying on suitable algorithms to identify galaxy clusters as galaxy concentrations in photometric redshift space (Euclid Collaboration: Adam et al. 2019), *Euclid* will provide a survey of optically identified clusters that is unique in terms of the depth and area covered, and that has the potential to provide tight constraints on cosmological models (Sartoris et al. 2016).

To fully exploit the cosmological potential of the *Euclid* Cluster Survey, several observational and theoretical systematic effects must be controlled. As for the latter, we concentrate in this paper on an accurate calibration of the HMF. To this purpose, we follow the approach of Ondaro-Mallea et al. (2021) and explicitly model the non-universality of the HMF. To guarantee the robustness of our results, we start by presenting a detailed study of the numerical systematic errors on *N*-body simulations and how they impact the HMF. We then employ a suite of *N*-body simulations and calibrate a model to predict the abundance of halos as a function of mass, cosmology, and redshift. Lastly, we forecast the impact of numerical and systematic uncertainties on *Euclid*'s number counts analysis.

This paper is organised as follows: in Sect. 2 we present our methodology. The study of the convergence of the simulation setup used in this work is presented in detail in Sect. 3. In Sect. 4, we present our modelling of the HMF. Our results are presented in Sect. 5. Finally, concluding remarks are provided in Sect. 6. Additionally, in Appendixes A, B, and C, we present further numerical convergence tests for our adopted setup, a comparison between the HMF calibrated in different halo finders, and the robustness of our calibration as a function of redshift, respectively.
## 2. Methodology

In this section, we present a short overview of the main concepts of the HMF (Sect. 2.1), a description of the sets of simulations carried out for our analysis (Sect. 2.2), the halo finders adopted (Sect. 2.3), the Bayesian framework used for the HMF calibration (Sect. 2.4), and the pipeline to forecast the effects of uncertainties on the HMF calibration on the cosmological constraints from the *Euclid* photometric survey of galaxy clusters (Sect. 2.5).

### 2.1. The halo mass function

The comoving number density of halos with mass in the range $[M, M + dM]$ is given by:

$$ \frac{dn}{dM} dM = \frac{\rho_m}{M} v f(v) d \ln v. \quad (1) $$

Equation (1) is known as the differential halo mass function, where $\rho_m$ is the comoving cosmic mean matter density, $v$ is the peak height, and the function $v f(v)$ is known as the multiplicity function. The peak height $v$ is a measure of how rare a halo is and is defined as $v = \delta_c / \sigma(M, z)$, where $\delta_c$ is the critical density for spherical halo collapse extrapolated to $z = 0$ (Peebles 2020) and $\sigma^2(M, z)$ is the mass variance at redshift $z$, which can be expressed in terms of the linear matter power spectrum $P_m(k, z)$ as:

$$ \sigma^2(M, z) = \frac{1}{2\pi^2} \int_0^\infty dk k^2 P_m(k, z) W^2(k R), \quad (2) $$

where $R(M) = [3 M / (4 \pi \rho_m)]^{1/3}$ is the Lagrangian radius of a sphere containing the mass $M$, and $W(kR)$ is the top-hat filter in $k$-space.

The multiplicity function is said to be universal if its cosmological dependence comes only through its dependence on the peak height $v$, with the functional form of this dependence being independent of cosmology. While this assumption of HMF universality holds to first approximation, a number of independent analyses based on $N$-body simulations have demonstrated that the HMF systematically deviates from universality at late times (Crocce et al. 2010; Courtin et al. 2011; Watson et al. 2013; Diemer 2020; Ondaro-Mallea et al. 2021).

The non-universality of the HMF has been observed to depend on both the halo definition (Watson et al. 2013; Despali et al. 2016; Diemer 2020; Ondaro-Mallea et al. 2021) and the residual dependence of $\delta_c$ on cosmology. For instance, although several works have used the mass variance $\sigma^2(M, z)$ as a characteristic parameter (Jenkins et al. 2001; Reed et al. 2003; Warren et al. 2006; Tinker et al. 2008; Crocce et al. 2010; Watson et al. 2013; Bocquet et al. 2016; Castro et al. 2020), it has been shown by Courtin et al. (2011) that correctly using the peak height $v(M, z) = \delta_c / \sigma(M, z)$, thus including the cosmological dependence of the linear collapse barrier $\delta_c$, provides more universal results. The effect is more prominent at high masses, which is the regime relevant for cluster cosmology and where we dedicate particular attention in this paper.

In our analysis, halos are defined as spherical overdensities with average enclosed density equal to $\Delta_{vir}(z)$ times the background density, where $\Delta_{vir}(z)$ is the non-linear density contrast of virialized structures as predicted by spherical collapse (see, e.g. Eke et al. 1996; Bryan & Norman 1998). In the following, we parameterize the HMF for a given cosmology at a given redshift according to the fitting function introduced by Bhattacharya

et al. (2011),

$$ v f(v) = A(p, q) \sqrt{\frac{2\alpha v^2}{\pi}} e^{-\alpha v^2/2} \left(1 + \frac{1}{(\alpha v^2)^p}\right) (v \sqrt{a})^{q-1}. \quad (3) $$

The dependence of the fitting parameters $\{a, p, q\}$ on cosmology is presented in Sect. 4. As for $\delta_c$, we use the fitting formula presented by Kitayama & Suto (1996).

According to the Halo Model (see e.g. Cooray & Sheth 2002, and references therein), all matter in the Universe is contained in halos. Thus, the integral of the HMF should be normalised to unity. In the case where $a, p$, and $q$ do not depend on $v$, this condition is satisfied by:

$$ A(p, q) = \left\{ \frac{2^{-1/2-p+q/2}}{\sqrt{\pi}} \left[ 2^p \Gamma\left(\frac{q}{2}\right) + \Gamma\left(-p+\frac{q}{2}\right) \right] \right\}^{-1}, \quad (4) $$

where $\Gamma$ denotes the Gamma function. However, we note that not all HMF models presented in the literature obey this condition. Furthermore, this normalisation is not likely to hold in our Universe (see e.g. Angulo & White 2010). As discussed below, we write the parameters $\{a, p, q\}$ explicitly as functions of the matter power spectrum shape and background evolution; thus, our model violates the condition of normalisation to unity even if it formally uses Eq. (4). In this sense, Eq. (4) is used here only to guarantee an asymptotic behaviour to the HMF when extrapolated to lower masses.

### 2.2. Simulations

Due to its intrinsically non-linear dynamics, the formation of virialized halos can only be studied in detail by resorting to cosmological simulations. The need for an accurate $N$-body simulation to numerically solve the evolution of billions of particles during many dynamical time scales entails the design and development of high performance and scalable algorithms and codes. Each gravity solver algorithm comes with different technical aspects characterising its validity, accuracy, and precision. A comprehensive comparison of the publicly available codes is beyond the scope of this paper (see e.g. Angulo & Hahn 2022, and references therein). Still, understanding and quantifying differences in the HMFs predicted by different, and widely used, $N$-body codes is key to establishing confidence in the convergence and robustness of the $N$-body solution of the HMF.

The HMF analysis that we present in this paper is based on three sets of simulations specifically designed for this purpose and which are summarised in Table 1. The *TEsting simulation SET* (TEASE) is used in Sects. 3 and 4 to quantify the aforementioned impact of numerical systematic effects on the HMF. The set consists of simulations of $500 h^{-1} \text{ Mpc}$ cosmological boxes run with different setups and different codes. We used the codes: Open-GADGET, GADGET-4,<fnref n="1" /> PKDGRAV-3,<fnref n="2" /> CONCEPT,<fnref n="3" /> and RAMSES,<fnref n="4" /> which respectively deploy the following gravity solver algorithms: tree-particle mesh (tree-PM; Xu 1995; Bagla 2002; Springel 2005), fast multipole method-particle mesh (FMM-PM; Springel et al. 2021), (FMM; Potter et al. 2017), particle particle - particle mesh (P³M; Dakin et al. 2022), and adaptive mesh refinement (AMR; Teyssier 2002). Initial conditions have been generated using either the Zeldovich

<footnote n="1">¹ https://www.mpa.mpa-garching.mpg.de/gadget4/</footnote>

<footnote n="2">² https://bitbucket.org/dpotter/pkdgrav3</footnote>

<footnote n="3">³ https://jmd-dk.github.io/concept</footnote>

<footnote n="4">⁴ https://www.ics.uzh.ch/~teyssier/ramses</footnote>
approximation (Zeldovich 1970) at $z = 99$ with MUSIC (Hahn & Abel 2011)<fnref n="5" /> or third-order Lagrangian perturbation theory (3LPT) at $z = 24$ with monofonIC (Michaux et al. 2020).<fnref n="6" /> For this simulation set, the background expansion history and linear matter power spectrum are both kept fixed, while resolution is varied (see Table 1).

The *scAle frEe simulaTions fOr cLuster cOsmoloGY* (AE-TIOLOGY) set of simulations is used in Sect. 4 to model the departure of the HMF from universality. The set consists of simulations with 1024³ particles in boxes of 1000 $h^{-1}$ Mpc on a side. The simulations and the initial conditions (2LPT at $z = 99$) were run with GADGET-4. This set is specifically designed to discriminate between the non-universality arising from the presence of characteristic timescales in the background expansion history, and from the presence of characteristic length scales in the matter power spectrum. For this purpose, the set combines simulations with either power-law or $\Lambda$CDM linear power spectrum on a background that is either Einstein-de Sitter – that is $(\Omega_m, \Omega_\Lambda) = (1, 0)$; EdS hereafter – or $\Lambda$CDM. The self-similarity of scale-free simulations has been extensively used by Leroy et al. (2021), Garrison et al. (2021a), Garrison et al. (2021b), Joyce et al. (2021), and Maleubre et al. (2022) to carry out controlled tests of numerical convergence for several cosmological algorithms. Here, instead, we focus on self-similarity to isolate the effect of the power spectrum shape from the effect of the background evolution on the HMF universality. We note that the idea of carrying out simulations with independent background and power spectrum parameters was also used by Ondaro-Mallea et al. (2021).

Lastly, the *PathfInder Cluster COsmoLOgy* (PICCOLO) set is used in Sect. 5 to calibrate our HMF model. The set comprises 26 cosmological boxes of $2000 h^{-1}$ Mpc in size and $4 \times 1280^3$ particles. These simulations were carried out with Open-GADGET. Initial conditions were generated using monofonIC according to 3LPT at $z = 24$. The set of cosmological parameters for this set are varied by uniformly drawing them from the 95 percent confidence level hyper-volume of the joint SPT and DES cluster abundance constraints presented by Costanzi et al. (2021) (see Table 2). Two realisations were simulated for each cosmology, with the exception of the reference C0 model, for which ten realisations were carried out.

In this paper, we deliberately ignored the effect of massive neutrinos in the simulations. Castorina et al. (2014) showed that the effects of massive neutrinos on the HMF are mostly absorbed into the computation of the mass variance once the cold-dark-matter power spectrum is considered instead of the total matter power (see also Costanzi et al. 2013). We will rely on this approach to extend our modelling to scenarios with massive neutrinos; however, as we recognise the importance of validating this approach ensuring it provides the accuracy required for the application to the *Euclid* cluster survey, in Sect. 5.2 we compare our model to two external sets of simulations that include neutrinos at the particle level.

We emphasise that, in order to guarantee a fair comparison between the different gravity solvers, we proceeded as follows.

- We carried out convergence tests for the adopted setup of Open-GADGET and GADGET-4 (see Sect. 3).

- We used the default setup of the CONCEPT code. The default setup of CONCEPT was extensively tested and compared with other codes by Dakin et al. (2022).

- We used the conservative PKDGRAV-3 example setup (D. Potter, 2021, priv. comm.).

- We used the recommended RAMSES setup suggested by the MUSIC and monofonIC initial conditions generators.

## 2.3. Halo finders

For the purpose of testing the HMF calibration against the choice of the halo finder, we used AHF,<fnref n="7" /> SUBFIND,<fnref n="8" /> VELOCIRaptor,<fnref n="9" /> DENHF, and ROCKSTAR<fnref n="10" /> to extract halo catalogues from the different simulation sets. While all these algorithms are based on a spherical overdensity (SO) method to define halo boundaries, they differ in the method applied to identify the centre from which spheres are grown. For each of these halo finders, we provide a very short description of their main characteristics below, while we refer to the original papers where they have been presented. We further refer to Knebe et al. (2011) for an extensive and detailed comparison of the different halo finders.

- AHF (Gill et al. 2004; Knollmann & Knebe 2009) deploys an AMR algorithm to identify prospective halos centres as density peaks on the end-leaves of the refined tree.

- DENHF (Despali et al. 2016) determines prospective halo centres as peaks of the density field estimated using the inverse of the cubic distance to the n-th closest neighbours.

- SUBFIND (Springel et al. 2001, 2021) determines halo centres by first executing a parallel implementation of the 3D friends-of-friends (FOF, see e.g. Davis et al. 1985) algorithm and then by assigning it to the position of the particle with the lowest potential.

- VELOCIRaptor (Elahi et al. 2019) can operate very similarly to SUBFIND if required. Additionally, the user can choose three group finders: 3DFOF, 6DFOF (see e.g. Diemand et al. 2006), and a 6DFOF with adaptive velocity scale. Lastly, halo centres are defined as the group's centre of mass.

- ROCKSTAR (Behroozi et al. 2013a) starts by partitioning the simulation volume into 3DFOF groups for straightforward parallelization. An adaptive and recursive 6DFOF algorithm is then run for each group, creating a hierarchical set of FOF subgroups. Halo centres are computed by averaging the positions of the particles belonging to the innermost subgroup in the hierarchy. We run the CONSISTENT<fnref n="11" /> algorithm on top of the previously extracted ROCKSTAR catalogues. CONSISTENT was shown in Behroozi et al. (2013b) to improve the consistency of the halo catalogues using a dynamical tracking of halo progenitors.

It is important to note that we employ the same halo mass definition; thus, the algorithms above will differ only in the centring definition and, more importantly, in their hierarchical conditions, that is, the requirements to discriminate between structures and substructures. Unless stated otherwise, we only consider inclusive SO masses, that is, including all particles inside the virial radius regardless of whether they are gravitationally bound.

All halo catalogues have been logarithmically binned according to the number of halo particles to minimise the effect of mass discretization. This effect and the adopted binning are further discussed in Sect. 4.1.

<footnote n="5">⁵ https://bitbucket.org/ohahn/music</footnote>

<footnote n="6">⁶ https://bitbucket.org/ohahn/monofonic</footnote>

<footnote n="7">⁷ http://popia.ft.uam.es/AHF/Download.html</footnote>

<footnote n="8">⁸ See n. 1.</footnote>

<footnote n="9">⁹ https://github.com/pelahi/VELOCIRaptor-STF</footnote>

<footnote n="10">¹⁰ https://bitbucket.org/gfcstanford/rockstar</footnote>

<footnote n="11">¹¹ https://bitbucket.org/pbehroozi/consistent-trees</footnote>
<table><caption>**Table 1.** The suite of simulations. The corresponding cosmology is given in Table 2.</caption><thead><tr><th rowspan="2">Set</th><th rowspan="2">L<sub>box</sub><br/>[h<sup>-1</sup> Mpc]</th><th rowspan="2">N<sub>p</sub></th><th rowspan="2">Background</th><th rowspan="2">P<sub>lin</sub>(k)</th><th colspan="3">Initial Conditions</th><th rowspan="2">Grav. Solver</th></tr><tr><th>Code</th><th>LPT Order</th><th>z</th></tr></thead><tbody><tr><td rowspan="5">TEASE</td><td rowspan="5">500</td><td>256<sup>3</sup></td><td rowspan="5">C0</td><td rowspan="5">ΛCDM</td><td rowspan="2">MUSIC</td><td rowspan="2">Zel.</td><td rowspan="2">99</td><td rowspan="2">tree-PM, FMM-PM<br/>FMM, P<sup>3</sup>M, AMR<sup>†</sup></td></tr><tr><td>512<sup>3</sup></td></tr><tr><td>1024<sup>3</sup></td><td rowspan="3">4 × 160<sup>3</sup></td><td rowspan="3">4 × 320<sup>3</sup></td><td rowspan="3">4 × 640<sup>3</sup></td><td rowspan="3">4 × 1280<sup>3</sup></td></tr><tr><td>4 × 1280<sup>3</sup></td></tr><tr><td>monofonIC</td><td>3LPT</td><td>24</td><td>tree-PM, FMM-PM<br/>FMM, P<sup>3</sup>M</td></tr><tr><td rowspan="2">AETIOLOGY</td><td rowspan="2">1000</td><td rowspan="2">1024<sup>3</sup></td><td>EdS</td><td>Power-law<br/>ΛCDM (C0)</td><td rowspan="2">GADGET-4</td><td rowspan="2">2LPT</td><td rowspan="2">99</td><td rowspan="2">FMM-PM</td></tr><tr><td>C0</td><td>Power-law<br/>ΛCDM</td></tr><tr><td>PICCOLO</td><td>2000</td><td>4 × 1280<sup>3</sup></td><td>C0 – C8</td><td>ΛCDM</td><td>monofonIC</td><td>3LPT</td><td>24</td><td>tree-PM</td></tr></tbody></table>

†Corresponding to the gravity solver algorithms deployed in Open-GADGET, GADGET-4, PKDGRAV-3, CONCEPT, and RAMSES, respectively.

<table><caption>**Table 2.** The cosmological parameters of the simulations presented in Table 1. The parameters have been uniformly drawn from the 95% confidence hyper-volume of the cluster abundances constraints presented in Costanzi et al. (2021).</caption><thead><tr><th>Name</th><th>Ω<sub>m</sub></th><th>h</th><th>Ω<sub>b</sub></th><th>n<sub>s</sub></th><th>σ<sub>8</sub></th></tr></thead><tbody><tr><td>C0</td><td>0.3158</td><td>0.6732</td><td>0.0494</td><td>0.9661</td><td>0.8102</td></tr><tr><td>C1</td><td>0.1986</td><td>0.7267</td><td>0.0389</td><td>0.9775</td><td>0.8590</td></tr><tr><td>C2</td><td>0.1665</td><td>0.7066</td><td>0.0417</td><td>0.9461</td><td>0.8341</td></tr><tr><td>C3</td><td>0.3750</td><td>0.6177</td><td>0.0625</td><td>0.9778</td><td>0.7136</td></tr><tr><td>C4</td><td>0.3673</td><td>0.6353</td><td>0.0519</td><td>0.9998</td><td>0.7121</td></tr><tr><td>C5</td><td>0.1908</td><td>0.6507</td><td>0.0527</td><td>0.9908</td><td>0.8971</td></tr><tr><td>C6</td><td>0.2401</td><td>0.8087</td><td>0.0357</td><td>0.9475</td><td>0.8036</td></tr><tr><td>C7</td><td>0.3020</td><td>0.5514</td><td>0.0674</td><td>0.9545</td><td>0.8163</td></tr><tr><td>C8</td><td>0.4093</td><td>0.7080</td><td>0.0446</td><td>0.9791</td><td>0.7253</td></tr></tbody></table>

## 2.4. HMF calibration: maximum likelihood approach

In the Press-Schechter formalism, halos are assumed to form due to the gravitational collapse of matter over-densities, filtered over a given mass scale $M$, above the linear collapse threshold $\delta_c$. The unconditional distribution of the number of excursions above a given threshold on a Gaussian field follows a Poisson distribution. This motivates us to write the likelihood $\mathcal{L}(N_i|\theta, z)$ for the number of halos $N_i^{\text{sim}}$ with masses $M_{\text{halo}} \in [M_i, M_{i+1})$ in the snapshot at redshift $z$ as:

$$ \ln \mathcal{L}(N_i|\theta, z) = N_i^{\text{sim}} \ln \left( \frac{N_i(\theta, z)}{N_i^{\text{sim}}} \right) - N_i(\theta, z) + N_i^{\text{sim}}, \quad (5) $$

where $N_i(\theta, z)$ is the theoretical expectation of halos computed by integrating Eq. (1) and multiplying it by the volume of the simulated cosmological box. Lastly, $\theta$ is the vector of parameters of Eq. (3) including their dependency on cosmology.

However, numerical systematic effects, such as round-off errors, can affect the distribution of halos introducing further scatter between the predicted HMFs. This problem is more prominent for low-mass halos as the abundance of halos quickly grows with decreasing mass and Poisson errors under-predict the total error budget. In this paper, we use the composite likelihood in-

stead:

$$ \ln \mathcal{L}(N_i|\theta, z) = \begin{cases} N_i^{\text{sim}} \ln \left( \frac{N_i(\theta, z)}{N_i^{\text{sim}}} \right) - N_i(\theta, z) + N_i^{\text{sim}} & N_i^{\text{sim}} \le 25 \\ \frac{1}{2} \ln (2\pi\sigma^2) + \frac{(N_i(\theta, z) - N_i^{\text{sim}})^2}{2\sigma^2} & N_i^{\text{sim}} > 25, \end{cases} \quad (6) $$

where $\sigma = \sqrt{N_i^{\text{sim}} + \sigma_{\text{sys}}^2}$ is the standard deviation resulting from the convolution of the Gaussian approximation to the Poisson distribution with a noise distribution assumed to be Gaussian with zero mean and variance $\sigma_{\text{sys}}^2$. We note that our likelihood presents a discontinuity at $N_i^{\text{sim}} = 25$. The impact of the discontinuity is two-fold. Firstly, the discontinuity amplitude depends on the chosen value for $\sigma_{\text{sys}}$. As we discuss here below, we assume $\sigma_{\text{sys}} \lesssim 1$ percent, thus minimising the impact on the transition regime. Secondly, there is the transition from the purely Poisson error to a symmetric Gaussian approximation. The chosen transition value for $N_i^{\text{sim}}$ guarantees that the transition is again smooth as the Poisson contribution largely dominates the total error budget while negative values for $N_i(\theta, z)$ are very rare ($\gtrsim 5\sigma$) fluctuations of a Gaussian distribution. We also checked a more complex composite likelihood function with a two-sided $\sigma$ following Watson et al. (2013):

$$ \sigma^2 = \left( \sqrt{N_i^{\text{sim}} + 0.25} \pm 0.5 \right)^2 + \sigma_{\text{sys}}^2. \quad (7) $$

The relative difference in the log-likelihood between the two functional forms is below $10^{-5}$ around the best-fit values and provides statistically indistinguishable results. We stick to the functional form presented in Eq. (6) for simplicity.

The final log-likelihood is computed by summing Eq. (6) over all redshifts, mass bins, and simulations. This amounts to assuming that different mass bins at fixed redshift and simulation outputs at different redshifts are independent of each other. However, we note that when using the output of the same simulation at different redshifts, the results are clearly not independent, as the binning in a given redshift will contain the progenitors or descendants of the object in another redshift. To minimise the
impact of the correlation, we follow Bocquet et al. (2016) and use a time-spacing of roughly 1.7 Gyr; such spacing is larger than the characteristic dynamical time of galaxy clusters.

## 2.5. Forecasting Euclid's cluster counts observations

To understand the impact of the HMF systematic errors on cosmological constraints it is important to realistically forecast the cosmological information to be extracted from the Euclid photometric cluster survey. Synthetic cluster abundance data are generated in the following way: we consider a Euclid-like light cone covering 15 000 deg$^2$, with redshift range $z = [0, 2]$ (Laurijs et al. 2011). As discussed in detail by Euclid Collaboration: Adam et al. (2019), clusters from the Euclid Wide Survey (Euclid Collaboration: Scaramella et al. 2022) will be identified as overdensities in photometric redshift space by applying two cluster finders, which have been demonstrated to be the most accurate in terms of completeness and purity among those considered, namely AMICO (Bellagamba et al. 2018) and PZWav (Gonzalez 2014). Once clusters are identified an optical richness is assigned to them.

The abundance of halos is sampled assuming our primary calibration for the HMF presented in Sect. 5. Optical richness $\lambda$ is assigned to the halos according to the richness-mass relation $\langle\lambda|M_{vir}, z\rangle$ (see e.g. Saro et al. 2015):

$$ \langle\ln \lambda|M_{\text{vir}}, z\rangle = \ln A_{\lambda} + B_{\lambda} \ln \left( \frac{M_{\text{vir}}}{3 \times 10^{14} h^{-1} M_{\odot}} \right) + C_{\lambda} \ln \left( \frac{E(z)}{E(z=0.6)} \right), \quad (8) $$

where $E(z)$ is the ratio of the Hubble parameter at redshift $z$ and 0. We assume a richness range $\lambda = [20, 2000]$ and a log-normal scatter given by:

$$ \sigma^2_{\ln \lambda|M_{\text{vir},z}} = D_{\lambda}^2. \quad (9) $$

We use the following fiducial values for the parameters of Eqs. (8) and (9) $A_{\lambda} = 37.8$, $B_{\lambda} = 1.16$, $C_{\lambda} = 0.91$, $D_{\lambda} = 0.15$. These parameter values were determined by converting the richness-mass relation presented by Saro et al. (2015) for $M_{500c}$ (presented in their Table 2) to the virial mass definition, assuming that halos have a Navarro-Frenk-White (NFW) profile (Navarro et al. 1997) and follow the mass-concentration relation given by Diemer & Joyce (2019). The adopted values are in agreement with the results presented by Castignani & Benoist (2016).

Lastly, Poisson and sample variance fluctuations are added through a multivariate Gaussian distribution with amplitude given by the covariance model of Hu & Kravtsov (2003), which was validated by Euclid Collaboration: Fumagalli et al. (2021).

# 3. Setup of N-body simulations

In this section, we define an accurate and precise numerical setup for our primary N-body code, Open-GADGET. We present convergence tests for the adopted configuration and parameter values (Sect. 3.1), and the simulation resolution (Sect. 3.2). Lastly, we compare the convergence of our configuration setup to other N-body solvers (Sect. 3.3).

## 3.1. Parameter and configuration setting

The chosen values for the internal key simulation parameters used within Open-GADGET are presented in Table 3. Param-

<table><caption>**Table 3.** Chosen values for the Open-GADGET runs. Parameter file parameters corresponds to the variables set at run time while configuration parameters are set at compilation time. From top to bottom they control: the time integration accuracy, the maximum time step allowed, the tree opening criterion, the angle threshold for tree opening, and the force accuracy, the maximum distance used to compute short-range forces, and the matching region between short-range and long-range forces.</caption><thead><tr><th>Variable</th><th>Value</th></tr></thead><tbody><tr><td colspan="2">Parameter File</td></tr><tr><td>ErrTolIntAccuracy</td><td>0.05</td></tr><tr><td>MaxSizeTimeStep</td><td>0.05</td></tr><tr><td>TypeOfOpeningCriterion</td><td>1</td></tr><tr><td>ErrTolTheta</td><td>0.4</td></tr><tr><td>ErrTolForceAcc</td><td>0.0025</td></tr><tr><td colspan="2">Configuration File</td></tr><tr><td>RCUT</td><td>6.0</td></tr><tr><td>ASMTH</td><td>1.25</td></tr></tbody></table>

eter file variables are set at run time, controlling time integration accuracy, the maximum time step allowed, the tree opening criterion, the value of the critical angle for tree opening, and the force accuracy. The configuration file parameters are set at compilation time. These latter control the maximum distance used to compute short-range forces (RCUT) and the matching region between short-range and long-range forces (ASMTH). See details of the implementation in Springel (2005) and in the official user guide.<fnref n="12" />

We checked that our parameter set provides sub-Poisson differences relative to higher precision sets (see Appendix A for details), deviating by less than a fraction of a percent in the abundance of halos more massive than $3 \times 10^{13} h^{-1} M_{\odot}$. In this test and the following ones, we use the cumulative HMF as a test for the convergence instead of the differential HMF for three reasons: firstly, the cumulative version presents less noise than the differential HMF; secondly, we do not want to assume a binning to compute the differential HMF at this stage before discussing the impact of binning (see Sect. 4.1 below); lastly, for cluster cosmology, most of the constraining power comes from measuring the total number of objects above a given mass threshold; the cumulative mass function is the limit where the mass-observable relation tends to a Dirac delta function. We assessed the numerical convergence of the results by directly comparing it to more precise setups where each parameter was set to to half the value presented in Table 3.

Testing the convergence of the simulation results to the configuration settings of parameters presented in Table 1 is slightly more delicate than the parameter settings as the configuration variables control the raw structure of the gravity solver algorithm. For instance, if one selects very large values for RCUT in Open-GADGET, the tree-PM algorithm results will converge to the tree. The tree algorithm struggles to provide accurate force calculations if the particle distribution is close to a regular uniform grid, as is commonly the case for initial conditions generated with low-order LPT. In this latter case, tweaking the parameters in order to search for a stable point does not guarantee convergence to accurate results (see e.g. Springel et al. 2021).

Instead, we test the accuracy of our configuration set through a comparison with GADGET-4. The rationale for this ap-

<footnote n="12">¹² https://www.mpa.mpa-garching.mpg.de/gadget/users-guide.pdf</footnote>
proach is that the fifth order FMM-PM algorithm deployed in GADGET-4 delivers more accurate force calculations than the tree-PM algorithm deployed in Open-GADGET at fixed accuracy parameters and has smaller error correlations due to the possibility of randomising the box centre at every domain decomposition step. We assessed the convergence of the results for three different initial conditions: $512^3$ particles displaced from a uniform regular grid according to 3LPT at $z = 24$, the same number of particles using Zeldovich at $z = 99$, and $4 \times 320^3$ particles displaced from a face centred cubic (FCC) grid according to 3LPT at $z = 24$. In all configurations, we verified that the agreement between Open-GADGET and GADGET-4 is better than a fraction of a percent for all masses of interest. However, due to the higher degree of symmetry, the FCC configuration shows even better agreement between the two codes. In this configuration, the standard tree algorithm delivers accurate forces from the simulation start, suppressing the spatial correlations of the force calculation errors.

## 3.2. Resolution convergence

Previous works (see e.g. Joyce et al. 2005; Marcos 2008; Garri son et al. 2016) showed that the mass element discretization and their departure from the fluid limit at initial conditions introduces transient systematic effects on N-body simulations. Michaux et al. (2020) showed that the impact of these transient effects is significantly suppressed if simulations start from initial condi tions generated with higher order Lagrangian perturbation the ory using a grid of elements with more planes of symmetry at the closest redshift prior to shell-crossing.

In Fig. 1, we present the convergence test for the cumulative mass function concerning both transients and modes sampled in the initial conditions. We consider five resolutions, correspond- ing to particle masses: $\{6.50 \times 10^{11}, 8.13 \times 10^{10}, 1.02 \times 10^{10}, 5.08 \times 10^9, 1.27 \times 10^9\} h^{-1} M_{\odot}$; these values for the particle masses cor- respond to $\{4 \times 160^3, 4 \times 320^3, 4 \times 640^3, 1280^3, 4 \times 1280^3\}$ parti- cles in a box of $500 h^{-1} \text{Mpc}$, respectively. The rationale behind such choices for the number of particles is that $4 \times 160^3$ dif- fers by ~ 2 percent from $256^3$, thus one can directly compare with other commonly used particle numbers at similar compu- tational cost. The gravitational softening is set to one-fortieth of the mean inter-particle distance in all simulations. Initial con- ditions are created by monofonIC using 3LPT at $z = 24$ and FCC grid for all simulations but the simulation with $1280^3$ par- ticles, which uses a standard equally spaced grid. The reference in all redshift panels is the simulation run at the highest reso- lution. For easier inspection, we add the moving average over five bins for each case. We note that the lower resolution simu- lation tends to suppress the formation of halos at all masses at $z = 0$, with a more severe effect at low masses. At $z = 0.14$, we observe good agreement with higher resolution simulations for objects more massive than $5 \times 10^{13} h^{-1} M_{\odot}$. For the other red- shifts, the lower resolution simulation presents either deviations or fluctuations comparable to the Poisson noise (shown in red). However, we note that the Poisson noise was computed assum- ing no correlation between the simulations, which is not true as the simulations share the same initial conditions. However, as assessing the correlation between the simulations would require many more realisations, we still decided to present the Poisson estimates, considering that they represent a very conservative es- timation of the true scatter.

From the next simulation in resolution order, we see a sub- percent agreement in the number of objects more massive than $5 \times 10^{13} h^{-1} M_{\odot}$. At the same computational cost, the simulation

with $4 \times 320^3$ (Fig. 1) particles presents a slightly better conver- gence than $512^3$ (Fig. A.3) at $z = 0$. Lastly, comparing the two most costly simulations at $z = 1.98$, we observe that the sim- ulation $4 \times 640^3$ has a more stable agreement with the higher resolution simulation than the simulation with $1280^3$ particles, despite the factor of two increase in the total number of parti- cles; this illustrates one of the advantages of using FCC grids instead of the standard ones for creating initial conditions.

## 3.3. Comparison of different N-body solvers

To gain insight into the impact of the different gravity solvers on the halo statistics, in Fig. 2 we present the matter density contrast of the 2D projection of a zoomed-in volume of $(15.37 h^{-1} \text{Mpc})^3$ around the position of the most massive halo found by AHF. The box size of this volume corresponds to six times the virial ra- dius of the central halo. The corresponding mass of the central halo in each simulation is $\{2.007, 1.904, 2.029, 1.887, 1.975\} \times 10^{15} h^{-1} M_{\odot}$ for Open-GADGET, GADGET-4, PKDGRAV-3, CONCEPT, and RAMSES. Silver and cyan circles denote the virial radius of halos and subhalos identified in this region, re- spectively. For the sake of clarity, a mass threshold of $(10^{12}) \times 10^{13} h^{-1} M_{\odot}$ was imposed to select the (sub)halos presented in Fig. 2. We observe that all N-body codes produce similar distri- butions of the most massive objects, however, due to slight dif- ferences in the evolution, the relation between a large halo and its surrounding depends on the N-body code as some structures are detected as an isolated halo in a simulation but as a subhalo in others.

Figure 2 also shows a stronger code signature on the dis- tribution of substructures identified by AHF. While in this pa- per we are not interested in any detailed analysis of substruc- tures, we trace them here for the sole purpose of comparing dif- ferent N-body solvers in detail. Open-GADGET, GADGET-4, PKDGRAV-3, and CONCEPT produce similar numbers of sub- structures in large objects, but a very heterogeneous spatial dis- tribution of them. On the other hand, we observe that RAMSES produces a smoother mass distribution than the other codes, sig- nificantly reducing the number of detected substructures. Simi- lar results were previously reported by Elahi et al. (2016). The tendency of RAMSES to give a smoother mass distribution is also confirmed by measuring the NFW concentration parameter $c$ of the central object: $c = \{5.587, 5.283, 6.170, 5.007, 4.610\}$ for Open-GADGET, GADGET-4, PKDGRAV-3, CONCEPT, and RAMSES.

The impact of different N-body solvers on the final results is expected to depend on the simulation resolution as critical pa- rameters are usually set as a function of the number of parti- cles; for example the softening length, the maximum refinement level, and the refinement strategy. Figure 3 presents the ratio of the cumulative mass function at $z = 0$ to that observed in Open- GADGET. The three panels correspond to three different mass resolutions of $\{6.35 \times 10^{11}, 7.94 \times 10^{10}, 9.92 \times 10^9\} h^{-1} M_{\odot}$; re- spectively $\{256^3, 512^3, 1024^3\}$ particles in a box of $500 h^{-1} \text{Mpc}$. For reference, we present the 68 percent confidence level for the CONCEPT case assuming that the number of halos observed in the two simulations is described by a Poisson realization of the mass function.

While the differences of the other codes with respect to Open-GADGET are stable and at the few-percent level for all resolutions for halos more massive than $3 \times 10^{13} h^{-1} M_{\odot}$, the dif- ference between RAMSES and Open-GADGET is quite sensi- tive to resolution. At the lowest resolution considered in Fig. 3, RAMSES agrees with the other codes for halos more massive
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Convergence test for the cumulative mass function with respect to the effects of both transients and modes sampled in the initial conditions, at fixed box size. In each panel, we report results for five particle masses: {6.50 × 10¹¹, 8.13 × 10¹⁰, 1.02 × 10¹⁰, 5.08 × 10⁹, 1.27 × 10⁹} h⁻¹ M⊙; respectively {4 × 160³, 4 × 320³, 4 × 640³, 1280³, 4 × 1280³} particles in a box of 500 h⁻¹ Mpc. The rationale for the chosen number of particles is that 4 × 160³ differs by ~ 2 percent from 256³, and therefore one can directly compare the results presented here with Fig. A.3 at similar computational cost. The filled region in grey represents the 1 percent region around unity and the filled region in red marks the 68 percent confidence level assuming Poisson distribution for the abundance of halos in both simulations. Each panel shows the comparison for a different redshift, as labelled.</figcaption></figure>

than ≈ 3 × 10¹⁴ h⁻¹ M⊙ while producing a significantly smaller number of less massive objects. At the highest resolution considered, this difference is partially reduced and RAMSES agrees to better than 1 percent with the other codes for halos more massive than ≈ 4 × 10¹³ h⁻¹ M⊙. The suppression in the halo abundance and the production of smoother density fields are known signatures of the RAMSES AMR gravity solver. As the adaptive refinement cannot be repeated indefinitely as it is bound to stop before producing empty voxels, the total number of particles in the simulation also sets the maximum force resolution achievable in RAMSES. Lastly, from Figs. 1 and 3, we conclude that the HMF convergence with respect to particle mass is achieved first for the other codes before RAMSES.

# 4. Modelling

In this section, we present our modelling for the HMF and the numerical and theoretical systematic effects that influence its assessment, including their dependence on initial conditions (Sect. 4.1) and the impact of the simulated volume (Sect. 4.2). We revisit the implications of assuming different halo definitions (Sect. 4.3); globally, comparing different halo finders (Sect. 4.3.1) and internally to a single halo finder (AHF), using different centring (Sect. 4.3.2) and different halo mass assignment (Sect. 4.3.3). Lastly, we present our modelling for the non-universality of the HMF (Sect. 4.4), modelling it as a func-

tion of the shape of the power spectrum (Sect. 4.4.1) and the background evolution Sect. 4.4.2).

## 4.1. Sensitivity of the HMF to initial conditions

In Fig. 4, we present the ratio of the cumulative mass function computed from simulations started with the same seed, but different LPT orders and redshifts. We consider the 3LPT and Zeldovich approximation and the following starting redshifts: z = 24, z = 49, and z = 99. The rationale for the chosen starting redshifts is two-fold: z = 99 and z = 49 have been extensively used in the literature to start simulations using Zeldovich and high-order LPT, respectively. Furthermore, Michaux et al. (2020) showed that starting the simulation at z = 24 using 3LPT is a good compromise between the convergence of the LPT (see, for instance, their Fig. 4) and the effect of particle discreteness on several summary statistics. While third-order LPT gives percentile accuracy on the cumulative mass function for objects more massive than 5 × 10¹³ h⁻¹ M⊙ independent of the starting redshift considered, setting initial conditions at z = 99 using the Zeldovich approximation suppresses the formation of structures by ≥ 1 percent with respect to 3LPT. Our results are in agreement with previous studies (Crocce et al. 2006; Reed et al. 2013; Michaux et al. 2020). We also note that, for objects less massive than 5 × 10¹³ h⁻¹ M⊙, the 3LPT initial conditions set at z = 24 slightly boosts the formation of structures. The reason for this is
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Matter density contrast in a 2D projection of a zoomed-in volume of $(15.37h^{-1}Mpc)^3$ around the position of the most massive halo found by AHF in the corresponding Open-GADGET simulation. The box size of this volume corresponds to six times the virial radius of the Open-GADGET central halo. Silver and cyan circles denote the virial radius of halos and subhalos identified in this region, respectively. For the sake of clarity, a mass threshold of $(10^{12}) 10^{13}h^{-1}M_{\odot}$ was imposed to select the (sub)halos presented.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Ratio of the cumulative mass function at $z = 0$ to that measured from the Open-GADGET simulation. From top to bottom, the panels correspond to a particle mass resolution of {6.35 × 10¹¹, 7.94 × 10¹⁰, 9.92 × 10⁹} $h^{-1}$ M☉; respectively, {256³, 512³, 1024³} particles in a box of 500 $h^{-1}$ Mpc. The red filled area represents the 68 percent confidence level for the CONCEPT case assuming that the number of halos observed in the two simulations follows an uncorrelated Poisson distribution.</figcaption></figure>

two-fold: firstly, as discussed in Sect. 3, there is a difference in the tree force accuracy calculations when the particle distribution is close to the initial grid. Secondly, shell crossing is known to artificially boost the formation of smaller objects (Power et al. 2016). In summary, although we have not run tests for 2LPT, from Fig. 4 we infer that the configuration used for the AETI-OLOGY set, that is, 2LPT at $z = 99$, provides results that agree to better than ≲ 2% with 3LPT at $z = 24$ as it should range between the green and red lines.

Besides the sensitivity to the LPT order used to generate initial conditions, structure formation is also sensitive to small perturbations in the initial positions of particles, such as those caused by round-off errors. This is due to the intrinsically chaotic dynamics obeyed by the several thousands of particles whose orbits are integrated by an N-body solver during many dynam-

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Ratio of the cumulative mass function measured from simulations started with same seed, but different LPT order (3LPT vs. Zel-dovich approximation) and starting redshifts ($z = 24, z = 49,$ and $z = 99$), as indicated by the different coloured lines in the legend.</figcaption></figure>

ical times as they follow the collapse of a halo. The variation of a simulation result on small perturbations in the initial conditions is dubbed the butterfly effect. Genel et al. (2019) thoroughly discuss this effect and how it is amplified in hydrodynamical simulations by thermal processes. Correspondingly, to assess the dependence of the HMF on small perturbations to the initial conditions in purely collisionless simulations, we ran ten simulations for which the initial positions of the particles were randomly displaced by a single unit in the last significant single-precision digit. For the box size of 500 $h^{-1}$ Mpc, this perturbation corresponds to a random displacement of ≲ 1 pc. We note that, for isolating the effect of perturbing initial conditions from round-off errors due to the use of single precision, those simulations were run using GADGET-4 with long integer (i.e., 64 bit) positions. The fluctuations in the HMF caused by the perturbation in initial positions are due to the increasing sensitivity of the non-linear structure formation to initial conditions. However, such small fluctuations cannot disrupt large groups. Instead, they can cause differences in the history of these objects that grow in time resulting in particles accreted by a given group in one simulation ending in a different group in another. In Fig. 5, we present the distribution of the mass of halos with similar mass matched by their position between different simulations at $z = 0$. In the top panel, we present the distribution of the relative mass difference for halos with masses residing in six intervals, each of 0.3 dex width, between 10¹³ and 10¹⁵ $h^{-1}$ M☉. We observe that the impact on halo masses depends strongly on the object mass; whereas halos with masses of a few times 10¹³ $h^{-1}$ M☉ can have their masses affected by several percent, the impact reduces to ≈ 1 percent for the most massive objects. In the bottom panel, we repeat the former exercise scaling the distribution by the number of particles $N_p$, corresponding to the object mass: $N_p = M_{vir}/m_p$, where $m_p$ is the particle mass. We note that the effect is rather universal if presented in these terms, and agrees with a Gaussian distribution of zero mean and standard deviation $\sigma \approx 0.4 \sqrt{m_p M_{vir}} \equiv 0.4 m_p \sqrt{N_p}$. Thus, the limited precision in the initial conditions propagates to a sub-Poisson fluctuation in the number of particles belonging to a given halo at low redshift. In relative terms, the effect is larger for objects with fewer particles and represents only a subpercent effect for objects with more than 1500 particles.

Figure 6 presents the results of the same analysis, when double precision (i.e., 64 bit floating-point) is used. For double precision the effect is not only strongly suppressed, but is also no
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Distribution of the relative difference between masses of halos at $z=0$ matched between different simulations with initial conditions perturbed by $\lesssim 1$ pc. **Top:** Distribution of the relative mass difference. **Bottom:** Same distribution scaled by the number of particles corresponding to the object mass.</figcaption></figure>

longer universal when scaled by the number of particles. The significant suppression in the scatter of the mass of the objects indicates that double precision should be used to setup initial conditions of cosmological simulations and internally by the gravity solver, assuming one can afford the factor two increase in the memory and storage. As storing double-precision outputs is not the default option for several codes, the propagation of double-precision round-off errors will not be further discussed in this work.

The mass fluctuations presented in Fig. 5 can be strongly amplified in binned statistics depending on the bin width. In Fig. 7, we show the r.m.s. variation in the HMF induced by the noise in the initial conditions, normalised to the expected Poisson noise, as a function of halo mass. Different curves correspond to different binnings of halo masses. As is commonly done for the calibration of the HMF, we considered logarithmically spaced mass bins. The results shown in Fig. 7 were obtained by creating a synthetic halo catalogue with masses distributed according to the HMF presented by Tinker et al. (2008). After that, several halo catalogues were created by perturbing the halo masses according to the distribution presented in Fig. 5. Lastly, the HMF was extracted from these catalogues binning the halos in mass and dividing by the volume times the mass bin width. Although we have to tacitly assume a volume in the previous exercise, we have observed negligible impact of the nominal volume assumed by considering two different volumes: $(2000 h^{-1} Mpc)^3$ and $(500 h^{-1} Mpc)^3$.

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Same as Fig. 5 but with double precision.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** The r.m.s. variation in the HMF induced by the noise in the initial conditions, normalised to the expected Poisson noise, as a function of halo mass. Curves with different colours correspond to different binnings of halo masses, as indicated in the legend.</figcaption></figure>

From Fig. 7 it is clear that the binning has to be carefully chosen to reduce the butterfly impact. Ideally, the bin width should be much larger than the scatter caused by the round-off errors. This leads to the condition:

$$ \delta \ln M \gg 1/\sqrt{N_{\min}}, \qquad (10) $$

where $N_{\min}$ is the number of particles of the smallest object of interest.

## 4.2. Impact of the simulated volume

The effect of the simulation volume on the HMF is two-fold: firstly, the computation of the mass variance $\sigma^2(M)$ presented in
Eq. (2) has to be truncated to the fundamental mode of the cosmological box; secondly, by construction, only a few independent modes are contained in the simulated volume for the first multiples of the fundamental mode, thus introducing an effect of sample variance into the computation of the mass variance at large halo masses.

Quantifying the effect of the simulated volume in the HMF directly with simulations would require a much more extensive set of simulations than the ones used in this paper. Instead, we assess the impact of the simulated volume through its impact on the mass variance calculation and propagate this effect to the HMF, assuming the analytical prescription.

In Fig. 8, we present the impact of the simulation box size $L$ on the HMF for three different cases $L = \{500, 1000, 2000\} h^{-1}$ Mpc shown in red, blue, and green, respectively. Solid lines represent the mean effect due to the truncation of the mass variance integration to the fundamental mode. The corresponding shaded regions correspond to the 68 percent interval due to the sample variance. In the top panel, we present the effect on the calculation of the mass variance itself, while in the bottom we propagate the impact on the mass variance to the differential HMF. The 68 percent regions were determined creating 1000 synthetic realisations of a matter power spectrum assuming the C0 cosmology for each box. The matter power spectrum was sampled between the fundamental and the 4096th modes. Sample variance was added to the power spectrum perturbing it with a Gaussian fluctuation of amplitude $\sigma_p(k)$ given by:

$$ \sigma_p(k) = \sqrt{\frac{1}{N(k)}} P_m(k), \qquad (11) $$

where $N(k) = 2\pi(\delta_k/k)^2$ is the number of independent modes inside the bin with width $\delta_k$ in k-space. Finally, the differential HMF in the bottom panel of Fig. 8 has been calculated using the model of Tinker et al. (2008).

In Fig. 8, we observe that the exponential dependence of the HMF on the mass variance at high-masses significantly amplify the impact of both the truncation and sample-variance impact on the mass variance. The absence of modes larger than $500 h^{-1}$ Mpc causes the suppression of the formation of objects more massive than $8 \times 10^{14} h^{-1} M_{\odot}$. No significant suppression on the formation of objects is observed for the two larger box sizes. On the other hand, the impact of sample variance on the HMF is below one percent for all boxes considered for halos lighter than $3 \times 10^{14} h^{-1} M_{\odot}$. For the $1000 h^{-1}$ Mpc and $2000 h^{-1}$ Mpc box sizes, the 1 percent threshold is exceeded for halos more massive than $8 \times 10^{14} h^{-1} M_{\odot}$ and $2 \times 10^{15} h^{-1} M_{\odot}$, respectively.

One way of taking into account the effect of the box size on the calibration of the HMF is to calculate the mass variance, using the matter power spectrum computed from the initial conditions (see e.g. Despali et al. 2016). This approach presents a few challenges: for instance, the computation of the power spectrum at initial conditions might be affected by the mesh used to create the initial displacement field – assuming one has not used a ‘glass’<fnref n="13" /> unstructured mesh from which to displace particles. Also, one has to choose the k-binning for the computation of the matter power spectrum wisely, as an overly thin shell would produce a noisy measurement, while an excessively broad shell would affect the mass variance integral accuracy. We instead advocate that to keep the impact of sample variance to a minimum,

<figure><img src="image_9.png" /><figcaption>Fig. 8. Impact of the simulation box size $L$ on the mass variance (top panel) and the differential HMF (bottom panel) for three different cases $L = \{500, 1000, 2000\} h^{-1}$ Mpc depicted in red, blue, and green, respectively. Solid lines represent the mean effect due to the truncation of the mass variance integration to the fundamental mode. The corresponding shaded regions correspond to the 68 percent interval due to the sample variance.</figcaption></figure>

one should use boxes larger than $1000 h^{-1}$ Mpc and produce initial conditions with fixed amplitudes of the initial conditions random field to the desired theoretical power spectrum, as presented in Angulo & Pontzen (2016). Following this approach one circumvents the above-mentioned challenges as the realised and theoretical power spectrum match exactly.

### 4.3. Impact of the halo definition

#### 4.3.1. Sensitivity to different halo finders

For a visual inspection of the effect of different halo finders on the HMF, Fig. 9 shows a comparison between AHF and SUBFIND halos identified within a ($26.65 h^{-1}$ Mpc)$^3$ volume extracted from a $500 h^{-1}$ Mpc box with $512^3$ particles started at $z = 99$ using the Zeldovich approximation. The mass of the largest halo located at the centre is $2 \times 10^{15} h^{-1} M_{\odot}$ while the smallest halo marked in each panel is $10^{13} h^{-1} M_{\odot}$. Whereas AHF and SUBFIND both find the same large groups, for smaller groups, we notice a non-negligible suppression on the number of objects, with SUBFIND tending to group together smaller objects into a larger one. The effect is more evident along the stream on the centre left of the larger object presented in Fig. 9. The tendency of pure FOF-based methods to merger smaller, dynamically distinct objects along tidal streams, building ‘particle bridges’ between structures, is well known (see e.g. Knebe et al. 2011, and references therein).

<footnote n="13">¹³ Glass pre-initial conditions are random samples that are dynamically evolved assuming a gravitational interaction with reversed sign until they settle in a quasi-equilibrium state.</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Comparison between AHF and SUBFIND halos extracted from the same (26.65 $h^{-1}$ Mpc)$^3$ volume. The largest and smallest halo masses depicted are $2 \times 10^{15} h^{-1}$ $M_{\odot}$ and $10^{13} h^{-1}$ $M_{\odot}$, respectively. The cyan circles denote the virial radius of halos identified in this region.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Ratio of the cumulative mass function extracted from AHF, SUBFIND, VELOCIRaptor, and DENHF to that extracted from ROCKSTAR. The filled regions in red and purple represent the 68% confidence regions for the AHF and VELOCIRaptor (6DFOF Adaptive), respectively, assuming that the number of objects in each catalogue follows a Poisson distribution.</figcaption></figure>

Figure 10 presents the ratio of the cumulative mass functions extracted from AHF, SUBFIND, VELOCIRaptor, and DENHF to that extracted from ROCKSTAR. The filled regions in red and purple represent the 68 percent confidence level region, for AHF and VELOCIRaptor (6DFOF Adaptive), respectively, assuming that the number of objects in each catalogue follows a Poisson distribution. Figure 10 clearly shows a separation of the algorithms considered here into two groups, with the 3DFOF and non-adaptive 6DFOF suppressing the number of halos less massive than $10^{14} h^{-1} M_{\odot}$.

We caution that it is not our goal here to model the specifics of each halo finder. The reason is two-fold; firstly, this is a complex task, as many parameters control the different algorithms. Secondly, investigating how the results from different finders compare to each other when changing the respective parameters, even in great detail, would not address the fundamental question of how does the definition of a halo adopted in simulations

compares with the observed clusters. The analysis presented is simply designed to deliver a flexible model for the HMF that can accommodate different halo definitions. With such a model calibrated against different halo catalogues, we will assess the impact of different definitions of halos in the HMF calibration on Euclid's cluster counts analysis.

### 4.3.2. Impact of the halo centre

We assessed the impact of the choice for the halo centre on the cumulative mass function using AHF on one of the TEASE simulations. AHF allows the user to choose the prospective halo centre alternatively as the geometrical centre of the refinement patch, the cell with the lowest potential, the cell with the highest density, or the centre of mass of the particles inside the refinement patch. The latter is our default choice, following the official AHF documentation.

We verified that the AHF cumulative mass function has a percent-level robustness to the choice of the halo centre. We remind that VELOCIRaptor (based on a 3DFOF) and SUBFIND differ only for the choice of the halo centre. Nevertheless, they differ from each other (see Fig. 10) by an amount that is larger than the differences between different halo centres in the AHF. This is, again, due to particle bridges that connect two dynamically distinct objects, which causes a stronger impact on the halo centre choice.

### 4.3.3. Total halo mass versus bound mass

All halo finders considered in this paper include all particles within a given radius, when computing the spherical overdensity masses, regardless of whether the particles are bound to a halo or not. In order to test the effect of this assumption, in Fig. 11, we present the difference in the cumulative HMF between including all particles inside the spherical overdensity and only the contribution from bound particles. AHF determines as bound the particles with velocities smaller than the local escape velocity multiplied by a constant parameter *VescTune* which we set to unity.
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Difference in the cumulative mass function when all particles inside the spherical overdensity are considered versus just the bound particles.</figcaption></figure>

The simulation used for this test is the same as the one used in Sect. 4.3.2. Not surprisingly, assigning to halos only bound particles reduces the cumulative mass function by about 5 per- cent. That is due to the fact that the gravitationally bound mass within a fixed halo radius is by construction smaller than the to- tal mass within the same radius. It is important to note that both are valid definitions to weigh a halo. The most adequate choice will depend on the method used to measure cluster masses from observations. Arguably, if one works with, for instance, masses obtained from the gravitational lensing, the total mass definition should be adopted, as this is the one that contributes to the lens- ing signal. In the remainder of this paper, we refer only to the total mass definition.

## 4.4. Non-universality of the virial mass function

Despali et al. (2016); Diemer (2020); Ondaro-Mallea et al. (2021) showed that the HMF preserves most of its universality when described as a function of the virial mass, as predicted by the spherical collapse model (Eke et al. 1996; Bryan & Norman 1998). Still, departure from universality can reach up to 20 per- cent for the high-end tail of the mass function (see e.g. Diemer 2020).

In the following sections, and for the specific purpose of trac- ing the origin of any departure from universality, we use the AE- TIOLOGY set of simulations to model the dependence of the virial HMF both as a function of the shape of the matter power spectrum and of the background evolution. Unless stated oth- erwise, the simulations used here were run using GADGET-4 with initial conditions generated at $z = 99$ according to 2LPT. Halo catalogues were extracted on the fly with the GADGET- 4SUBFIND implementation. Halos were binned according to their mass using 50 log-spaced intervals in the number of par- ticles.

For the calibration of the HMF, we impose a mass cut that corresponds to a minimum number of 300 particles so as to minimise the noise in the identification of low-mass halos (see e.g. Leroy et al. 2021). We assume the likelihood presented in Sect. 2.4 with $\sigma_{\text{sys}} = 0.01 N_i^{\text{sim}}$. This choice for the systematic error means that the relative error has a floor of 1 percent in the total error budget, thus avoiding over-weighting mass bins with many halos, for which the purely Poisson error under-predicts the true data variance, as discussed in Sects. 4.1 and 4.2. Its value was chosen so that the best-fit neither over- nor under-

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** The 68% and 95% confidence level contours of the calibration carried out independently in several simulations assuming a power-law linear power spectrum in an EdS background, i.e., $P_m(k) \propto k^{n_s}$ and $\Omega_m = 1$.</figcaption></figure>

fits the simulations over the entire redshift range considered, that is, keeping the fit-quality constant. We use the snapshots at $z = \{2.00, 1.25, 0.90, 0.52, 0.29, 0.14, 0.0\}$.

### 4.4.1. Dependence on the power spectrum slope

In order to understand the dependence of the HMF parameters on the power spectrum shape, in Fig. 12, we present the 68 and 95 percent confidence level contours of the calibration carried out independently in several simulations assuming a power-law linear power spectrum in an EdS background, that is, $P_m(k) \propto k^{n_s}$ and $(\Omega_m, \Omega_\Lambda) = (1, 0)$.

The first interesting result is that, even for self-similar cos- mologies, the HMF is not universal against changing the spectral index. While the parameters {$p, q$} of Eq. (3) seem to have a lin- ear dependence on $n_s$, the parameter $a$ exhibits a non-monotonic dependence with a local minimum around $n_s = -1.75$. The over- all simple and well-behaved dependence of the parameters of the HMF on the spectral index over a range of values covering all the relevant regimes for structure formation motivates our approach to select a flexible fitting function that precisely accommodates the HMF shape on a self-similar cosmology. At the same time, the non-universality is modelled through the explicit dependence of the HMF parameters on cosmology.

### 4.4.2. Dependence on the background evolution

In order to understand the dependence of the HMF parameters on the background evolution we follow an approach similar to that presented in Ondaro-Mallea et al. (2021). In Fig. 13, we present the 68 and 95 percent confidence level contours of the calibration carried out independently in simulations with $n_s = \{-2.0, -2.5\}$ and background evolution given by either a cosmology in agree- ment with the Planck 2018 (P18) results for the cosmological parameters (Planck Collaboration VI. 2020) or EdS. For the spe-
<figure><img src="image_14.png" /><figcaption>Fig. 13. The 68% and 95% confidence level contours of the calibration carried out independently in simulations with $n_s = \{-2.0, -2.5\}$ and background evolution given by either a cosmology in agreement with Planck Collaboration VI. (2020) or EdS.</figcaption></figure>

cific purpose of this exercise, the chosen values for the spectral index are not important, and other combinations produce similar trends; still, we choose $n_s = \{-2.0, -2.5\}$ as this range corresponds roughly to the slope of the ΛCDM matter power spectrum on cluster scales.

At fixed $n_s$, we note that $p$ is consistent with being independent of the background evolution. On the other hand, the parameters $\{a, q\}$ consistently show a dependence on the chosen background, with departures from the EdS scenario providing smaller values of such parameters.

### 4.4.3. The HMF model

From the above results, we define the following model to capture the dependence of the multiplicity function on redshift and power spectrum shape as:

$$a = a_R \Omega_m^a(z) \qquad (12)$$

$$p = p_1 + p_2 \left( \frac{d \ln \sigma}{d \ln R} + 0.5 \right) \qquad (13)$$

$$q = q_R \Omega_m^q(z). \qquad (14)$$

where:

$$a_R = a_1 + a_2 \left( \frac{d \ln \sigma}{d \ln R} + 0.6125 \right)^2 \qquad (15)$$

$$q_R = q_1 + q_2 \left( \frac{d \ln \sigma}{d \ln R} + 0.5 \right). \qquad (16)$$

The growth rate, that is, $d \ln G/d \ln(1+z)^{-1}$ where $G(z)$ is the linear growth of density perturbations, was used by Ondaro-Mallea et al. (2021) to characterise the non-universality of the HMF instead of $\Omega_m(z)$. For the cosmological models considered here this

is well approximated by (see e.g. Amendola & Tsujikawa 2010):

$$\frac{d \ln G}{d \ln(1+z)^{-1}} = \Omega_m^\gamma(z), \qquad (17)$$

with $\gamma = 0.55$. Therefore, Eqs. (12) to (16) produce exactly the same results for the cosmological models studied here once one substitutes:

$$\Omega_m(z) \rightarrow \left( \frac{d \ln G}{d \ln(1+z)^{-1}} \right)^{1/\gamma}. \qquad (18)$$

Whether such a substitution leads to universal extensions of the model to cosmologies with growth histories given by modified gravity theories is left for further investigation. We concentrate here on the description as a function of $\Omega_m$, as it is more straightforward to compute and produces the same results.

Regarding the characterisation of the dependence of the multiplicity function on the shape of the power spectrum, the first obvious choice is to use its logarithmic slope $d \ln P_m(k)/d \ln k$ which reduces to the spectral index $n_s$ for power-law cosmologies. However, in more realistic cosmologies, the logarithmic slope of the power spectrum contains fluctuations due to the baryonic acoustic oscillations (BAOs) for the characteristic length of the most massive halos. One can circumvent this problem by smoothing the BAO before the computation of the slope, as was done in Diemer & Kravtsov (2015). Conversely, we propose the logarithmic slope of the mass variance $d \ln \sigma(R)/d \ln R$. For the power-law cosmologies this reduces to:

$$\frac{d \ln \sigma}{d \ln R} = - \frac{(n_s + 3)}{2}. \qquad (19)$$

## 5. Results

In this section, we present the calibration of the HMF model (Sect. 5.1), compare it with previous works (Sect. 5.3), and forecast the impact of numerical and statistical uncertainties related to our calibration on Euclid's cluster counts analysis (Sect. 5.4).

### 5.1. Calibration of the HMF

First, we provide a short recap of the rationale behind the setup of the PICCOLO simulations. Before presenting the calibration of the HMF model validated in the AETIOLOGY suite, we present a summary of the numerical and theoretical systematic effects on the HMF and how the PICCOLO suite was designed to address them. The PICCOLO simulations were run with OpenGADGET set according to the parameters presented in Table 1. This choice of the parameters is shown in Sect. 3.1 and 3.2 to produce results of subpercent accuracy and the robustness of the results is assessed by code-comparison in Sect. 3.3. This suite of simulations uses FCC grids as pre-initial conditions to mitigate the impact of transients due to fluid discretization (see Michaux et al. 2020). The FCC pre-initial conditions and initial displacements computed with 3LPT at $z = 24$ also reduce the impact of correlated errors on the force computations at early redshifts. We refer to Fig. A.2, where we compare Open-GADGET and GADGET-4, and to the discussion of force accuracy presented by Springel et al. (2021). To mitigate the effect of round-off errors on the HMF, we use the binning condition presented at the end of Sect. 4.1. To minimise the impact of sample-variance, box sizes are $2000 h^{-1}$ Mpc and initial conditions have been created with fixed amplitudes of initial density perturbations (Angulo & Pontzen 2016) as discussed in Sect. 4.2.
<table><caption>**Table 4.** Parameters of the multiplicity function $f(v, z)$ presented in Eq. 3 for the PICCOLO set of simulations.</caption><thead><tr><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>a<sub>z</sub></th><th>p<sub>1</sub></th><th>p<sub>2</sub></th><th>q<sub>1</sub></th><th>q<sub>2</sub></th><th>q<sub>z</sub></th></tr></thead><tbody><tr><td colspan="8">ROCKSTAR</td></tr><tr><td>0.7962</td><td>0.1449</td><td>-0.0658</td><td>-0.5612</td><td>-0.4743</td><td>0.3688</td><td>-0.2804</td><td>0.0251</td></tr><tr><td colspan="8">AHF</td></tr><tr><td>0.7937</td><td>0.1119</td><td>-0.0693</td><td>-0.5689</td><td>-0.4522</td><td>0.3652</td><td>-0.2628</td><td>0.0376</td></tr><tr><td colspan="8">SUBFIND</td></tr><tr><td>0.7953</td><td>0.1667</td><td>-0.0642</td><td>-0.6265</td><td>-0.4907</td><td>0.3215</td><td>-0.2993</td><td>0.0330</td></tr><tr><td colspan="8">VELOCIRaptor (adaptive 6DFOF)</td></tr><tr><td>0.7987</td><td>0.1227</td><td>-0.0523</td><td>-0.5912</td><td>-0.4088</td><td>0.3634</td><td>-0.2732</td><td>0.0715</td></tr></tbody></table>

In order to calibrate the parameters in Eqs. (12) to (16), we use the PICCOLO set of simulations presented in Table 1. The halo catalogues were binned in 50 logarithmically spaced bins in number of particles, corresponding to roughly $\delta \ln M \approx 0.15$. We assumed the likelihood presented in Sect. 2.4 with $\sigma_{\text{sys}} = 0.005 N_i^{\text{sim}}$. As for the AETIOLOGY fit, its value was chosen such that the best-fit neither over- nor under-fits the simulations and this value is in agreement with the expected non-Poisson scatter caused by round-off errors discussed in Sect. 4.1 for the most populated bin considered in the calibration. Its impact is rapidly diluted as the Poisson contribution grows quickly with mass and redshift. Lastly, we imposed a mass cut below which a halo is not considered as resolved, which corresponds to a minimum of 300 particles. Finally, we analyzed the outputs of simulations at redshifts $z = \{2.00, 1.25, 0.90, 0.52, 0.29, 0.14, 0.0\}$.

In Table 4, we present the best fitting parameters of our calibration for our primary halo finder (ROCKSTAR) and three other auxiliary halo finders (AHF, SUBFIND, and VELOCIRaptor). Details of the difference between the calibrations for different halo finders are further discussed in Appendix B. The AHF, SUBFIND, and VELOCIRaptor fits were performed using one simulation for each cosmology presented in Table 2. For the ROCKSTAR fit, we used two simulations for each cosmology including C0 with which we carried out ten realisations. We decided not use the remaining eight C0 simulations for calibration in order to avoid over-weighting this cosmology in our derivation of the HMF fit; we only use them to reduce the Poisson fluctuations of rare halos, which allows us to assess our uncertainty in the calibration at this regime.

In Fig. 14, we present the HMF for ROCKSTAR halos and the prediction of the respective best fit from the PICCOLO simulations at $z = 0$. The vertical dotted line represents the cut in mass corresponding to the cut in the minimum number of 300 particles for the C0 runs. Figure 15 presents the ratio of ROCKSTARbest-fit to the mean abundance of halos extracted from the simulations at $z = \{0.0, 0.5, 1.0, 2.0\}$. As in Fig. 14, the vertical dotted line represents the cut in mass corresponding to the cut in the minimum number of 300 particles for the C0 runs. The regions in grey represent the relative 1 percent and 2.5 percent regions. For the sake of illustration, we present the Poisson error bars corresponding to the C0 and C3 (our cosmology with fewer halos) – the former counts with ten realisations while the latter counts with two.

Quite remarkably, our model shows a percent-level agreement with the simulations over the entire mass and redshift ranges considered in the calibration despite the much larger difference between the cosmological models presented in Fig. 14. This result demonstrates that our model HMF is capable of accurately describing the cosmology dependence of the deviations

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** HMF at $z = 0$ for the PICCOLO cosmologies and the ROCKSTARbest-fit prediction. The vertical dotted line represents the cut in mass corresponding to the minimum number of 300 particles for the C0 runs. Error bars are shown assuming a Poisson distribution of the expected number of objects.</figcaption></figure>

from universality, for a given set of simulations and for a given choice of the halo finder. In Appendix C, we present further tests of the robustness of our calibration as a function of redshift.

## 5.2. Cosmologies with massive neutrinos

The response of the HMF to neutrinos $\mathcal{R}$ is defined as

$$ \mathcal{R}(x) = \frac{dn/dm_{\sum M_v=x}}{dn/dm_{\sum M_v=0}}, \qquad (20) $$

where $\sum M_v$ is the sum of the masses of the three neutrino species.

In Figs. 16 and 17, we present the response of the HMF to massive neutrinos and the accuracy of the model presented in Castorina et al. (2014) in order to take this latter into account. Briefly, the model assumes that neutrinos impact the HMF passively, that is, through its effect on the background evolution. Operationally, implementation of the model only requires that we replace the total matter power spectrum $P_m$ in the computation of the mass variance (Eq. 2) with the cold dark matter plus baryons power spectrum and that we ignore the contribution of neutrinos to the mass of the Lagrangian patch corresponding to a given halo. To assess the accuracy of the model, we compare the response $\mathcal{R}$ predicted by the model and from simulations. The rationale for using $\mathcal{R}$ to assess the accuracy of the model is to mitigate the impact of systematic effects due to differences in the setup between the external simulations and the ones used for the calibration of the HMF in this work.

We used two independent external simulations: the DEMNUni (see Carbone et al. 2016, for details of the simulation)<fnref n="14" /> and the Open-GADGET subset of the Euclid’s neutrino code comparison simulations (Euclid Collaboration: Adamek et al. 2022). Halos were extracted in both sets using SUBFIND. Thus, the results presented in this subsection rely on our SUBFIND calibration presented in Table 4.

The DEMNUni set consists of three simulations of 2000 $h^{-1}$ Mpc boxes assuming a cosmology in agreement with Planck 2013 results for the cosmological parameters (Planck Collaboration XVI. 2014).<fnref n="15" /> The reference simulation considers massless

<footnote n="14">¹⁴ https://www.researchgate.net/project/DEMN-Universe-DEMNUni</footnote>

<footnote n="15">¹⁵ https://github.com/jmd-dk/nucodecomp-data</footnote>
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Ratio of ROCKSTARbest-fit to the mean abundance of halos extracted from the simulations at $z = \{0, 0.5, 1.0, 2.0\}$. The vertical dotted line represents the cut in mass corresponding to the minimum number of 300 particles for the C0 runs. The regions in grey represent the relative 1% and 2.5% regions. The Poisson error bars correspond to the C0 and C3 (our cosmology with fewer halos), which count with ten and two realisations, respectively. We only show mass bins with more than 50 halos for improved readability.</figcaption></figure>

neutrinos, while the other two simulations use the particle-based implementation of neutrinos, assuming a total neutrino masses of 0.16 and 0.32 eV, respectively. The reference simulations include 2048³ dark matter particles, while the latter includes the same number of dark matter particles and as many 2048³ neutrino particles.

The Open-GADGET subset of the *Euclid*’s neutrino code comparison shares the same implementation, a similar mass resolution, and a similar cosmological parameters as the DEMUni but in $1000 h^{-1}$ Mpc boxes. Instead, the total neutrino masses simulated are 0.0, 0.15, 0.30, and 0.60 eV. These simulations have been extensively compared with different neutrino implementations. The agreement for the power spectrum, bispectrum, and HMF was observed to be better than 1 percent for the range of interest.

In Fig. 16, we observe the well-established suppression in the cluster abundance caused by massive neutrinos and its dependence on the total neutrino mass. In Fig. 17, we observe that, for total neutrino masses smaller than 0.32 eV and for both simulations, the model of Castorina et al. (2014) agrees with simulations to better than 1 percent (dark shaded area) over the entire mass range of validity of our calibration; as in Fig. 15, the dotted line represents the 300-particle mass cut used for the calibration. We only note that at lower masses, below ~ $3 \times 10^{13} h^{-1} M_{\odot}$ the model HMF tends to underestimate the effect of massive neutrinos, an effect that increases with neutrino masses. While we are not interested in the HMF calibration in this mass range, we ascribe this difference to the effect of neutrino free-streaming on the growth of CDM perturbations, which cannot be captured by simply ignoring the neutrino component in the linear matter power spectrum.

For the most massive neutrinos considered, we only rely on the simulations carried out for the code comparison. Due to its smaller volume, larger fluctuations are observed with respect to the DEMUni, but at the mass threshold of validity of our calibration, we nevertheless observe that the model agrees with the simulation within a few percent, and for lower masses the model starts to over-predict the impact of neutrinos. The agreement for masses larger than $4 \times 10^{13} h^{-1} M_{\odot}$ is again consistent with 1 percent scatter.

### 5.3. Comparison to other works

In Fig. 18, we compare our HMF model to the works of Tinker et al. (2008), Despali et al. (2016), and Ondaro-Mallea et al. (2021). We use our ROCKSTAR calibration as a reference and present our SUBFIND calibration as our SUBFIND setup matches those used by Ondaro-Mallea et al. (2021), making the comparison easier. We note that, in Fig. 10, we present a comparison of the cumulative mass function produced with ROCKSTAR with that produced by DENHF. The latter, used by Despali et al. (2016), differs from the former by roughly 1 percent. Lastly, Behroozi et al. (2013a) compared ROCKSTAR catalogues with the prediction of Tinker et al. (2008) and presented an agreement within 5 percent at $z = 0$. We also present the mean of the multiplicity function $\nu f(\nu)$ measured from the ten PICCOLO C0 runs. The grey areas depict the 1 and 5 percent regions. To embed the significant differences with respect to the other works, we adopted a symmetric log scale on the y-axis, where the region between the dotted lines is presented in a linear scale.

As already shown in Fig. 15, our model accurately reproduces results from simulations over fairly large ranges of masses and redshifts. Globally, the differences are larger at both large
<figure><img src="image_17.png" /><figcaption>Fig. 16. Impact of massive neutrinos on the HMF. We present the ratio of the HMF with massive neutrinos to the corresponding massless neutrinos counterpart as observed in two independent external simulations DEMNUni (top) and the Open-GADGET subset of Euclid's neutrino code comparison simulations (bottom). See the definition of $R$ in Eq. (20).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Accuracy of the model presented in Castorina et al. (2014) in accounting for the impact of massive neutrinos on the HMF. We compare the ratio of the HMF with massive neutrinos to the corresponding massless neutrinos counterpart, assuming the model of Castorina et al. (2014) divided by the same quantity as observed in two independent external simulations: DEMNUni (top) and the Open-GADGET subset of Euclid's neutrino code comparison simulations (bottom). See the definition of $R$ in Eq. (20).</figcaption></figure>

masses and high redshifts, where the statistics are poorer. At $z=0$, the model of Tinker et al. (2008) differs from ours by a maximum of 3 percent for halos less massive than $10^{15} h^{-1} M_{\odot}$; it crosses the 5 percent threshold at $2 \times 10^{15} h^{-1} M_{\odot}$, and beyond that it deviates from our model, predicting significantly fewer halos. The model of Despali et al. (2016) crosses the 5 percent threshold over a narrower mass range; it deviates from our model by more than 5 percent beyond $4 \times 10^{14} h^{-1} M_{\odot}$, and over-predicts the number of halos more massive than $10^{15} h^{-1} M_{\odot}$ by more than 20 percent. Comparing the model by Ondaro-Mallea et al. (2021) to our HMF calibration based on SUBFIND we note that they start differing by more than 5 percent for halos more massive than $7 \times 10^{14} h^{-1} M_{\odot}$, above which the model of these latter starts to follow the same trend as the model of Despali et al. (2016). The picture at $z=1$ and $z=2$ are similar; at smaller masses, all models tend to predict fewer halos than observed in our simulations and predicted by our model. For Despali et al. (2016) and Ondaro-Mallea et al. (2021), the trend flips at intermediate masses and starts to over-predict the abundance of the most massive halos. Similar results were obtained by comparing our calibration directly with the simulations used and made available by Ondaro-Mallea et al. (2021), which reassures us of the robustness of our calibration.

## 5.4. Impact on cluster counts analysis

In this subsection, we forecast the impact of the uncertainties on the HMF calibration on the cosmological constraints to be derived from the Euclid cluster counts.

### 5.4.1. Impact of the halo finder

In Table 5, we summarise the impact of the different halo finders on the inference of the marginalised cosmological parameters $\Omega_m$ and $\sigma_8$. To do so, we compare the results obtained when one assumes VELOCIRAPTOR, SUBFIND or AHF calibrations to create the synthetic catalogue while the fiducial ROCKSTAR calibration is used for the analysis. We perform this forecast following the methodology described in Sect. 2.5. For the likelihood analysis, we assume flat priors on the cosmological parameters and Gaussian priors with amplitudes on the mass-observable parameters of 1, 3, and 5 percent. The likelihood sampling is performed with ZEUS (see Karamanis & Beutler 2021; Karamanis et al. 2021). The impact of the different halo finders' calibrations is quantified using the index of inconsistency (IOI; Lin & Ishak 2017), which is calculated as

$$ \text{IOI} = \frac{\delta^t \Sigma^{-1} \delta}{2}, \qquad (21) $$
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Comparison of our HMF model to the works of Tinker et al. (2008), Despali et al. (2016), and Ondaro-Mallea et al. (2021). We use our ROCKSTAR calibration as the reference and present our SUBFIND calibration for easy comparison with Ondaro-Mallea et al. (2021). We also present the mean of the multiplicity function $\nu f(\nu)$ measured from the ten PICCOLO C0 runs. The grey regions depict the 1 and 5 percent regions. We adopted a symmetric-log scale on the y-axis, where the region between the dashed lines is presented in a linear scale.</figcaption></figure>

where $\delta$ is the two-dimensional difference vector between the best-fit and the assumed cosmological parameters $\{\Omega_m, \sigma_8\} = \{0.30711, 0.8288\}$. Additionally, $\Sigma$ is the covariance matrix between these parameters, which we assume to be Gaussian distributed. In all cases, the tension in the ($\Omega_m, \sigma_8$) plane increases monotonically as the priors on the richness-mass relation tightens. For both VELOCIRaptor and SUBFIND the tension goes from $\lesssim 1 \sigma$ for 5 and 3 percent priors to $\lesssim 2 \sigma$ when the prior tightens to 1 percent. The tension for the AHF case is $\lesssim 1 \sigma$ for all priors considered, which is a result of the similarity observed in Figs. 10 and B.1 between AHF and ROCKSTAR calibrations.

Therefore, we conclude that differences in the HMF calibration associated with different choices of the halo finder propagate into systematic effects in the measurements of cosmological parameters that are comparable to the formal uncertainties on such parameters. For instance, if the cluster richness-mass relation from *Euclid* data could be calibrated at $< 3$ percent precision, then a crucial factor in deriving cosmological constraints from cluster number counts would become the way in which a halo is defined and identified in simulations in the process of the HMF calibration. Lastly, if one increases the error budget for the HMF calibration until it comprises the different halo-finder prescriptions studied here, the IOI presented in Table 5 also provides the expected reduction factor in the FOM of the cosmological constraining power of the *Euclid* cluster counts.

### 5.4.2. Impact of the calibration error

In Table 5, we also summarise the impact of the systematic and statistical errors of our main (ROCKSTAR) calibration and one of the auxiliary calibrations (VELOCIRaptor) on the marginalised uncertainties in the cosmological parameters $\Omega_m$ and $\sigma_8$. We only consider one calibration of each group of halo finders, as this test is dominated by the number of simulations used in the HMF calibration. As VELOCIRaptor, SUBFIND, and AHF all use the same number of simulations — equal to half of the set used for the ROCKSTAR calibration — they present very similar results. We compare the FOM change in the ($\Omega_m, \sigma_8$) plane obtained by fixing the halo mass function parameters to the calibrated values with the ones obtained by marginalizing over such parameters using a multi-variate Gaussian with covariance given by the fit uncertainties. We consider again 1, 3,

and 5 percent priors on the richness-mass scaling relations. For ROCKSTAR the statistical uncertainty only marginally affects the cosmological inference. For VELOCIRaptor, we observe that the only significant impact is seen for the 1 percent prior, where the FOM is over-estimated by ~ 10 percent when the HMF parameters are left fixed. Therefore, from this test, we conclude that the residual uncertainties in the HMF parameters have a negligible impact on the corresponding cosmological constraints.

## 6. Conclusions

In this paper, we carried out a detailed analysis to assess the numerical robustness of the halo mass function (HMF) predicted by N-body simulations, and to quantify and model its deviation from universality. The variety of tests that we carried out include changing the prescription for generating initial conditions, the effect of resolution and of round-off errors in particle positions in the initial conditions, the N-body integrator, the definition of a halo, and the halo finder. While our reference analysis were carried out assuming a vanilla $\Lambda$CDM cosmology, we also simulated the effect on the HMF of including massive neutrinos. Furthermore, in order to trace the origin of departures from universality, we also ran simulations with a purely scale-free power spectrum, both assuming Einstein-de Sitter and $\Lambda$CDM expansion histories. Finally, with the resulting high-resolution calibration of the HMF, we assess the impact of systematic effects in the cosmological parameters inference from an idealized *Euclid* cluster number counts experiment. Our main conclusions can be summarised as follows.

- The different gravity solvers considered in this paper agree better than 1 percent on the HMF of cluster-sized halos for particle masses smaller than $10^{10} h^{-1} M_{\odot}$. Interestingly, the N-body integrator of RAMSES, the only code based on adaptive mesh refinement among those considered here, seems to systematically predict a lower number of halos with $M \lesssim 10^{14} h^{-1} M_{\odot}$ an effect that is more apparent at lower mass resolutions;

- Our adopted setup for the PICCOLO set, which includes simulations for nine different cosmological models, provides a percent-level convergence on the HMF model when compared to higher resolution simulations.
<table><caption>**Table 5.** Summary statistics for the forecast of the impact of different halo finders and calibration errors on *Euclid* cluster counts cosmological constraints on $\Omega_m$ and $\sigma_8$. The IOI quantifies the tension in the posteriors if one uses the ROCKSTAR calibration to create the synthetic data while either the VELOCIRAPTOR, SUBFIND, or AHF calibration is used for the analysis. The relative difference of the FOM assesses the attenuation of the constraining power of cluster counts if one marginalizes over the HMF parameters assuming the calibration chain as a prior. The latter statistics are presented only for the ROCKSTAR and VELOCIRAPTOR calibrations as VELOCIRAPTOR, SUBFIND, and AHF use half the simulations used for the ROCKSTAR calibration and present very similar results. Errors for both statistics were estimated using bootstrap resampling.</caption><thead><tr><th>Summary statistics</th><th>richness-mass relation priors</th><th>Analysis</th><th>Synthetic catalogue</th><th>Value</th></tr></thead><tbody><tr><td rowspan="12">IOI</td><td>1%</td><td rowspan="3">VELOCIRAPTOR<br>(Fixed)</td><td rowspan="3">ROCKSTAR</td><td>1.66 ± 0.01</td></tr><tr><td>3%</td><td>0.77 ± 0.01</td></tr><tr><td>5%</td><td>0.65 ± 0.01</td></tr><tr><td>1%</td><td rowspan="3">SUBFIND<br>(Fixed)</td><td rowspan="3">ROCKSTAR</td><td>1.70 ± 0.02</td></tr><tr><td>3%</td><td>0.84 ± 0.01</td></tr><tr><td>5%</td><td>0.61 ± 0.01</td></tr><tr><td>1%</td><td rowspan="3">AHF<br>(Fixed)</td><td rowspan="3">ROCKSTAR</td><td>0.90 ± 0.02</td></tr><tr><td>3%</td><td>0.61 ± 0.01</td></tr><tr><td>5%</td><td>0.47 ± 0.00</td></tr><tr><td rowspan="12">ΔFOM/FOM</td><td>1%</td><td rowspan="3">ROCKSTAR<br>(marginalised)</td><td rowspan="3">ROCKSTAR</td><td>0.04 ± 0.05</td></tr><tr><td>3%</td><td>0.06 ± 0.04</td></tr><tr><td>5%</td><td>-0.01 ± 0.02</td></tr><tr><td>1%</td><td rowspan="3">VELOCIRAPTOR<br>(marginalised)</td><td rowspan="3">VELOCIRAPTOR</td><td>-0.09 ± 0.05</td></tr><tr><td>3%</td><td>0.00 ± 0.03</td></tr><tr><td>5%</td><td>-0.02 ± 0.03</td></tr></tbody></table>

- Numerical artifacts, such as round-off errors, add non-Poisson fluctuations to the mass-binned distribution of halos. Choosing the mass binning accordingly mitigates the problem.

- The differences in the abundance of halos coming from different halo finders largely dominate all the other numerical systematic errors considered here. The final impact on cosmological constraints depends on how well the mass-observable relation is kept under control, which highlights the need of a better understanding of how halos identified in simulations are related to clusters identified in an optical/NIR photometric survey, such as those provided by *Euclid*.

- The HMF non-universality of virial halos depends on both the background evolution and the power spectrum shape, confirming the results of Ondaro-Mallea et al. (2021).

- Our HMF model was calibrated against nine ΛCDM cosmologies evenly covering the 95 percent confidence level constraints on cosmological parameters from DES+SPT cluster counts (Costanzi et al. 2021) and four different halo finders (ROCKSTAR, AHF, SUBFIND, VELOCIRAPTOR). Our HMF calibration reproduces the abundance of virial halos more massive than $3 \times 10^{13} h^{-1} M_{\odot}$ with a precision of better than 1 percent for the range of cosmological parameters studied here. However, our calibration is expected to retain its accuracy beyond this range, as the HMF modelling was validated on a more extreme set of simulations (AETI-OLOGY).

- Using two external sets of simulations that include massive neutrinos, we validated the model presented in Castorina et al. (2014). Jointly with our calibration, the model by Castorina et al. (2014) reproduces the abundance of virial halos more massive than $4 \times 10^{13} h^{-1} M_{\odot}$ for a total neutrino mass in the range [0.0 – 0.6] eV with a precision of better than 1 percent.

- The statistical uncertainties on the HMF calibration presented in our analysis are significantly smaller than the expected accuracy for the mass-observational relation of *Euclid*. However, the difference between the HMF of the halo finders studied here is comparable to the expected accuracy

for the mass-observational relation of *Euclid* and, as such, could lead to a biased inference of cosmological parameters.

One of the results of our analysis is that the main source of uncertainty in the calibration of the HMF from N-body simulations is related to the definition of halos and to the finder used to identify them. Reassuringly, the differences we found by applying four different halo finders do not compromise the accuracy of the HMF required by the *Euclid* cluster survey. Still, our analysis does not provide information as to the the correspondence between a halo identified in an N-body simulation and a cluster identified in a photometric survey, and how uncertain knowledge of this correspondence will impact on the derivation of cosmological posteriors remains an open question. Indeed, uncertainties in the relation between richness and mass enclosed within a suitably defined (and cosmology-dependent) radius, projection effects in the selection function of the cluster sample, and miscentering, are all effects that need to be controlled and convolved with the predicted HMF. While all such issues need to be addressed one by one, the analysis presented here demonstrates that the precision in the definition of a fitting function for the HMF predicted by N-body simulations of $\Lambda(\nu)$CDM cosmologies is not a limiting factor for cluster cosmology with *Euclid*. In forthcoming analyses we will verify whether a similar precision can be maintained when including uncertainties related to the astrophysics of baryons and departures from the standard $\Lambda(\nu)$CDM framework.

*Acknowledgements.* It is a pleasure to thank Valerio Marra for constructive comments during the production of this work, Renate Mauland-Hus and Hans Winter for the support with RAMSES simulations, Peter Berhoozi for the support with ROCKSTAR, Oliver Hahn for the support with monofonIC and MUSIC, and Douglas Potter for the support with PKDGRAV-3 set up. Lourdes Ondaro-Mallea and Matteo Zennaro for sharing their HMF data. TC, SB, and AS are supported by the INFN INDARK PD51 grant. TC and AS are also supported by the FARE MIUR grant ‘ClustersXEuclid’ R165SBKTMA. AS is also supported by the ERC ‘ClustersXCosmo’ grant agreement 716762. KD acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC-2094 - 390783311 as well as support through the COMPLEX project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679. AR is supported by the PRIN-MIUR 2017 WSCC32 ZOOMING grant. We acknowledge the computing centre
of CINECA and INAF, under the coordination of the “Accordo Quadro (MoU) per lo svolgimento di attività congiunta di ricerca Nuove frontiere in Astrofisica: HPC e Data Exploration di nuova generazione”, for the availability of computing resources and support. We acknowledge the use of the HOTCAT computing infrastructure of the Astronomical Observatory of Trieste – National Institute for Astrophysics (INAF, Italy) (see Bertocco et al. 2020; Taffoni et al. 2020). The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Academy of Finland, the Agenzia Spaziale Italiana, the Belgian Science Policy, the Canadian Euclid Consortium, the French Centre National d'Etudes Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt, the Danish Space Research Institute, the Fundação para a Ciência e a Tecnologia, the Ministerio de Ciencia e Innovación, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Nederlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Romanian Space Agency, the State Secretariat for Education, Research and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (http://www.euclid-ec.org).

## References

Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2020, PRD, 102, 023509

Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409

Amendola, L. & Tsujikawa, S. 2010, Dark energy: theory and observations (Cambridge University Press)

Angulo, R. E. & Hahn, O. 2022, Living Reviews in Computational Astrophysics, 8, 1

Angulo, R. E. & Pontzen, A. 2016, MNRAS, 462, L1

Angulo, R. E. & White, S. D. M. 2010, MNRAS, 401, 1796

Aricò, G., Angulo, R. E., Contreras, S., et al. 2021, MNRAS, 506, 4070

Artis, E., Melin, J.-B., Bartlett, J. G., & Murray, C. 2021, A&A, 649, A47

Bagla, J. S. 2002, J. Astrophys. Astron., 23, 185

Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013a, ApJ, 762, 109

Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., et al. 2013b, ApJ, 763, 18

Bellagamba, F., Roncarelli, M., Maturi, M., & Moscardini, L. 2018, MNRAS, 473, 5221

Bertocco, S., Goz, D., Tornatore, L., et al. 2020, in Astronomical Society of the Pacific Conference Series, Vol. 527, 303

Bhattacharya, S., Heitmann, K., White, M., et al. 2011, ApJ, 732, 122

Bocquet, S., Dietrich, J. P., Schrabback, T., et al. 2019, ApJ, 878, 55

Bocquet, S., Heitmann, K., Habib, S., et al. 2020, ApJ, 901, 5

Bocquet, S., Saro, A., Dolag, K., & Mohr, J. J. 2016, MNRAS, 456, 2361

Bocquet, S., Saro, A., Mohr, J. J., et al. 2015, ApJ, 799, 214

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

Bryan, G. L. & Norman, M. L. 1998, The Astrophysical Journal, 495, 80

Carbone, C., Petkova, M., & Dolag, K. 2016, JCAP, 07, 034

Castignani, G. & Benoist, C. 2016, A&A, 595, A111

Castorina, E., Sefusatti, E., Sheth, R. K., Villaescusa-Navarro, F., & Viel, M. 2014, JCAP, 2014, 049

Castro, T., Borgani, S., Dolag, K., et al. 2020, MNRAS, 500, 2316

Cooray, A. & Sheth, R. K. 2002, Phys. Rept., 372, 1

Corasaniti, P. S. & Achitou, I. 2011a, PRD, 84, 023009

Corasaniti, P. S. & Achitou, I. 2011b, PRL, 106, 241302

Costanzi, M., Saro, A., Bocquet, S., et al. 2021, PRD, 103, 043522

Costanzi, M., Villaescusa-Navarro, F., Viel, M., et al. 2013, JCAP, 12, 012

Courtin, J., Rasera, Y., Alimi, J. M., et al. 2011, MNRAS, 410, 1911

Crocce, M., Fosalba, P., Castander, F. J., & Gaztanaga, E. 2010, MNRAS, 403, 1353

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

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

Dakin, J., Hannestad, S., & Tram, T. 2022, MNRAS, 513, 991

Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371

Despali, G., Giocoli, C., Angulo, R. E., et al. 2016, MNRAS, 456, 2486

Diemand, J., Kuhlen, M., & Madau, P. 2006, ApJ, 649, 1

Diemer, B. 2020, ApJ, 903, 87

Diemer, B. & Joyce, M. 2019, ApJ, 871, 168

Diemer, B. & Kravtsov, A. V. 2015, ApJ, 799, 108

Eke, V. R., Cole, S., & Frenk, C. S. 1996, MNRAS, 282, 263

Elahi, P. J., Cañas, R., Poulton, R. J. J., et al. 2019, Publ. Astron. Soc. Austral., 36, e021

Elahi, P. J., Knebe, A., Pearce, F. R., et al. 2016, MNRAS, 458, 1096

Euclid Collaboration: Adam, R., Vannier, M., Maurogordato, S., et al. 2019, A&A, 627, A23

Euclid Collaboration: Adamek, J. et al. 2022 [arXiv:2211.12457]

Euclid Collaboration: Fumagalli, A., Saro, A., Borgani, S., et al. 2021, A&A, 652, A21

Euclid Collaboration: Scaramella, R., Amiaux, J., Mellier, Y., et al. 2022, A&A, 662, A112

Garcia, R. & Rozo, E. 2019, MNRAS, 489, 4170

Garrison, L. H., Abel, T., & Eisenstein, D. J. 2021a, MNRAS, 509, 2281

Garrison, L. H., Eisenstein, D. J., Ferrer, D., Metchnik, M. V., & Pinto, P. A. 2016, MNRAS, 461, 4125

Garrison, L. H., Joyce, M., & Eisenstein, D. J. 2021b, MNRAS, 504, 3550

Genel, S., Bryan, G. L., Springel, V., et al. 2019, ApJ, 871, 21

Gill, S. P. D., Knebe, A., & Gibson, B. K. 2004, MNRAS, 351, 399

Gonzalez, A. 2014, in Building the Euclid Cluster Survey - Scientific Program, 7

Hahn, O. & Abel, T. 2011, MNRAS, 415, 2101

Hasselfield, M., Hilton, M., Marriage, T. A., et al. 2013, JCAP, 07, 008

Holder, G., Haiman, Z., & Mohr, J. 2001, ApJ Lett., 560, L111

Hu, W. & Kravtsov, A. V. 2003, ApJ, 584, 702

Jenkins, A., Frenk, C. S., White, S. D. M., et al. 2001, MNRAS, 321, 372

Joyce, M., Garrison, L., & Eisenstein, D. 2021, MNRAS, 501, 5051

Joyce, M., Marcos, B., Gabrielli, A., Baertschiger, T., & Labini, F. S. 2005, PRL, 95, 011304

Karamanis, M. & Beutler, F. 2021, Stat. Comput., 31, 61

Karamanis, M., Beutler, F., & Peacock, J. A. 2021, MNRAS, 508, 3589

Kitayama, T. & Suto, Y. 1996, ApJ, 469, 480

Knebe, A., Knollmann, S. R., Muldrew, S. I., et al. 2011, MNRAS, 415, 2293

Knollmann, S. R. & Knebe, A. 2009, ApJ Suppl., 182, 608

Kravtsov, A. V. & Borgani, S. 2012, ARA&A, 50, 353

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv:1110.3193

Leroy, M., Garrison, L., Eisenstein, D., Joyce, M., & Maleubre, S. 2021, MNRAS, 501, 5064

Lesci, G. F., Marulli, F., Moscardini, L., et al. 2022, A&A, 659, A88

Lin, W. & Ishak, M. 2017, PRD, 96, 023532

Maggiore, M. & Riotto, A. 2010a, ApJ, 711, 907

Maggiore, M. & Riotto, A. 2010b, ApJ, 717, 515

Maggiore, M. & Riotto, A. 2010c, ApJ, 717, 526

Maleubre, S., Eisenstein, D., Garrison, L. H., & Joyce, M. 2022, MNRAS, 512, 1829

Mantz, A. B., von der Linden, A., Allen, S. W., et al. 2015, MNRAS, 446, 2205

Marcos, B. 2008, PRD, 78, 043536

McClintock, T., Rozo, E., Becker, M. R., et al. 2019, ApJ, 872, 53

Michaux, M., Hahn, O., Rampf, C., & Angulo, R. E. 2020, MNRAS, 500, 663

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Nishimichi, T., Takada, M., Takahashi, R., et al. 2019, ApJ, 884, 29

Ondaro-Mallea, L., Angulo, R. E., Zennaro, M., Contreras, S., & Aricò, G. 2021, MNRAS, 509, 6077

Peebles, P. J. E. 2020, The large-scale structure of the universe, Vol. 98 (Princeton university press)

Planck Collaboration VI. 2020, A&A, 641, A6, [Erratum: A&A 652, C4 (2021)]

Planck Collaboration XVI. 2014, A&A, 571, A16

Planck Collaboration XX. 2014, A&A, 571, A20

Planck Collaboration XXIV. 2016, A&A, 594, A24

Potter, D., Stadel, J., & Teyssier, R. 2017, Computational Astrophysics and Cosmology, 4, 2

Power, C., Robotham, A. S. G., Obreschkow, D., Hobbs, A., & Lewis, G. F. 2016, MNRAS, 462, 474

Press, W. H. & Schechter, P. 1974, ApJ, 187, 425

Reed, D., Gardner, J., Quinn, T. R., et al. 2003, MNRAS, 346, 565

Reed, D. S., Smith, R. E., Potter, D., et al. 2013, MNRAS, 431, 1866

Rozo, E., Wechsler, R. H., Rykoff, E. S., et al. 2010, ApJ, 708, 645

Salvati, L., Douspis, M., & Aghanim, N. 2020, A&A, 643, A20

Saro, A., Bocquet, S., Rozo, E., et al. 2015, MNRAS, 454, 2305

Sartoris, B., Biviano, A., Fedeli, C., et al. 2016, MNRAS, 459, 1764

Schneider, A. & Teyssier, R. 2015, JCAP, 12, 049

Sheth, R. K., Mo, H. J., & Tormen, G. 2001, MNRAS, 323, 1

Sheth, R. K. & Tormen, G. 1999, MNRAS, 308, 119

Springel, V. 2005, MNRAS, 364, 1105

Springel, V., Pakmor, R., Zier, O., & Reinecke, M. 2021, MNRAS, 506, 2871

Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726

Taffoni, G., Becciani, U., Garilli, B., et al. 2020, in Astronomical Society of the Pacific Conference Series, Vol. 527, 307

Teyssier, R. 2002, A&A, 385, 337

Tinker, J. L., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

Velliscig, M., van Daalen, M. P., Schaye, J., et al. 2019, MNRAS, 442, 2641

Warren, M. S., Abazajian, K., Holz, D. E., & Teodoro, L. 2006, ApJ, 646, 881

Watson, W. A., Iliev, I. T., D'Aloisio, A., et al. 2013, MNRAS, 433, 1230

Xu, G.-H. 1995, ApJ Suppl., 98, 355

Zeldovich, Y. B. 1970, A&A, 5, 84
<footnote n="1">¹ INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, I-34143 Trieste, Italy</footnote>

<footnote n="2">² INFN, Sezione di Trieste, Via Valerio 2, I-34127 Trieste TS, Italy</footnote>

<footnote n="3">³ IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2, I-34151 Trieste, Italy</footnote>

<footnote n="4">⁴ Dipartimento di Fisica - Sezione di Astronomia, Università di Trieste, Via Tiepolo 11, I-34131 Trieste, Italy</footnote>

<footnote n="5">⁵ Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal, 4, 20018, Donostia-San Sebastián, Guipuzkoa, Spain</footnote>

<footnote n="6">⁶ IKERBASQUE, Basque Foundation for Science, 48013, Bilbao, Spain</footnote>

<footnote n="7">⁷ University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, 81679 Munich, Germany</footnote>

<footnote n="8">⁸ INAF-IASF Milano, Via Alfonso Corti 12, I-20133 Milano, Italy</footnote>

<footnote n="9">⁹ Department of Physics and Astronomy, University of Aarhus, Ny Munkegade 120, DK-8000 Aarhus C, Denmark</footnote>

<footnote n="10">¹⁰ Universitäts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstrasse 1, 81679 München, Germany</footnote>

<footnote n="11">¹¹ Max-Planck-Institut für Astrophysik, Karl-Schwarzschild Str. 1, 85741 Garching, Germany</footnote>

<footnote n="12">¹² Istituto Nazionale di Astrofisica (INAF) - Osservatorio di Astrofisica e Scienza dello Spazio (OAS), Via Gobetti 93/3, I-40127 Bologna, Italy</footnote>

<footnote n="13">¹³ Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Via Irnerio 46, I-40126 Bologna, Italy</footnote>

<footnote n="14">¹⁴ Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, I-40129 Bologna, Italy</footnote>

<footnote n="15">¹⁵ Laboratoire Univers et Théorie, Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, F-92190 Meudon, France</footnote>

<footnote n="16">¹⁶ Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK</footnote>

<footnote n="17">¹⁷ Institut für Theoretische Physik, University of Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany</footnote>

<footnote n="18">¹⁸ INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93/3, I-40129 Bologna, Italy</footnote>

<footnote n="19">¹⁹ INFN-Sezione di Bologna, Viale Berti Pichat 6/2, I-40127 Bologna, Italy</footnote>

<footnote n="20">²⁰ Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1, D-85748 Garching, Germany</footnote>

<footnote n="21">²¹ Dipartimento di Fisica, Università di Genova, Via Dodecaneso 33, I-16146, Genova, Italy</footnote>

<footnote n="22">²² INFN-Sezione di Roma Tre, Via della Vasca Navale 84, I-00146, Roma, Italy</footnote>

<footnote n="23">²³ Department of Physics "E. Pancini", University Federico II, Via Cinthia 6, I-80126, Napoli, Italy</footnote>

<footnote n="24">²⁴ INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, I-80131 Napoli, Italy</footnote>

<footnote n="25">²⁵ Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, I-10125 Torino, Italy</footnote>

<footnote n="26">²⁶ INFN-Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy</footnote>

<footnote n="27">²⁷ INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, I-10025 Pino Torinese (TO), Italy</footnote>

<footnote n="28">²⁸ Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain</footnote>

<footnote n="29">²⁹ Port d'Informació Científica, Campus UAB, C. Albareda s/n, 08193 Bellaterra (Barcelona), Spain</footnote>

<footnote n="30">³⁰ INAF-Osservatorio Astronomico di Roma, Via Frascati 33, I-00078 Monteporzio Catone, Italy</footnote>

<footnote n="31">³¹ INFN section of Naples, Via Cinthia 6, I-80126, Napoli, Italy</footnote>

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

<footnote n="33">³³ INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125, Firenze, Italy</footnote>

<footnote n="34">³⁴ Centre National d'Etudes Spatiales, Toulouse, France</footnote>

<footnote n="35">³⁵ Institut national de physique nucléaire et de physique des particules, 3 rue Michel-Ange, 75794 Paris Cédex 16, France</footnote>

<footnote n="36">³⁶ Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="37">³⁷ ESAC/ESA, Camino Bajo del Castillo, s/n., Urb. Villafranca del Castillo, 28692 Villanueva de la Cañada, Madrid, Spain</footnote>

<footnote n="38">³⁸ European Space Agency/ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma, Italy</footnote>

<footnote n="39">³⁹ Univ Lyon, Univ Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822, F-69622, Villeurbanne, France</footnote>

<footnote n="40">⁴⁰ Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland</footnote>

<footnote n="41">⁴¹ Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, Edificio C8, Campo Grande, PT1749-016 Lisboa, Portugal</footnote>

<footnote n="42">⁴² Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, PT-1749-016 Lisboa, Portugal</footnote>

<footnote n="43">⁴³ Department of Astronomy, University of Geneva, ch. d'Écogia 16, CH-1290 Versoix, Switzerland</footnote>

<footnote n="44">⁴⁴ Université Paris-Saclay, CNRS, Institut d'astrophysique spatiale, 91405, Orsay, France</footnote>

<footnote n="45">⁴⁵ Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="46">⁴⁶ Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK</footnote>

<footnote n="47">⁴⁷ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, Astro-physique, Instrumentation et Modélisation Paris-Saclay, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="48">⁴⁸ Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193 Barcelona, Spain</footnote>

<footnote n="49">⁴⁹ Institut d'Estudis Espacials de Catalunya (IEEC), Carrer Gran Capità 2-4, 08034 Barcelona, Spain</footnote>

<footnote n="50">⁵⁰ INAF-Osservatorio Astronomico di Padova, Via dell'Osservatorio 5, I-35122 Padova, Italy</footnote>

<footnote n="51">⁵¹ Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, N-0315 Oslo, Norway</footnote>

<footnote n="52">⁵² von Hoerner & Sulger GmbH, SchloßPlatz 8, D-68723 Schwetzingen, Germany</footnote>

<footnote n="53">⁵³ Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark</footnote>

<footnote n="54">⁵⁴ Institut d'Astrophysique de Paris, 98bis Boulevard Arago, F-75014, Paris, France</footnote>

<footnote n="55">⁵⁵ Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="56">⁵⁶ Aix-Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France</footnote>

<footnote n="57">⁵⁷ Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK</footnote>

<footnote n="58">⁵⁸ Université de Genève, Département de Physique Théorique and Centre for Astroparticle Physics, 24 quai Ernest-Ansermet, CH-1211 Genève 4, Switzerland</footnote>

<footnote n="59">⁵⁹ Department of Physics and Helsinki Institute of Physics, Gustaf Hällströminkatu 2, 00014 University of Helsinki, Finland</footnote>

<footnote n="60">⁶⁰ NOVA optical infrared instrumentation group at ASTRON, Oude Hoogeveensedijk 4, 7991PD, Dwingeloo, The Netherlands</footnote>

<footnote n="61">⁶¹ Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany</footnote>

<footnote n="62">⁶² Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, via Piero Gobetti 93/2, I-40129 Bologna, Italy</footnote>

<footnote n="63">⁶³ European Space Agency/ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands</footnote>

<footnote n="64">⁶⁴ AIM, CEA, CNRS, Université Paris-Saclay, Université de Paris, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="65">⁶⁵ Space Science Data Center, Italian Space Agency, via del Politecnico snc, 00133 Roma, Italy</footnote>

<footnote n="66">⁶⁶ Institute of Space Science, Bucharest, Ro-077125, Romania</footnote>

<footnote n="67">⁶⁷ Instituto de Astrofísica de Canarias, Calle Vía Láctea s/n, E-38204, San Cristóbal de La Laguna, Tenerife, Spain</footnote>

<footnote n="68">⁶⁸ Departamento de Astrofísica, Universidad de La Laguna, E-38206, La Laguna, Tenerife, Spain</footnote>

<footnote n="69">⁶⁹ Dipartimento di Fisica e Astronomia "G.Galilei", Universitá di Padova, Via Marzolo 8, I-35131 Padova, Italy</footnote>

<footnote n="70">⁷⁰ INFN-Padova, Via Marzolo 8, I-35131 Padova, Italy</footnote>
71 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, 91109, USA

72 Departamento de Física, FCFM, Universidad de Chile, Blanco En-
calada 2008, Santiago, Chile

73 Centro de Investigaciones Energéticas, Medioambientales y Tec-
nológicas (CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain

74 Universidad Politécnica de Cartagena, Departamento de Electrónica y Tecnología de Computadoras, 30202 Cartagena, Spain

75 Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands

76 Infrared Processing and Analysis Center, California Institute of Tech-
nology, Pasadena, CA 91125, USA

77 INAF-Osservatorio Astronomico di Brera, Via Brera 28, I-20122 Mi-
lano, Italy

78 INAF-Istituto di Astrofisica e Planetologia Spaziali, via del Fosso del Cavaliere, 100, I-00100 Roma, Italy

79 Junia, EPA department, F 59000 Lille, France

80 SISSA, International School for Advanced Studies, Via Bonomea 265, I-34136 Trieste TS, Italy

81 Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara, Via Giuseppe Saragat 1, I-44122 Ferrara, Italy

82 Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara, Via Giuseppe Saragat 1, I-44122 Ferrara, Italy

83 Institut de Physique Théorique, CEA, CNRS, Université Paris-Saclay F-91191 Gif-sur-Yvette Cedex, France

84 Institut de Recherche en Astrophysique et Planétologie (IRAP), Uni-
versité de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, F-
31400 Toulouse, France

85 NASA Ames Research Center, Moffett Field, CA 94035, USA

86 INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, I-40129 Bologna, Italy

87 INFN-Bologna, Via Irnerio 46, I-40126 Bologna, Italy

88 Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Lab-
oratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice cedex
4, France

89 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Tapada da Ajuda, PT-1349-018 Lisboa, Portugal

90 Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, D-52056 Aachen, Germany

91 Department of Physics & Astronomy, University of California Irvine, Irvine CA 92697, USA

92 University of Lyon, UCB Lyon 1, CNRS/IN2P3, IUF, IP2I Lyon, France

93 INFN-Sezione di Genova, Via Dodecaneso 33, I-16146, Genova, Italy

94 Université Paris Cité, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France

95 Instituto de Física Teórica UAM-CSIC, Campus de Cantoblanco, E-
28049 Madrid, Spain

96 Department of Physics, P.O. Box 64, 00014 University of Helsinki, Finland

97 Ruhr University Bochum, Faculty of Physics and Astronomy, Astro-
nomical Institute (AIRUB), German Centre for Cosmological Lensing
(GCCL), 44780 Bochum, Germany

98 Department of Physics, Lancaster University, Lancaster, LA1 4YB, UK

99 Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France

100 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK

101 Astrophysics Group, Blackett Laboratory, Imperial College London, London SW7 2AZ, UK

102 Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, Av-
enue des Martyrs, 38000, Grenoble, France

103 Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 2, I-00185 Roma, Italy

104 University of Applied Sciences and Arts of Northwestern Switzer-
land, School of Engineering, 5210 Windisch, Switzerland

105 Aix-Marseille Univ, CNRS, CNES, LAM, Marseille, France

106 Zentrum für Astronomie, Universität Heidelberg, Philosophenweg 12, D-69120 Heidelberg, Germany

107 Department of Mathematics and Physics E. De Giorgi, University of Salento, Via per Arnesano, CP-193, I-73100, Lecce, Italy

108 INAF-Sezione di Lecce, c/o Dipartimento Matematica e Fisica, Via per Arnesano, I-73100, Lecce, Italy

109 INFN, Sezione di Lecce, Via per Arnesano, CP-193, I-73100, Lecce, Italy

110 Observatoire de Sauverny, Ecole Polytechnique Fédérale de Lau-
sanne, CH-1290 Versoix, Switzerland

111 Institute for Computational Science, University of Zurich, Win-
terthurerstrasse 190, 8057 Zurich, Switzerland

112 Higgs Centre for Theoretical Physics, School of Physics and Astron-
omy, The University of Edinburgh, Edinburgh EH9 3FD, UK

113 Université St Joseph; Faculty of Sciences, Beirut, Lebanon

114 Department of Astrophysical Sciences, Peyton Hall, Princeton Uni-
versity, Princeton, NJ 08544, USA

115 Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland

116 Department of Mathematics and Physics, Roma Tre University, Via della Vasca Navale 84, I-00146 Rome, Italy
<figure><img src="image_20.png" /></figure>

**Fig. A.1.** Sensitivity of the cumulative mass function to the accuracy parameters of Open-GADGET. Each line corresponds to the ratio of the cumulative mass function obtained setting one of the accuracy parameters to half the value presented in Table 3 with respect to the fiducial set. For comparison, we add the 68% confidence level in red for the ErrTolForceAcc assuming that halos in the simulations are distributed according to a Poisson distribution.

## Appendix A: Convergence tests

In Fig. A.1, we test the sensitivity and convergence of the cumulative HMF extracted from AHF catalogues to the accuracy parameters of Open-GADGET for one of the TEASE simulations of 512³ particles displaced from an equally spaced grid according to the Zeldovich approximation at $z = 99$. Each line corresponds to the ratio of the cumulative HMF obtained setting one of the accuracy parameters to half the value presented in Table 3. For comparison, we also add the 68 percent confidence level in red for the ErrTolForceAcc assuming that the number of halos in the simulations are distributed according to a Poisson distribution. Fig. A.1 shows that our parameter set provide sub-Poisson differences to higher precision sets, deviating by less than a fraction of a percent in the abundance of halos more massive than $3 \times 10^{13} h^{-1} M_{\odot}$.

In Fig. A.2, we present the comparison between the cumulative mass function extracted with AHF from cosmological simulations run with Open-GADGET and GADGET-4. We tested the convergence of the results for three different initial conditions: 512³ particles displaced from an equally spaced grid according to 3LPT at $z = 24$, the same number of particles using Zeldovich at $z = 99$, and $4 \times 320^3$ particles displaced from a face centred cubic (FCC) grid according to 3LPT at $z = 24$. The box size is $500 h^{-1} \text{ Mpc}$. In all configurations, the agreement between Open-GADGET and GADGET-4 is better than a fraction of percent for all masses considered. Due to the higher degree of planes of symmetry, the FCC configuration shows an even better agreement between the two codes, as in this configuration the usual tree algorithm delivers more accurate forces from the simulation start.

In Fig. A.3, we present a convergence test for the cumulative mass function computed from AHF halo catalogue for different particle masses and spatial resolutions. We consider three mass resolutions, corresponding to the following particle masses: $\{6.35 \times 10^{11}, 7.94 \times 10^{10}, 9.92 \times 10^9\} h^{-1} M_{\odot}$; respectively $\{256^3, 512^3, 1024^3\}$ particles in a box of $500 h^{-1} \text{ Mpc}$. For the gravitational softening we consider: $\epsilon = \{1/20, 1/40, 1/80\}$, where *l* is the mean interparticle distance. For reference, the 1 percent region is shown with the grey-shaded area, while the 68 percent confidence level assuming Poisson distribution for

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

**Fig. A.2.** Comparison between the cumulative mass function extracted from Open-GADGET and GADGET-4 for three different initial conditions: 512³ displaced from a equally spaced grid according to 3LPT at $z = 24$, the same number of particles using Zeldovich at $z = 99$, and $4 \times 320^3$ particles displaced from a FCC grid at $z = 24$.

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

**Fig. A.3.** Convergence test for the cumulative mass function for different particle mass ($\{6.35 \times 10^{11}, 7.94 \times 10^{10}, 9.92 \times 10^9\} h^{-1} M_{\odot}$, respectively, $\{256^3, 512^3, 1024^3\}$ particles in a box of $500 h^{-1} \text{ Mpc}$) and spatial resolution set by the gravitational softening ($\epsilon = \{1/20, 1/40, 1/80\}$, where *l* is the mean interparticle distance). The filled region in grey represents the 1 percent region around unity and the filled region in red represents the 68 percent confidence level assuming Poisson distribution for the abundance of halos in both simulations.

the halo distribution of the simulations with $1024^3$ particles with $\epsilon = 1/40$ and $\epsilon = 1/80$ is depicted in red. We note that a particle mass of few times $10^{10} h^{-1} M_{\odot}$ is enough to have a cumulative mass function that agrees to within 1 percent with the higher resolution simulation for halos more massive than $10^{14} h^{-1} M_{\odot}$. Lighter objects are more strongly affected by both the mass and spatial resolutions. For the simulation with $256^3$ particles, even the formation of very massive objects is suppressed by $\approx 4$ percent.

## Appendix B: Comparison between the different halo finder calibrations

The ratio of the multiplicity function as a function of the halo virial mass for the different calibrations with respect to ROCKSTAR is shown in Fig. B.1 colour coded by the redshift. Fig. B.1 confirms the trends at low redshift presented in Fig. 10, but also shows their dependence with redshift. Despite the very different algorithm, AHF and ROCKSTAR provide a multiplicity func-
<figure><img src="image_23.png" /></figure>

**Fig. B.1.** Ratio of the multiplicity function as a function of the halo virial mass for the different calibrations to ROCKSTAR colour coded according to redshift. From *top* to *bottom*, we present the results for VELOCIRaptor, SUBFIND, and AHF.

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

**Fig. C.1.** Dependence of ROCKSTAR fit quality as a function of redshift. The fit quality (p-value) was estimated from sampling from the likelihood used for the calibration assuming the best-fit parameters. The p-value is the fraction of random catalogues that had a residual higher than the original catalogue. On the right panel, we present the histogram projection of the left panel.

## Appendix C: Assessment of the fit-quality dependence on redshift

In Fig. C.1, we present the ROCKSTAR fit quality dependence with redshift. The fit-quality (p-value) in Fig. C.1 was estimated from sampling 1000 points from the likelihood used for the calibration assuming the best-fit parameters. On the right panel, we present the histogram projection of the left panel. The p-value is the fraction of random catalogues that had a residual higher than the original catalogue fit. While all the PICCOLO snapshots were used to produce the left panel of the Fig. C.1, only the specific snapshots used for the calibration were used to present the histogram on the right panel so as to avoid using strongly correlated snapshots. Figure C.1 shows that our model provides robust results for the redshift range [0, 2] and confirms that our choice of $\sigma_{\text{sys}}$ is reasonable.

tion that agrees within 10 percent up to $M_{\text{vir}} \sim 10^{15} h^{-1} M_{\odot}$ for the redshift range considered.

ROCKSTAR's agreement with VELOCIRaptor is worse than the agreement with AHF, depending more strongly on both mass and redshifts. At low redshift, VELOCIRaptor tends to suppress the number of halos with mass $M_{\text{vir}} \lesssim 10^{14} h^{-1} M_{\odot}$ by few percents.

ROCKSTAR's worst agreement is with respect to SUBFIND, with SUBFIND suppressing the number of objects with mass $M_{\text{vir}} \lesssim 10^{15} h^{-1} M_{\odot}$ at all redshifts.