Euclid preparation. XXIV. Calibration of the halo mass function in Λ ( ν ) CDM cosmologies

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 e ﬀ ects in the simulation. First, we test the convergence of HMF predictions from di ﬀ erent N -body codes, by using initial conditions generated with di ﬀ erent orders of Lagrangian Perturbation theory, and adopting di ﬀ erent simulation box sizes and mass resolution. Then, we quantify the e ﬀ ect of using di ﬀ erent halo ﬁnder algorithms, and how the resulting di ﬀ erences 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 di ﬀ erent spectral indexes, assuming both Einstein–de Sitter and standard Λ CDM expansion histories. Based on these results, we construct a ﬁtting function for the HMF that we demonstrate to be sub-percent accurate in reproducing results from 9 di ﬀ erent variants of the Λ 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 e ﬀ ect due to the halo ﬁnder, that could lead to biased cosmological inference.


Introduction
Structure formation in the Universe follows a hierarchical process, with larger structures forming from the collapse of smaller Article number, page 1 of 25 arXiv:2208.02174v2[astro-ph.CO] 16 Mar 2023 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 excursionset 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 , 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 Nbody 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.

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).

The halo mass function
The comoving number density of halos with mass in the range [M, M + dM] is given by: Equation ( 1) is known as the differential halo mass function, where ρ m is the comoving cosmic mean matter density, ν is the peak height, and the function ν f (ν) is known as the multiplicity function.The peak height ν is a measure of how rare a halo is and is defined as ν = δ c /σ(M, z), where δ c is the critical density for spherical halo collapse extrapolated to z = 0 (Peebles 2020) and σ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: where R(M) = 3 M / (4 π ρ m )1/3 is the Lagrangian radius of a sphere containing the mass M, and W(k R) 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 ν, 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 δ c on cosmology.For instance, although several works have used the mass variance σ 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 ν(M, z) = δ c /σ(M, z), thus including the cosmological dependence of the linear collapse barrier δ 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 ∆ vir (z) times the background density, where ∆ 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), (3) The dependence of the fitting parameters {a, p, q} on cosmology is presented in Sect. 4. As for δ 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 ν, this condition is satisfied by: where Γ 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.

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 sim-ulAtion 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 Mpc cosmological boxes run with different setups and different codes.We used the codes: Open-GADGET, GADGET-4, 1 PKDGRAV-3, 2 CONCEPT, 3 and RAMSES, 4 which respectively deploy the following gravity solver algorithms: tree-particle mesh (tree-PM; Xu 1995;Bagla 2002;Springel 2005), fast multipole methodparticle mesh (FMM-PM; Springel et al. 2021), (FMM;Potter et al. 2017), particle particle -particle mesh (P 3 M; Dakin et al. 2022), and adaptive mesh refinement (AMR; Teyssier 2002).Initial conditions have been generated using either the Zeldovich approximation (Zeldovich 1970) at z = 99 with MUSIC (Hahn  & Abel 2011) 5 or third-order Lagrangian perturbation theory (3LPT) at z = 24 with monofonIC (Michaux et al. 2020). 6For 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 3 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 ΛCDM linear power spectrum on a background that is either Einstein-de Sitter -that is (Ω m , Ω Λ ) = (1, 0); EdS hereafter -or ΛCDM.The selfsimilarity 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), andMaleubre 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 × 1280 3 particles.These simulations were carried out with Open-GADGET.Initial conditions were generated using mono-fonIC 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-darkmatter 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.

Halo finders
For the purpose of testing the HMF calibration against the choice of the halo finder, we used AHF,7 SUBFIND, 8 VELOCIraptor,9 DENHF, and ROCKSTAR10 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(Springel et al. , 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 CONSISTENT11 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.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 δ 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 L(N i |θ θ θ, z) for the number of halos where N i (θ θ θ, z) is the theoretical expectation of halos computed by integrating Eq. ( 1) and multiplying it by the volume of the simulated cosmological box.Lastly, θ θ θ 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: where σ = N sim i + σ 2 sys 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 σ 2 sys .We note that our likelihood presents a discontinuity at N sim i = 25.The impact of the discontinuity is two-fold.Firstly, the discontinuity amplitude depends on the chosen value for σ sys .As we discuss here below, we assume σ sys 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 sim i guarantees that the transition is again smooth as the Poisson contribution largely dominates the total error budget while negative values for N i (θ θ θ, z) are very rare ( 5 σ) fluctuations of a Gaussian distribution.We also checked a more complex composite likelihood function with a two-sided σ following Watson et al. (2013): 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.

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] (Laureijs 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 λ is assigned to the halos according to the richness-mass relation λ|M vir , z (see e.g.Saro et al. 2015): where E(z) is the ratio of the Hubble parameter at redshift z and 0. We assume a richness range λ = [20, 2000] and a log-normal scatter given by: We use the following fiducial values for the parameters of Eqs. ( 8) and ( 9) A λ = 37.8, B λ = 1.16,C λ = 0.91, D λ = 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).

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).

Parameter and configuration setting
The chosen values for the internal key simulation parameters used within Open-GADGET are presented in Table 3. Param-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.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. 12e 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×10 13 h −1 M .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-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 × 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.

Resolution convergence
Previous works (see e.g.Joyce et al. 2005;Marcos 2008;Garrison 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 conditions generated with higher order Lagrangian perturbation theory 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, corresponding to particle masses: {6.50×10 11 , 8.13×10 10 , 1.02×10 10 , 5.08× 10 9 , 1.27 × 10 9 } h −1 M ; these values for the particle masses correspond to {4 × 160 3 , 4 × 320 3 , 4 × 640 3 , 1280 3 , 4 × 1280 3 } particles in a box of 500 h −1 Mpc, respectively.The rationale behind such choices for the number of particles is that 4 × 160 3 differs by ∼ 2 percent from 256 3 , thus one can directly compare with other commonly used particle numbers at similar computational cost.The gravitational softening is set to one-fortieth of the mean inter-particle distance in all simulations.Initial conditions are created by monofonIC using 3LPT at z = 24 and FCC grid for all simulations but the simulation with 1280 3 particles, which uses a standard equally spaced grid.The reference in all redshift panels is the simulation run at the highest resolution.For easier inspection, we add the moving average over five bins for each case.We note that the lower resolution simulation 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 × 10 13 h −1 M .For the other redshifts, 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 assuming 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 estimation of the true scatter.
From the next simulation in resolution order, we see a subpercent agreement in the number of objects more massive than 5 × 10 13 h −1 M .At the same computational cost, the simulation with 4 × 320 3 (Fig. 1) particles presents a slightly better convergence than 512 3 (Fig. A.3) at z = 0. Lastly, comparing the two most costly simulations at z = 1.98, we observe that the simulation 4 × 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 particles; this illustrates one of the advantages of using FCC grids instead of the standard ones for creating initial conditions.

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 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 radius 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} × 10 15 h −1 M 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, respectively.For the sake of clarity, a mass threshold of (10 12 ) 10 13 h −1 M was imposed to select the (sub)halos presented in Fig. 2. We observe that all N-body codes produce similar distributions of the most massive objects, however, due to slight differences 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 distribution of substructures identified by AHF.While in this paper we are not interested in any detailed analysis of substructures, we trace them here for the sole purpose of comparing different N-body solvers in detail.Open-GADGET, GADGET-4, PKDGRAV-3, and CONCEPT produce similar numbers of substructures in large objects, but a very heterogeneous spatial distribution of them.On the other hand, we observe that RAMSES produces a smoother mass distribution than the other codes, significantly reducing the number of detected substructures.Similar 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 parameters are usually set as a function of the number of particles; 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 × 10 11 , 7.94 × 10 10 , 9.92 × 10 9 } h −1 M ; respectively {256 3 , 512 3 , 1024 3 } particles in a box of 500 h −1 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 × 10 13 h −1 M , the difference between RAMSES and Open-GADGET is quite sensitive to resolution.At the lowest resolution considered in Fig. 3  than ≈ 3 × 10 14 h −1 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 13 h −1 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.

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).

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 percentlevel accuracy on the cumulative mass function for objects more massive than 5 × 10 13 h −1 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 13 h −1 M , the 3LPT initial conditions set at z = 24 slightly boosts the formation of structures.The reason for this is 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-  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 hydrodynamic 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 singleprecision 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 13 and 10 15 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 13 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 σ 0.4 m p M vir ≡ 0.4 m p 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 Article number, page 10 of 25

PDF
M vir = 3.0 × 10 13 h 1 M M vir = 6.9 × 10 13 h 1 M M vir = 1.8 × 10 14 h 1 M M vir = 6.6 × 10 14 h 1 M ( = 0.0, = 0.4) 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 doubleprecision 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 .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: where N min is the number of particles of the smallest object of interest.

Impact of the simulated volume
The effect of the simulation volume on the HMF is two-fold: firstly, the computation of the mass variance σ 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 σ P (k) given by: where N(k) = 2π (δ k /k) 2 is the number of independent modes inside the bin with width δ 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×10 14 h −1 M .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 × 10 14 h −1 M .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 × 10 14 h −1 M and 2 × 10 15 h −1 M , 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' 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, 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.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.

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 SUB-FIND 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 × 10 15 h −1 M while the smallest halo marked in each panel is 10 13 h −1 M .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).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 .
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.

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.

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.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 percent.That is due to the fact that the gravitationally bound mass within a fixed halo radius is by construction smaller than the total 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 lensing signal.In the remainder of this paper, we refer only to the total mass definition.(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 percent for the high-end tail of the mass function (see e.g.Diemer 2020).
In the following sections, and for the specific purpose of tracing 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 otherwise, 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 particles.
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 σ sys = 0.01 N sim i .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- 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}.

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, The first interesting result is that, even for self-similar cosmologies, the HMF is not universal against changing the spectral index.While the parameters {p, q} of Eq. (3) seem to have a linear dependence on n s , the parameter a exhibits a non-monotonic dependence with a local minimum around n s = −1.75.The overall 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.

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 agreement with the Planck 2018 (P18) results for the cosmological parameters (Planck Collaboration VI. 2020) or EdS.For the spe- 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.

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: where: 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 Ω m (z).For the cosmological models considered here this is well approximated by (see e.g.Amendola & Tsujikawa 2010): with γ = 0.55.Therefore, Eqs. ( 12) to ( 16) produce exactly the same results for the cosmological models studied here once one substitutes: 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 Ω 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 σ(R)/d ln R. For the power-law cosmologies this reduces to:

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).

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 Open-GADGET 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 4. Parameters of the multiplicity function f (ν, z) presented in Eq. 3 for the PICCOLO set of simulations.
a 1 a 2 a z p 1 p 2 q 1 q 2 q z ROCKSTAR 0.7962 0.1449 −0.0658 −0.5612 −0.4743 0.3688 −0.2804 0.0251 AHF 0.7937 0.1119 −0.0693 −0.5689 −0.4522 0.3652 −0.2628 0.0376 SUBFIND 0.7953 0.1667 −0.0642 −0.6265 −0.4907 0.3215 −0.2993 0.0330 VELOCIraptor (adaptive 6DFOF) 0.7987 0.1227 −0.0523 −0.5912 −0.4088 0.3634 −0.2732 0.0715 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 δ ln M ≈ 0.15.We assumed the likelihood presented in Sect.2.4 with σ sys = 0.005 N sim i .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 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.

Cosmologies with massive neutrinos
The response of the HMF to neutrinos R is defined as where M ν 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 R predicted by the model and from simulations.The rationale for using 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 DEM-NUni (see Carbone et al. 2016, for details of the simulation) 14and 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). 15The reference simulation considers massless 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 3 dark matter particles, while the latter includes the same number of dark matter particles and as many 2048 3 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 DEMNUni 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 × 10 13 h −1 M 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 DEMNUni, 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 × 10 13 h −1 M is again consistent with 1 percent scatter.

Comparison to other works
In Fig. 18, we compare our HMF model to the works of Tinker et al. (2008), Despali et al. (2016), andOndaro-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 ν f (ν) 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 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 ; it crosses the 5 percent threshold at 2 × 10 15 h −1 M , 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×10 14 h −1 M , and over-predicts the number of halos more massive than 10 15 h −1 M 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 × 10 14 h −1 M , 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.

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.

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 Ω m and σ 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  2021).We also present the mean of the multiplicity function ν f (ν) 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.
where δ is the two-dimensional difference vector between the best-fit and the assumed cosmological parameters {Ω m , σ 8 } = {0.30711,0.8288}.Additionally, Σ is the covariance matrix between these parameters, which we assume to be Gaussian distributed.In all cases, the tension in the (Ω m , σ 8 ) plane increases monotonically as the priors on the richness-mass relation tightens.For both VELOCIraptor and SUBFIND the tension goes from 1 σ for 5 and 3 percent priors to 2 σ when the prior tightens to 1 percent.The tension for the AHF case is 1 σ 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 .

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 Ω m and σ 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 (Ω m , σ 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.

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 Λ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 Λ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 .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 10 14 h −1 M 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.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 massobservable 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 × 10 13 h −1 M 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 × 10 13 h −1 M 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 miscentring, 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 Λ(ν)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 Λ(ν)CDM framework.

Fig. 1 .
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 11 , 8.13 × 10 10 , 1.02 × 10 10 , 5.08 × 10 9 , 1.27 × 10 9 } h −1 M ; respectively {4 × 160 3 , 4 × 320 3 , 4 × 640 3 , 1280 3 , 4 × 1280 3 } particles in a box of 500 h −1 Mpc.The rationale for the chosen number of particles is that 4 × 160 3 differs by ∼ 2 percent from 256 3 , 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.

mFig. 2 .Fig. 3 .
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 was imposed to select the (sub)halos presented.

Fig. 4 .
Fig. 4. Ratio of the cumulative mass function measured from simulations started with same seed, but different LPT order (3LPT vs. Zeldovich approximation) and starting redshifts (z = 24, z = 49, and z = 99), as indicated by the different coloured lines in the legend.

Fig. 5 .
Fig. 5. Distribution of the relative difference between masses of halos at z = 0 matched between different simulations with initial conditions perturbed by 1 pc.Top: Distribution of the relative mass difference.Bottom: Same distribution scaled by the number of particles corresponding to the object mass.

Fig. 8 .
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.

Fig. 9 .Fig. 10 .
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 × 10 15 h −1 M and 10 13 h −1 M , respectively.The cyan circles denote the virial radius of halos identified in this region.

Fig. 11 .
Fig. 11.Difference in the cumulative mass function when all particles inside the spherical overdensity are considered versus just the bound particles.

Fig. 12 .
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) ≈ k ns and Ω m = 1.

Fig. 13 .
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.

Fig. 14 .
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.

Fig. 15 .
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.

Fig. 16 .
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).

Fig. 17 .
Fig. 17.Accuracy of the model presented inCastorina 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 ofCastorina 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).

Fig. 18 .
Fig. 18.Comparison of our HMF model to the works ofTinker et al. (2008),Despali et al. (2016), andOndaro-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 ν f (ν) 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.

Table 1 .
The suite of simulations.The corresponding cosmology is given in Table2.

Table 2 .
The cosmological parameters of the simulations presented in
, RAMSES agrees with the other codes for halos more massive Castro et al.: Euclid preparation.XXIV.Calibration of the halo mass function in Λ(ν)CDM cosmologies Ondaro-Mallea et al. (2021)dOndaro-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 byOndaro-Mallea et al. (2021), which reassures us of the robustness of our calibration.

Table 5 .
Summary statistics for the forecast of the impact of different halo finders and calibration errors on Euclid cluster counts cosmological constraints on Ω m and σ 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.Numerical artifacts, such as round-off errors, add non-Poisson fluctuations to the mass-binned distribution of halos.