Non-universality of the mass function: dependence on the growth rate and power spectrum shape
Lurdes Ondaro-Mallea, Raul E. Angulo, Matteo Zennaro, Sergio Contreras, Giovanni Aricò
MMNRAS , 1–14 (2021) Preprint 19 February 2021 Compiled using MNRAS L A TEX style file v3.0
Non-universality of the mass function: dependence on thegrowth rate and power spectrum shape
Lurdes Ondaro-Mallea, , (cid:63) Raul E. Angulo, , † Matteo Zennaro, Sergio Contreras, and Giovanni Aric`o , . Donostia International Physics Center (DIPC), Manuel Lardizabal Ibilbidea, 4, 20018 Donostia, Gipuzkoa, Spain Universidad Aut´onoma de Madrid (UAM), C/ Francisco Tom´as y Valiente, 7, 28049 Madrid, Spain IKERBASQUE, Basque Foundation for Science, 48013, Bilbao, Spain Universidad de Zaragoza, Pedro Cerbuna 12, 50009 Zaragoza, Spain
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
The abundance of dark matter haloes is one of the key probes of the growth of structure and expansion historyof the Universe. Theoretical predictions for this quantity usually assume that, when expressed in a certain form, itdepends only on the mass variance of the linear density field. However, cosmological simulations have revealed thatthis assumption breaks, leading to 10-20% systematic effects. In this paper we employ a specially-designed suite ofsimulations to further investigate this problem. Specifically, we carry out cosmological N -body simulations wherewe systematically vary growth history at a fixed linear density field, or vary the power spectrum shape at a fixedgrowth history. We show that the halo mass function generically depends on these quantities, thus showing a clearsignal of non-universality. Most of this effect can be traced back to the way in which the same linear fluctuationgrows differently into the nonlinear regime depending on details of its assembly history. With these results, wepropose a parameterization with explicit dependence on the linear growth rate and power spectrum shape. Using anindependent suite of simulations, we show that this fitting function accurately captures the mass function of haloesover cosmologies spanning a vast parameter space, including massive neutrinos and dynamical dark energy. Finally,we employ this tool to improve the accuracy of so-called cosmology-rescaling methods and show they can deliver 2%accurate predictions for the halo mass function over the whole range of currently viable cosmologies. Key words: cosmology: theory – large-scale structure of Universe – methods: statistical – methods: numerical
Collapsed dark matter structures, a.k.a. haloes, offer an im-portant way to constrain fundamental properties of the Uni-verse. The abundance of haloes is sensitive to the growth ofstructure and the statistics of primordial fluctuations, thus itcan be employed to, for instance, constrain the value of cos-mic parameters including dark energy and the sum of neu-trino masses (Weinberg et al. 2013).In the next decades, up to hundreds of thousands of haloeswith mass above ∼ h − M (cid:12) will be detected by upcomingobservational surveys (e.g. eROSITA (Hofmann et al. 2017), EUCLID (Sartoris et al. 2016),
LSST (Ivezi´c et al. 2019),
Simons Observatory (Ade et al. 2019),
CMB-S4 (Abazajianet al. 2019), and
J-PAS (Bonoli et al. 2020)). These futuresurveys will employ various observables over different wave-lengths to identify haloes, such as their Sunyaev-Zeldovicheffect, X-ray emission, gravitational lensing, or number of (cid:63)
E-mail: [email protected] (LO) † E-mail: [email protected] (REA) optically-detected galaxies. Despite these differences, a nec-essary ingredient for all such analyses is accurate predictionsfor the abundance of haloes of a given mass as a function ofcosmological parameters.In the Press-Schechter formalism (Press & Schechter 1974)(hereafter, PS), the abundance of dark matter haloes of mass M is fundamentally given by the relative abundance of peaksof different types in a Gaussian random field. Specifically, thehalo mass function reads: n ( M ) d log( M ) = − ρ b M d log σd log M νf ( ν ) (1)where ρ b is the background matter density of the Uni-verse; f ( ν ) = (cid:112) /π exp( − . ν ); ν is the so-called “peakheight” associated to a halo of mass M and is defined as ν ≡ δ c ( z ) /σ ( M, z ); δ c is the critical overdensity for collapse;and σ ( M, z ) is the rms linear variance extrapolated at theredshift of interest, z .In this approach, cosmological parameters and the shapeof the power spectrum of fluctuations, P ( k ), are considered © a r X i v : . [ a s t r o - ph . C O ] F e b L. Ondaro-Mallea et al. only through modifications to σ : σ ( R, z ) = D ( z )2 π (cid:90) ∞ d kP ( k ) ˜ W ( k, R ) (2)where D ( z ) is the linear growth factor, and ˜ W ( k, R ) is theFourier transform of the top-hat window function and M = π ρ b R . On top of this, n ( M ) is affected by the cosmologcalparameters through ρ b . Since f ( ν ) is cosmology-independent,the halo mass function is predicted to be “universal”.The “universality” of the mass function is a key propertybecause it allows for accurate predictions even if PS itselfis inaccurate. For instance, if the mass function is universal,a single simulation is needed to measure f ( ν ), and then useEq. 1 to make predictions for any cosmological model. There-fore, computational resources can be focused on accuratelymeasuring f ( ν ) using high force and high mass resolutionsimulations of large cosmic volumes, rather than requiringlarge ensembles of simulations spanning the full range of cos-mological parameters of interest.In fact, several early studies found that the PS halo massfunction describes only qualitatively the abundance of darkmatter haloes in N -body simulations. Motivated by the uni-versality of the mass function, these works have providedmuch more precise fitting functions for f ( ν ), usually em-ploying functional forms inspired by ellipsoidal collapse, butstill assuming that all cosmology dependence can be capturedthrough σ ( M ) (e.g. Sheth & Tormen 1999; Jenkins et al. 2001;Sheth & Tormen 2002; Reed et al. 2003; Warren et al. 2006;Reed et al. 2007; Crocce et al. 2010; Bhattacharya et al. 2011;Angulo et al. 2012; Watson et al. 2013; Bocquet et al. 2016;Seppi et al. 2020).More recently, various authors pointed out and quantifiedthe “non-universality” of the mass function (Tinker et al.2008; Courtin et al. 2010; Despali et al. 2015; Diemer 2020).They have found that the amplitude and shape of f ( ν ) doesdepend on redshift and cosmology in a complicated manner,which depends on the halo definition, and can modify by upto 10% the expected abundance of haloes of a given mass.This can be easily the leading theory systematic error in thecosmological analysis of future cluster catalogues (Artis et al.2021).One of the main goals of this paper is to explore the non-universality of the halo mass function. That is, the depen-dence of the abundance of dark matter haloes on cosmologyand/or redshift in addition to that on the linear rms varianceof fluctuations, σ ( M ). For this, we will consider cosmologieswith identical values for σ ( R ) at z = 0, but with very differentgrowth histories. In this way, any signal of non-universalitycan be attributed to the way in which haloes grow, since thestatistics of the initial fluctuation field are identical. This canshed light on the origin of the mass function non-universalityand allow for a more accurate modelling. In addition, we willconsider simulations with fixed growth history but varyingthe power spectrum of primordial fluctuations.Indeed we will show that by explicitly accounting for thedependence of f ( ν ) on the growth rate and power spectrumslope, we are able to predict the halo mass function with a2-3% accuracy over essentially the whole currently viable cos-mological parameter space, including dynamical dark energy.Moreover, this modelling allows to improve the accuracy withwhich cosmology-rescaling algorithms predict the abundanceof haloes. Figure 1.
Linear properties of the cosmological models we con-sider and simulate throughout this paper.
Left panel:
Growth fac-tor, D ( a ) and growth rate, f ( a ), as a function of the expansionfactor a . Right panel:
Linear mass variance at z = 0 as a functionof the Lagrangian radius of haloes of mass M . The outline of this paper is as follows. In § N -body simula-tion we carry out. In § § § § § In this section we will describe our set of cosmological sim-ulations and our measurements of the halo mass function.Specifically, in § § § We will consider 9 cosmological models given by a combina-tion of 3 different growth histories and 3 linear power spectra.In this way we can explore the effect of the growth historyat a fixed linear mass variance, and of the power spectrumshape at a fixed growth history. We note that, in practice, weobtain varying growth histories by defining them with differ-ent values of the matter density parameter, Ω m , and vary thepower spectrum shape by considering different values of theprimordial spectral index n s (see Table 1).In the left panel of Figure 1 we show the linear growthfactor, D ( a ), and growth rate, f ≡ d log Dd log a , as a function ofexpansion factor a of the models we will consider. By con-struction, at z = 0 all models have the same linear ampli-tude, however, they show very different values for the lineargrowth rate. At one extreme (green lines) we have a cos-mology where structure initially grew very quickly and thenstalled, where we expect very little mass accretion today. Atthe other extreme (blue lines) is a cosmology where structure MNRAS , 1–14 (2021) on-universality of the mass function Table 1.
The cosmological parameters that we vary to obtainthe 9 cosmological models we simulate. We keep the rest of thecosmological parameters fixed assuming flat cosmology and Ω b =0 . σ = 0 . h = 0 . ν = 0 . , w = 0 . w a = 0 . n s m
1. 0.307 0.148Ω Λ
0. 0.693 0.852 has been growing at the same pace through the history of theUniverse, and in particular, we expect it to yield the highestpresent-day accretion rates onto dark matter haloes.In the right panel of Figure 1 we show the 3 different σ ( R )at z = 0 we consider. The respective power spectra are givenby linear predictions for a cosmology consistent with recentobservational constraints (c.f. Table 1), for three different val-ues of the primordial spectral index, n s = { . , . , . } .Although these values are clearly inconsistent with currentdata, they will allow us to clearly identify the role of theshape of fluctuations at a fixed growth history. Specifically,the cosmology with n s = 0 .
75 displays a very flat power spec-trum, thus the density field has more similar fluctuations onall scales. We expect this to yield to similar collapse redshiftsamong different halo masses. On the other hand, the casewith n s = 1 .
25 features stronger small scales fluctuations,thus we expect small haloes collapsing at high redshifts andlarge haloes forming at progressively later time. N -body simulations For each cosmological model described in the previous sub-section, we have carried out a suite of cosmological simu-lations with N = 1024 particles and 4 different box sizes, L = { , , , } h − Mpc. This allows us to com-pute the halo mass function over a broad range of halo masseswith a sufficient statistical accuracy at a moderate computa-tional cost. Therefore, in total we have a suite of 36 simula-tions. The details of the simulations are listed in Table 2.Each of our simulations is initialized at z = 49 usingsecond-order Lagrangian perturbation theory. As recentlypointed out by Michaux et al. (2020), this configuration isexpected to be accurate at the 2% level for the abundancehaloes resolved with more than 100 particles.We carry out our simulations with an updated version ofthe L-Gadget3 code (Angulo et al. 2020), employing 48 MPITasks. In all cases, we set the Plummer-equivalent softeninglength to a 2% of the mean interparticle separation. Eachof our simulations took approximately 1 to 3 thousand CPUhours, depending on the mass resolution of the simulation.
We construct halo catalogues employing a Friends-of-Friends(FoF) algorithm with a linking length parameter b = 0 . z = { , . , } simulation outputs. Additionally, foreach FoF halo we compute the spherical-overdensity masses M ∆ = π ∆ r , for ∆ = { ρ c , ρ b , ∆ vir } , where ρ b is the Table 2.
The main numerical parameters of our simulations. L isthe box size; (cid:15) the gravitational softening length; and m p the massof each N -body particle. L [ h − Mpc] 200 600 1200 2400 (cid:15) [ h − Mpc] 0.004 0.012 0.023 0.047 m p [ h − M (cid:12) / e
10] extreme1 0.207 5.58 44.67 357.3central 0.064 1.71 13.7 109.7extreme2 0.031 0.82 6.62 52.9
Figure 2.
The differential abundance of M haloes in one ofour cosmological model (Ω m = 0 . n s = 0 . z = 0 and z = 1, as estimated in 4 cosmological simulations of various sizes.The middle and bottom panels display the measurements rela-tive to the expectations of the fitting function with dependenceon n eff and α eff developed in this work. The shaded regions corre-spond to the Poisson uncertainty of the measurements. Moreover,for comparison, dotted lines display the fitting functions developedin Despali et al. (2015) and Tinker et al. (2008), as indicated bythe legend. mean matter density of the Universe, and ∆ vir ≡ ρ c { π − − Ω m ( z )] − − Ω m ( z )] } is the virial overdensity ex-pected at each cosmological model. For each mass definition,we compute the halo mass function by considering haloeswith more than 32 particles in equally-spaced logarithmicbins, ∆ log M = 0 . MNRAS000
The differential abundance of M haloes in one ofour cosmological model (Ω m = 0 . n s = 0 . z = 0 and z = 1, as estimated in 4 cosmological simulations of various sizes.The middle and bottom panels display the measurements rela-tive to the expectations of the fitting function with dependenceon n eff and α eff developed in this work. The shaded regions corre-spond to the Poisson uncertainty of the measurements. Moreover,for comparison, dotted lines display the fitting functions developedin Despali et al. (2015) and Tinker et al. (2008), as indicated bythe legend. mean matter density of the Universe, and ∆ vir ≡ ρ c { π − − Ω m ( z )] − − Ω m ( z )] } is the virial overdensity ex-pected at each cosmological model. For each mass definition,we compute the halo mass function by considering haloeswith more than 32 particles in equally-spaced logarithmicbins, ∆ log M = 0 . MNRAS000 , 1–14 (2021)
L. Ondaro-Mallea et al.
Thus, we have followed a conservative approach and imposea cut of 200 particles per halo without any additional correc-tion. This limit, as shown by Ludlow et al. (2019) is enoughto keep all the numerical effects below 5% for all mass defi-nitions. Consequently, we will add in quadrature to Poissonerrors this 5% to account for possible systematic errors in themeasurement of our mass functions.
In order to span a broad mass range we have combined simu-lations of many box sizes. For different box sizes, however, theoutput times can vary slightly since in
L-Gagdet3 we choosethem to coincide with a global timesteps which, in turn, canvary from simulation to simulation. In Appendix B we de-scribe and validate a simple model with which we accountfor this effect in our measurements.In addition, the lack of modes larger than the simulatedbox will induce systematic differences among different boxsizes (e.g. Power & Knebe 2006; Luki´c et al. 2007; Reed et al.2007). However, we have checked that for all the boxes theseeffects are sub-percent at the relevant masses.In Figure 2 we display our measured M halo mass func-tion at z = 0 for the cosmological model with Ω m = 0 . n s = 0 . × h − M (cid:12) up to 10 h − M (cid:12) . The agreement in the overlapping regionsis always better than 5%, consistent with our systematic er-ror estimate. Although not shown here, we have checked thatthis also holds for the other 8 cosmological models.For comparison, in the bottom panels we also display thefitting functions of Despali et al. (2015) and Tinker et al.(2008). Although some differences among our data and thesemodels are expected due to differences in the group finder, thecomparison readily highlights the impact of non-universalityof the mass function. At z = 0 our model and that of Despaliet al. (2015) and Tinker et al. (2008) are in reasonable agree-ment. However, at z = 1, these fits overestimate by morethan 10-15% the abundance of haloes in our simulations. Insubsequent sections we will explore this issue in greater de-tail. In this section we will compare how the same linear fluctu-ation turns into collapsed objects of different mass for dif-ferent cosmologies. We will then explore the dependence ofthe mass function on both growth rate and the slope of thepower spectrum.
Figure 3.
The projected simulated density field normalised by themean background density at z = 0 for two cosmological modelsthat share the same linear density field but differ significantly intheir current growth rate. Top panels show the full simulated box,L=200 h − Mpc, whereas the middle panels and bottom panel showzooms into regions of 30 h − Mpc and 7 h − Mpc a side centered ina halo of normalised mass M /ρ b ∼ . × h Mpc . In theleft column we plot the cosmological simulation with the lowestmatter density Ω m = 0 . m = 1. The solid, dashedand dotted circles show r b , r c and r vir radii of the halo. The universality of the mass functions assumes that the massfunction is completely described by the linear density field.In order to test this assumption, we have run simulationswith very different growth histories that share the same lin-ear density field at z = 0. In Figure 3 we show the simu-lated density field at z = 0 for our L = 200 h − Mpc sim-ulations with the most dissimilar growth histories for the n s = 0 . h − Mpc wide, whereas in the middle and bottom pan-els we zoom on a massive dark matter halo of normalisedmass M /ρ b ∼ . × h Mpc . We can see that, althoughboth cases corresponding to identical z = 0 linear densitypeaks, their nonlinear counterparts are different. Specifically,in the case with the highest Ω m value, and thus, highest cur-rent growth rate (in the right), haloes are significantly lessdense in its center, which is consistent with its expected lower MNRAS , 1–14 (2021) on-universality of the mass function Figure 4.
The profiles of the crossmatched haloes with the samelinear density field for two ν bins. In the first row we display thedensity profiles, in the second row the cumulative mass and in thethird row the logarithmic slope of the density profile. The verticallines indicate the values of r , r vir and r radii. The shadedarea represents r < . (cid:15) , where (cid:15) is the softening length. formation redshift and thus lower concentration parameters.Various definitions of halo radii are displayed by white cir-cles in each case. By comparing r and r vir radii in the twocosmologies, we see that they identify very different regionsof the halo, unlike r which is similar in both cases.To explore this further, we will compare haloes of the samepeak height in different simulations. In particular, we havecross-matched halo catalogues among simulations that sharethe same linear power spectrum at z = 0. For this, we as-sociate two haloes based on their position and peak height.From Figure 3 we expect that the same fluctuation in thelinear density field will end up having a different mass de-pending on its nonlinear evolution.We will characterise each halo by an ”effective growth rate”and an effective ”local power spectrum slope”, which we de-fine respectively as: α eff ( a ) ≡ d log( D ) d log a | a = a ev , (3)where a ev is defined implicitly via D ( a ev ) = γD ( a ) with γ =4 /
5, and
Figure 5.
The mass ratio of the cross-matched haloes at z = 0respect to the self-similar cosmology, i.e. α eff = 1. The coloursdepict the current growth rate value of the given cosmology. Ineach panel we display the ratios for M FoF , M , M vir , and M mass definitions. n eff ≡ − − d log σ ( R ) d log R | κR L ( M ) (4)where κ = 1, and R L is the Lagrangian radius of a haloof mass M . Physically, these two parameters will be captur-ing how quickly haloes have recently grown and the densityprofile of the collapsing region, which can be considered asa proxy for the full mass accretion and merger history of agiven halo.Note that the effective growth rate is not evaluated at theredshift in which we identify a halo, but it is evaluated inthe past, i.e. γ <
1. By this we seek to capture not the rateof current mass accretion, but instead the amount of massthat has been accreted recently. We have tried different defi-nitions of α eff and found that this distinction was particularlyimportant for models with dynamical dark energy. We chosethe numerical value for γ as that which provided the mostaccurate and simplest model for the halo mass functions, aswe will show in Section 6.In Figure 4 we display the spherically averaged mass dis-tribution around crossmatched haloes in 2 bins of the peakheight, ν ∼ ∼
2. We display the average density pro-file, the cumulative mass profile, and the logarithmic slopeof the halo density profiles. Vertical lines indicate the radiusat which the average enclosed density reaches a value equalto 200 times the background, virial and critical density, asindicated in the legend. Coloured lines indicate three growthhistories for the simulations with n s = 0 . α eff , at all values of ν . The higher the growthrate the lower the enclosed mass respect to the backgrounddensity at a given physical radius. We emphasise that all MNRAS000
2. We display the average density pro-file, the cumulative mass profile, and the logarithmic slopeof the halo density profiles. Vertical lines indicate the radiusat which the average enclosed density reaches a value equalto 200 times the background, virial and critical density, asindicated in the legend. Coloured lines indicate three growthhistories for the simulations with n s = 0 . α eff , at all values of ν . The higher the growthrate the lower the enclosed mass respect to the backgrounddensity at a given physical radius. We emphasise that all MNRAS000 , 1–14 (2021)
L. Ondaro-Mallea et al.
Figure 6.
The mass ratio of the cross-matched haloes at z = 0with respect to the cosmological model with self-similar growth,i.e. α eff = 1, and identical linear density field. The colours repre-sent the n eff values of the cross-matched haloes. In each panel wedisplay the ratios for M FoF , M , M vir , and M mass defini-tions. these objects share the same shape and amplitude of theirlinear overdensity field at z = 0. Thus all changes necessarilyare caused by the different growth history.The different growth histories are expected to influencethe internal structure of haloes. In particular, lower growthrates are expected to cause lower current accretion rates ontohaloes, which implies haloes formed earlier and thus are ex-pected to have higher concentrations. In the first row of pan-els we see that this is indeed the case. Inner regions of haloesappear more concentrated. However, the changes are not lim-ited to the concentration, as external parts are also modifiedincreasing their density the higher the current growth rate. Infact, there seems to be an inflection point located at around r radius.Despite the systematic dependence on α eff , the profiles arevery similar when expressed in physical units. However, as aconsequence of the pseudo-evolution of the halo boundaries(Diemer et al. 2013), when expressed in r ∆ units the profilesbecome very different. The pseudo-evolution of the bound-aries is clear in Figure 4. In the α eff = 1 cosmology r isalmost three times larger than in the α eff = 0 .
343 case, while r radii remain roughly constant. Thus, depending on howwe define the boundary of our halo, the mass differences willbe enhanced or suppressed.In the lower panels of Figure 3 this can be appreciated visu-ally. The panels show the most massive crossmatched halo at z = 0 in α eff = 0 .
343 (left) and α eff = 1 (right) cosmologies.The dashed, dotteed and solid lines represent r , r vir and r radii of the halo. While r defines a halo boundaryroughly at the same physical location, r compares verydifferent regions of the density field. This effect is less impor-tant for r vir , which might explain why Despali et al. (2015)found that M vir is the mass definition that leads to the mostuniversal behavior.In order to explore this effect more systematically, in Fig- ure 5 we show how the masses of the crossmatched haloesdiffer depending on the growth history and the mass defini-tion. We display the ratio of the masses of the crossmatchedhaloes respect to the α eff = 1 cosmology as a function of ν .We see that in the cosmology with the lowest current growthrate value M masses are around 30% smaller than in ourreference cosmology. However, for the same cosmologies andhaloes, M masses are around 20% more massive.This effect has two contributions. On the one hand, r radii lie in the inner parts where the effect of the growth his-tory on the mass profile is larger. On the other hand, becauseof the pseudo evolution of r , we compare the masses en-closed in different physical radii. As a consequence, even ifat a given physical radius the enclosed mass is always largerfor haloes in low growth rate cosmologies, when comparing M masses it seems that haloes in high growth rate cos-mologies are more massive. Note that this is solely because wecompare masses enclosed in different physical regions. Thus,the non-universality of M mass function is in a big partdue to the evolution of the boundary itself.Finally, we want to explore the effect of the local slopeof the power spectrum in the crossmatched haloes. For agiven power spectrum, redshift and ν , n eff is completely de-termined. Therefore, in order to see the effect of this variableon the mass of the haloes, we crossmatch the cosmologies with α eff = 0 .
52 and α eff = 1 for the three power spectra definedwith the three n s values we have considered in this work.Note that we only crossmatch cosmologies with the same lin-ear power spectrum. However, if the local power spectrumslope affects the mass of the haloes, we expect the depar-tures of the masses of α eff = 0 .
52 from α eff = 1 cosmologyto be different in the three linear density fields. In Figure 6we show the results following Figure 5, coloured by the n eff values of the crossmatched haloes. Indeed, we see that for agiven ν , the departures of α eff = 0 .
52 haloes from α eff = 1haloes are different depending on the n eff value of the halo.Nevertheless, these differences are much smaller than the dif-ferences that haloes with different α eff show.In summary, the whole density profile of the halo is affectedby the growth history in a non trivial way. This effect will bereflected in the mass function in a different fashion dependingon how masses are defined. Specifically, we expect the non-universality of the mass function to change with the massdefinition. From the numerical simulations described in the previous sec-tion, we have obtained measurements of the halo mass func-tion in 9 cosmologies – 3 growth histories and 3 differentpower spectrum slopes – covering a broad range of masses atvarious redshifts.In order to compare these measurements and estimate theimpact of the non-universality in the mass function, we havecomputed νf ( ν ), where ν is δ c /σ ( M ). Operationally, we firstmeasure the mass function n ( M ) and then estimate νf ( ν )inverting Equation 1.The critical density for collapse, δ c ( z ), can be estimatedas 3 / . π ) / Ω m ( z ) . (Kitayama & Suto 1996), with anexplicit dependence on Ω m . In this work we approximate itwith the value that corresponds to a universe where there isonly matter, δ c = 1 . MNRAS , 1–14 (2021) on-universality of the mass function Figure 7.
The measured M mass functions of the 9 cosmo-logical models considered in this work at z = 0, z = 0 . z = 1.The bottom panel shows the ratio relative to the mean value ineach ν bin. We show in gray the z > z = 0 mass functions according to their α eff value. Theshaded area corresponds to measurements with ν < .
7, which willbe excluded when developing a fitting function for f ( ν ). and cosmology dependence and make it easier to model. Wehave checked that the deviations from universal behaviour ofthe mass function are stronger (up to ± ±
10% (see Appendix A).In Figure 7 we display the measurements of νf ( ν ) from M mass function in the cosmologies listed in the Table 1at z = 0, z = 0 . z = 1. The mass functions with z > z = 0 mass functions arecoloured according to their α eff value. The average value ineach ν bin is displayed as a black solid line, and it is used as areference for the ratio displayed in the bottom panel. To avoidpossible biases due to differential coverage of our models ,we will restrict our subsequent analysis to the range 0 . <ν <
5. For the cosmology most consistent with observationalconstraint, this implies a mass range of 10 < M/ [ h − M (cid:12) ] < at z = 0.In this figure we can clearly see deviations from an uni-versal behavior. For ν values above unity, haloes of a givenpeak height in one cosmology can be up to 70% more abun-dant than in others. By construction, the origin of this non-universal behavior must be in a combination of the differentstatistics of the initial Gaussian random fields and the differ-ent growth histories. Indeed, we can already see that this is At a fixed volume and number of particles, the differences in thepower spectrum shape and Ω m lead to differences in the range of ν -peaks that our simulations are able to resolve. the case for the cosmologies that share the same linear den-sity field at z = 0. At a given ν , haloes seem to be moreabundant the lower the growth rate value. Recall that, as wesaw in the previous section, this is a consequence of the samefluctuation being more massive for low growth rate values.We now explore how these deviations correlate the valueof the effective growth rate and power spectrum slope at anyredshift. In Figure 8 we display the deviations from the mean νf ( ν ), at a fixed ν as a function of n eff and α eff (horizontaland vertical axis). Each panel shows the result for a differentmass definition.For all mass definitions we can see that the non-universalityclearly correlates with these properties – despite them beingonly a proxy of very different merger and assembly histories.In the second panel we see that for M , deviations aroundthe mean can reach ± νf ( ν ) typically have lower growth rate values,whereas those with higher growth rate values lead to lowerabundances. At fixed α eff , deviations from universality aremuch smaller, about 15%, and they correlate with n eff . Notethat here we are plotting measurements of many redshifts,therefore, we expect that the redshift evolution of the massfunction could be described through the dependence on thesephysically-motivated variables.However, the non-universality of the mass function de-pends on the mass definition. Among those considered in thiswork, M ( M FoF ) mass functions are the most (least) non-universal with deviations up to ±
30% ( ± M vir mass functions are the most universal amongthe SO mass functions. Furthermore, the dependence on thegrowth rate is inverted in M and M vir cases with respectto M case (in agreement with Diemer 2020).In summary, haloes of a given peak height, ν , are moreabundant the lower the growth rates and the shallower thepower spectrum slope. In other words, a halo that forms earlyand has grown mostly trough minor mergers, will be moremassive than another that has recently formed and has expe-rienced a lot of major mergers, even if both have an identicalpeak height in the linearly extrapolated initial field.There could be different paths to follow at this point. Onecould be to find the halo boundary definition that minimisesthe non-universality of the mass function. In fact, we haveseen that mass definitions based on the critical density in-duce strong pseudo-evolution in the mass function driven bythe change of the boundary of the halo. One quantity that hasbeen argued is more physical is the turnaround radius, whichby definition encloses the outermost shell that has collapsed.In the same direction, the first explorations of the splash-back mass functions have been done (Diemer 2020). Otheralternatives have been recently proposed, which are claimedseparate better the linear and nonlinear regimes (e.g. Garciaet al. 2020; Fong & Han 2020).Any universal mass definition, as we saw earlier, would cor-respond to very large scales, which although perhaps bettersuited for describing the mass distribution, might not de-scribe equally well, and thus it might display less correlationwith the properties of collapsed gas and of the galaxies hostedby the halo. In addition, many of the proposed halo defini-tions are ambiguous to implement numerically.Another option would be to develop a model for thechanges of the full density profile as a function of the mass ac- MNRAS , 1–14 (2021)
L. Ondaro-Mallea et al.
Figure 8.
The deviation of the mass functions respect to the mean value (computed in each ν bin) plotted according to the α eff and n eff values. The panels correspond to different mass definitions. cretion history. For instance, in an analogous manner to themodels developed for the relationship between the concen-tration and the expected mass accretion history in ExtendedPress Schechter (Ludlow et al. 2016, 2019), it is perhaps pos-sible to develop a model for the outer regions of a halo, whichwould then predict the changes in halo mass at any radius.The option we will follow here is to adopt a standard halodefinition but calibrate the predictions for the halo abun-dance to be a function of the peak height but also of theproperties of the cosmological model. We will show that witha simple parameterization in terms of α eff and n eff , we canaccurately describe the halo mass function for a large regionin cosmological parameters space. In the previous sections we showed how the halo mass func-tion varies systematically with growth rate and slope of thepower spectrum. In this section, we will model this depen-dence explicitly.We will focus on the M mass function, because as ithas been discussed in the previous Section, it is the mostphysically motivated choice and presents the least pseudo-evolution among the overdensity mass definitions. FoF massfunction could also be a good candidate, but the fact that ithas no clear observational counterpart and that the massesare very subject to numerical effects make it a less interest-ing candidate. However, in Appendix C we show that ourapproach is valid to model the mass function of any of theother mass definitions considered in this work.We have employed the following functional form f ( ν, n eff , α eff ) to model each of our measurements f ( ν, n eff , α eff ) = f ( ν ) f ( n eff ) f ( α eff ) (5) f ( ν ) = 2 A mp (1 + ( a ν ) − p ) (cid:114) a ν π e − . a ν (6) f ( n eff ) = n n + n n eff + n (7) f ( α eff ) = a α eff + a (8) where { a, p, A mp , n , n , n , a , a } are the free parameters ofthe model, where we have used the same functional form for f ( ν ) as Despali et al. (2015). Notice that the contributionsof the variables are separable. Thus, in principle one couldcalibrate f ( ν ) separately, or reuse previously ran simulations.In each case, we find the best fitting parameters by mini-mizing the χ of the quantity νf ( ν ). The uncertainty in eachmeasurement is given by the Poisson statistics and a system-atic uncertainty of 5%, as discussed in Section 2. We imposea limit of 200 particles per halo and 400 haloes per mass bin.The minimization is done with the optimize.minimize pack-age of scipy , imposing bounds on the possible values that theparameters may take. Other than our main model (Eq. 5), wehave found the best fit parameters for the functional formsthat depend only on ν (Eq. 6), ν and n eff (Eq. 6 ×
7) and ν and α eff (Eq. 6 × νf ( ν ), in all our simulations at all three redshifts, to theircorresponding predictions of Equation 5 (red). In the left,middle, and right panels we display the residuals as a functionof ν , n eff , and α eff , respectively. In all cases, solid and dashedlines display the mean and the median, whereas the shadedareas denote the regions enclosing 90% of the measurements.Thus, this plot quantifies the overall accuracy of each modelin describing the mass function diversity we measured.We can see that indeed, for the full model, f ( ν, α eff , n eff ),the residuals are smaller than 10% over the whole range ofvalues explored, with no noticeable remaining dependencewith either parameter. For comparison, we display also resid-uals with respect to a version of Equation 6 where we havemeasured their parameters to our whole dataset but onlyadopting dependence with respect to ν (blue). As expected,in this case, the residuals are significantly larger, reachingvariations of ± α eff and n eff is achieved by construction, it is in principle not guaranteedthat the amplitude of these residuals decrease significantly.For instance, the mass function could have shown dependenceon many more details of the assembly history of haloes andthe statistics of peaks than simply on the effective growth rateand power spectrum slope. It is, therefore, remarkable thatthe residuals in the mass function are all contained within aregion of ± MNRAS , 1–14 (2021) on-universality of the mass function Table 3.
A table listing the best fit parameters of our M b fitting functions. a p A mp n n n a a f ( ν ) 0.769 0.0722 0.3173 – – – – – f ( ν ) f ( n eff ) 0.7741 0.1746 0.3038 -0.1912 -0.4211 0.9859 – – f ( ν ) f ( α eff ) 0.772 0.0308 0.3069 – – – -0.6255 1.5654 f ( ν, n eff , α eff ) 0.7691 0.1309 0.3092 -0.1178 -0.3389 0.3022 -1.0785 2.97 Figure 9.
Deviations between νf ( ν ) measured in our N -body sim-ulations and the predictions of the fitting functions developed inthis work. Left, middle, and right panels display these deviationsas a function of the peak height, ν ; the effective power spectrumslope, n eff ; and the effective growth rate, α eff , respectively. In eachpanel, shaded regions indicate the region that contains 90% of oursimulated results when employing a fitting function calibrated only as a function of ν (blue) or additionally including dependence withrespect to n eff (green), α eff (orange), or both of them (red). In the next section we will explore whether our approachis actually able to describe accurately the mass function inmultiple cosmologies currently allowed by observational data.
To asses the accuracy of our description for the halo massfunction, we will compare its predictions against a suite ofsimulations with 30 different cosmologies. Each of our simu-lations evolved 1536 particles inside a box of approximately L = 512 h − Mpc. The initial conditions where created using2nd-order Lagrangian Perturbation theory at z start = 49 andfixing the amplitude of Fourier modes (Angulo & Pontzen2016). The cosmologies were chosen so that they cover a re-gion of approximately 10 σ around Planck’s best fit values.Specifically, they cover the following parameters ranges: σ ∈ [0 . , . m ∈ [0 . , . b ∈ [0 . , . n s ∈ [0 . , .
99] (9) h [100 km s − Mpc − ] ∈ [0 . , . M ν [eV] ∈ [0 . , . w ∈ [ − . , − . w a ∈ [ − . , . M ν , us-ing the linear response approach of Ali-Ha¨ımoud & Bird Figure 10.
Comparison between M halo mass functions inmultiple cosmologies as measured in N -body simulations relativeto the fit developed in this work (solid lines) and the model devel-oped in Despali et al. (2015)(dotted lines), at z = 0. We display R | data , fit = dndlnM (cid:0) datafit (cid:1) . The vertical lines display the limits ofthe bins with at least 400 haloes resolved with more than 200 par-ticles. Each row, from top to bottom, display variations in n s , Ω m ,Ω b and σ for the first column and h , Ω ν , w and w a for thesecond column. The gray lines display the result for the fiducialcosmology. The shaded areas denote regions of ±
5% and ± (2013); and dynamical dark energy with an equation of state w ( z ) = w + (1 + z ) w a . The cosmology of each simulation isobtained by changing one cosmological parameter of a fiducialcosmology while keeping the rest fixed. The fiducial cosmol-ogy assumes flat geometry, massless neutrinos ( M ν = 0), adark energy equation of state with w = − w a = 0,an amplitude of matter fluctuations σ = 0 .
9, cold dark mat-ter density Ω cdm = 0 . b = 0 .
05, andnormalised Hubble constant h = 0 . MNRAS000
05, andnormalised Hubble constant h = 0 . MNRAS000 , 1–14 (2021) L. Ondaro-Mallea et al.
Figure 11.
Same as Figure 10 at z = 1. show the limits of the bins with at least 400 haloes with re-solved with more than 200 particles as vertical lines.In each row we show the mass functions of the cosmologieswhere we vary one cosmological parameter keeping the restfixed. In Figure 10 we display the results at z = 0, and inFigure 11 the results at z = 1. For , we also show the resid-uals respect the model developed in Despali et al. (2015) asdotted lines, which assumes universality of the mass function.We recall that the functional dependence on f ( ν ) is the samein both models, while in our model we have added extra de-pendences on n eff and α eff in order to capture the effect ofgrowth history on the mass function.At z = 0, our model describes the low mass end of the massfunction at an accuracy of 3%, while Despali et al. (2015)predicts that haloes are 10% more abundant. This may bea consequence of using different group finders, SO and FoFrespectively. For haloes with masses above M > h − M (cid:12) ,there seems to be an underprediction of our fitting function.To investigate this, we have compared our predictions againstthe simulations of Angulo et al. (2020), which feature thesame mass resolution as our test suite but on a volume 27times larger. Although not shown here, in such case we findan agreement to better than 5% up to 10 h − M (cid:12) . Combinedwith the good agreement of our predictions with those ofDespali et al. (2015), we speculate that there is a systematicover prediction of the abundance of haloes in our test simsfor M > h − M (cid:12) , which could be caused by finite-volumeeffects.At z = 1, the redshift evolution of the mass function isevident. Even if at z = 0 Despali et al. (2015) is a good description to the mass function, at z = 1 it overpredictsthe abundances for more than 10%. Our model captures thisand yields results that are accurate within 5% at all massesconsidered.Compared to the redshift evolution, the cosmology depen-dence of the mass function seems to be weak. However, themass functions of w cosmologies present strong deviationsfrom universality. The scatter of the ratio with respect toDespali et al. (2015) is of ∼
5% among cosmologies with dif-ferent w values at both redshifts. After taking into accountthe dependences on α eff and n eff , this scatter vanishes. Weemphasise that we have only used ΛCDM cosmolgies to cal-ibrate the fit, and so α eff and n eff are physically meaningfulproxies of the non-universality of the mass function.It is interesting to notice that the largest part of the im-provement is obtained when adding α eff to the universal de-scription. This is expected, because as discussed in Section3, the deviations from universality correlate much strongerwith α eff than with n eff . To approach an optimal exploitation of current and futureobservations of the abundance of dark matter haloes and theclustering of galaxies, very accurate theoretical predictionsfor these quantities are required. Although fitting functionsand calibrated recipes are indeed extremely valuable, they fallshort in providing correlations among different observables orthe full three dimensional distribution of clusters of galaxies.One option to obtain those predictions is to employ cosmo-logical N -body simulations together with cosmology rescalingalgorithms.The basic idea of cosmology rescaling is to employ a fewsimulations carried out adopting specific cosmological param-eters, and then manipulate their outputs so that representnonlinear structure in any other set of cosmologies. Thesealgorithms have been extensively discussed and tested in An-gulo & White (2010); Angulo & Hilbert (2015); Ruiz et al.(2011); Mead & Peacock (2014b,a); Renneby et al. (2018). Inparticular, Zennaro et al. (2019) showed these are applica-ble to cases of massive neutrinos, and Contreras et al. (2020)showed that the clustering of dark matter and dark matterhaloes and subhaloes can be obtained to better than 3% accu-racy from large to very small scales (0 . < k/h Mpc − < MNRAS , 1–14 (2021) on-universality of the mass function Figure 12.
Comparison between M halo mass functions inmultiple cosmologies as measured in N -body simulations relativeto that in cosmology-rescaled simulations at z = 0. Specifically, wedisplay R | scaled , target = dndlnM (cid:0) scaledtarget (cid:1) .The vertical lines displaythe limits of the bins with at least 400 haloes resolved with morethan 200 particles. Each row, from top to bottom, display varia-tions in n s , Ω m , Ω b and σ for the first column and h , Ω ν , w and w a for the second column. In each panel we show results before(dotted lines) and after (solid lines) applying our additional correc-tion accounting for dependence on growth history, as indicated bythe legend (see text for details). The shaded regions denote regionsof ±
5% and ± by minimizing the difference in the linear mass variance inthe target and rescaled cosmologies. In this operation thenumber of particles in each halo is left invariant. Using thisrecipe the scaling algorithm sets the same linear density fieldin the rescaled and target simulations, which is equivalentto assuming the universality of the mass function. However,in this work we have shown that the mass function dependsnot only on the linear density field but also on the entiregrowth history. In fact, we can see this effect in the dottedlines of Figure 12, where we compare the mass functions fromthe rescaled simulations with the target ones. By assuminguniversality of the mass function, the rescaled M massfunctions differ from the target mass functions up to 10% insome cosmologies at z = 0.Our model for the dependence of the halo mass functionon growth history gives us the possibility to construct anadditional correction for cosmology rescaling by taking intoaccount the different growth histories the target and rescaledcosmologies have gone through. Specifically, we first predict Figure 13.
Same as Fig.12 but at z = 1. the rescaled and target mass functions. In the target cosmol-ogy case, the prediction is straightforward. In the rescaledoriginal cosmology case, we compute the expected mass func-tion of the original cosmology once we have applied the cor-responding mass and length scalings, i.e, once we have set thelinear density field equal to the target cosmology’s linear den-sity field. As we have seen in Section 3.1, the growth-historyaffects the masses of the haloes, not the abundance of them.Thus, the difference between target and rescaled mass func-tions is given by a change in mass, which we find by mappingthe rescaled masses to the target masses where the abun-dances are the same. Therefore, for a given pair of original-target cosmologies, we can predict a halo-by-halo mass cor-rection that describes the effect of the non-universality of thehalo mass function.We show the results after applying this correction as solidlines in Figure 12. We can see that in all cosmologies, theaccuracy of the predictions improve in a clear manner. At z = 0 the differences are in most of the cases smaller than2%. Note that our model, calibrated on simulations where weonly vary Ω m and n s , is able to capture the non-universalityof general cosmologies, even in beyond ΛCDM cosmologieswith massive neutrinos and dynamical dark energy included.As seen in Figure 13, at z = 1 the accuracy is as good as for z = 0, reaching ± −
2% over the full range of masses wherewe can measure the mass function accurately.
MNRAS000
MNRAS000 , 1–14 (2021) L. Ondaro-Mallea et al.
In this paper we have studied the non-universality of thehalo mass function. We have run simulations with very ex-treme cosmologies to maximise the deviations from universalbehaviour and we have shown that the halo masses are af-fected by the entire growth history. As a consequence, giventhe same linear density field in two different cosmologies, thehalo mass functions are different.In order to shed light in the origin of the non-universalityof the mass functions, we have crossmatched haloes of differ-ent cosmologies that share the same linear density field andwe have compared their density profiles. Generally, we haveobserved that all the density profiles up to very large radiiare affected by the growth history of the haloes. Further-more, the physical boundaries of haloes selected with densitycriteria are subject to pseudo evolution, and correspond todifferent physical radii for different cosmologies. This effectis more pronounced for overdensities defined with respect tothe critical density of the universe. Therefore, different massdefinitions yield mass functions with different dependenceson redshift and cosmology.We have modelled the non-universality of the mass func-tion adding two additional parameters other than the peak-height ν : the effective growth rate, α eff , and the local slopeof the power spectrum, n eff . Using a total of 8 free parame-ters, our model captures the non-universality and can lowerthe scatter on the halo mass functions in all the cosmologiesconsidered from ±
20% to about ±
5% up to z = 1. In the liter-ature, the redshift evolution of the mass function is typicallyparametrised explicitly within a fiducial cosmology (see e.g.Tinker et al. 2008). On the contrary, here we have modelledsimultaneously the cosmology and redshift non-universalityof the mass function by using physically motivated parame-ters.We have tested our model on an independent set of simula-tions of 30 different cosmologies, including massive neutrinosand dynamical dark energy. By considering the α eff and n eff dependences, we have been able to reproduce the halo massfunctions within a 5% accuracy in all the cosmologies up to M ∼ × h − M (cid:12) until z = 1. We emphasise that thesimulations that we have used to calibrate the models havebeen run with ΛCDM cosmologies. Thus, it is not a trivialresult that our model is able to describe the halo mass func-tions within, for instance, cosmologies that include massiveneutrinos or dynamical dark energy.As an application of our model, we have applied it togetherwith the cosmology rescaling method presented in Angulo& White (2010). We have found that the accuracy in thescaling of the halo mass function improves from 10% to 2% inall cosmologies including dark energy and massive neutrinos,mostly because of the dependency on the growth rate.There are many paths that we would like to explore in fu-ture works. It is well known that baryonic processes alter innon-trivial way the halo mass function. In particular, astro-physical feedback ejects a large amount of gas outside thehaloes boundaries, and therefore haloes result less massive,even by when including a baryonic modelling (see e.g. Castroet al. 2021; Debackere et al. 2020). We plan to extend our for-malism to include the effect of baryons on the halo mass func-tion, by using the so-called baryonification technique (Schnei-der & Teyssier 2015; Aric`o et al. 2019). By combining it with cosmological-rescaling algorithms, we will construct an emu-lator of the halo mass functions, as a function of cosmologicaland astrophysical parameters.In the near future, more precise and accurate predictionsof the halo mass function will be necessary in order to fullyexploit the data of the future surveys. As an example, Artiset al. (2021) estimated that, only considering the precisionof the parameters of the fitting functions in the analysis (i.eassuming universality of the mass function), an improvementfrom 30% to 70% is required. This framework provides uswith a very accurate fit of the halo mass function, which canbe eventually exploited to directly compare against observedclusters count, from optical, X-ray or Sunyaev-Zel’dovich sur-veys. Thus, we anticipate that this model will be of greatvalue value in constraining the cosmological parameters ofthe Universe. ACKNOWLEDGMENTS
LO acknowledges the Summer Internship Program of theDonostia International Physics Center. The authors acknowl-edge the support of the ERC-StG number 716151 (BACCO).SC acknowledges the support of the “Juan de la Cierva For-maci´on” fellowship (FJCI-2017-33816). The authors acknowl-edge computing resources at MareNostrum-IV and the tech-nical support provided by Barcelona Supercomputing Cen-ter (RES-AECT-2019-2-0012, RES-AECT-2020-3-0014). Theauthors thank Jens St¨uecker for the visualization routine usedin this work.
DATA AVAILABILITY
The data underlying this article will be shared on reasonablerequest to the corresponding author.
REFERENCES
Abazajian K., et al., 2019, arXiv e-prints, p. arXiv:1907.04473Ade P., et al., 2019, J. Cosmology Astropart. Phys., 2019, 056Ali-Ha¨ımoud Y., Bird S., 2013, MNRAS, 428, 3375Angulo R. E., Hilbert S., 2015, MNRAS, 448, 364Angulo R. E., Pontzen A., 2016, MNRAS, 462, L1Angulo R. E., White S. D. M., 2010, MNRAS, 405, 143Angulo R. E., Springel V., White S. D. M., Jenkins A., BaughC. M., Frenk C. S., 2012, MNRAS, 426, 2046Angulo R. E., Zennaro M., Contreras S., Aric`o G., Pellejero-Iba˜nezM., St¨ucker J., 2020Aric`o G., Angulo R. E., Hern´andez-Monteagudo C., ContrerasS., Zennaro M., Pellejero-Iba˜nez M., Rosas-Guevara Y., 2019,arXiv e-prints, p. arXiv:1911.08471Aric`o G., Angulo R. E., Contreras S., Ondaro-Mallea L.,Pellejero-Iba˜nez M., Zennaro M., 2020, arXiv e-prints, p.arXiv:2011.15018Artis E., Melin J.-B., Bartlett J. G., Murray C., 2021, Impact ofthe calibration of the Halo Mass Function on galaxy clusternumber count cosmology ( arXiv:2101.02501 )Bhattacharya S., Heitmann K., White M., Luki´c Z., Wagner C.,Habib S., 2011, ApJ, 732, 122Bocquet S., Saro A., Dolag K., Mohr J. J., 2016, MNRAS, 456,2361Bonoli S., et al., 2020, arXiv e-prints, p. arXiv:2007.01910MNRAS , 1–14 (2021) on-universality of the mass function Castro T., Borgani S., Dolag K., Marra V., Quartin M., Saro A.,Sefusatti E., 2021, MNRAS, 500, 2316Contreras S., Zennaro R. E. A. M., Aric´o G., Pellejero-Iba˜nez M.,2020, arXiv e-prints, p. arXiv:2001.03176Courtin J., Rasera Y., Alimi J.-M., Corasaniti P.-S., Boucher V.,F¨uzfa A., 2010, Monthly Notices of the Royal AstronomicalSociety, pp no–noCrocce M., Fosalba P., Castander F. J., Gazta˜naga E., 2010, MN-RAS, 403, 1353Debackere S. N. B., Schaye J., Hoekstra H., 2020, MNRAS, 492,2285Despali G., Giocoli C., Angulo R. E., Tormen G., Sheth R. K.,Baso G., Moscardini L., 2015, Monthly Notices of the RoyalAstronomical Society, 456, 2486Diemer B., 2020, The Astrophysical Journal, 903, 87Diemer B., More S., Kravtsov A. V., 2013, The Astrophysical Jour-nal, 766, 25Fong M., Han J., 2020, arXiv e-prints, p. arXiv:2008.03477Garcia R., Rozo E., Becker M. R., More S., 2020, arXiv e-prints,p. arXiv:2006.12751Hofmann F., et al., 2017, A&A, 606, A118Ivezi´c ˇZ., et al., 2019, ApJ, 873, 111Jenkins A., Frenk C. S., White S. D. M., Colberg J. M., Cole S.,Evrard A. E., Couchman H. M. P., Yoshida N., 2001, MNRAS,321, 372Kitayama T., Suto Y., 1996, The Astrophysical Journal, 469, 480Leroy M., Garrison L., Eisenstein D., Joyce M., Maleubre S., 2021,MNRAS, 501, 5064Ludlow A. D., Bose S., Angulo R. E., Wang L., Hellwing W. A.,Navarro J. F., Cole S., Frenk C. S., 2016, MNRAS, 460, 1214Ludlow A. D., Schaye J., Bower R., 2019, MNRAS, 488, 3663Luki´c Z., Heitmann K., Habib S., Bashinsky S., Ricker P. M., 2007,ApJ, 671, 1160Luki´c Z., Reed D., Habib S., Heitmann K., 2009, The AstrophysicalJournal, 692, 217Mead A. J., Peacock J. A., 2014a, MNRAS, 440, 1233Mead A. J., Peacock J. A., 2014b, MNRAS, 445, 3453Michaux M., Hahn O., Rampf C., Angulo R. E., 2020, MonthlyNotices of the Royal Astronomical Society, 500, 663683Mo H., van den Bosch F., White S., 2010, Galaxy For-mation and Evolution. Cambridge University Press,doi:10.1017/CBO9780511807244More S., Kravtsov A. V., Dalal N., Gottl¨ober S., 2011, The Astro-physical Journal, 195, 4Power C., Knebe A., 2006, MNRAS, 370, 691Press W. H., Schechter P., 1974, ApJ, 187, 425Reed D., Gardner J., Quinn T., Stadel J., Fardal M., Lake G.,Governato F., 2003, MNRAS, 346, 565Reed D. S., Bower R., Frenk C. S., Jenkins A., Theuns T., 2007,MNRAS, 374, 2Renneby M., Hilbert S., Angulo R. E., 2018, MNRAS, 479, 1100Ruiz A. N., Padilla N. D., Dom´ınguez M. J., Cora S. A., 2011,MNRAS, 418, 2422Sartoris B., et al., 2016, MNRAS, 459, 1764Schneider A., Teyssier R., 2015, J. Cosmology Astropart. Phys.,2015, 049Seppi R., et al., 2020, arXiv e-prints, p. arXiv:2008.03179Sheth R. K., Tormen G., 1999, MNRAS, 308, 119Sheth R. K., Tormen G., 2002, MNRAS, 329, 61Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M.,Yepes G., Gottl¨ober S., Holz D. E., 2008, ApJ, 688, 709Warren M. S., Abazajian K., Holz D. E., Teodoro L., 2006, TheAstrophysical Journal, 646, 881Watson W. A., Iliev I. T., D’Aloisio A., Knebe A., Shapiro P. R.,Yepes G., 2013, MNRAS, 433, 1230Weinberg D. H., Mortonson M. J., Eisenstein D. J., Hirata C.,Riess A. G., Rozo E., 2013, Physics Reports, 530, 87
Figure A1.
The difference of the deviations of the mass func-tion respect to the mean computed in each ν bin between themass functions obtained with the critical density for collapse pre-sented in (Kitayama & Suto 1996) and with the critical densitycorresponding to a universe with only matter. ∆ Kitayama , EdS =( f ( ν ) Kitayama − f ( ν ) EdS ) /f ( ν ) EdS where f ( ν ) = νf ( ν ) / <νf ( ν ) > . Figure A2.
Same as Figure A1 but for the critical density com-puted following (Mo et al. 2010).Zennaro M., Angulo R. E., Aric`o G., Contreras S., Pellejero-Ib´a˜nezM., 2019, MNRAS, 489, 5938Zennaro M., Angulo R. E., Pellejero-Ib´a˜nez M., St¨ucker J., Contr-eras S., Aric`o G., 2021, arXiv e-prints, p. arXiv:2101.12187
APPENDIX A: COSMOLOGY AND REDSHIFTDEPENDENT CRITICAL DENSITY
In this appendix we show the effect of taking into accountthe redshift dependence of the critical density for collapseon the non-universality of the mass function. Specifically, wecompute the relative difference of the deviations in each ν bin between the mass functions with δ c ( z ) and δ c = 1 . δ c ( z ). We see that,for most of the cases they do not exceed the ± MNRAS000
In this appendix we show the effect of taking into accountthe redshift dependence of the critical density for collapseon the non-universality of the mass function. Specifically, wecompute the relative difference of the deviations in each ν bin between the mass functions with δ c ( z ) and δ c = 1 . δ c ( z ). We see that,for most of the cases they do not exceed the ± MNRAS000 , 1–14 (2021) L. Ondaro-Mallea et al.
Figure B1.
The ratio of the differential mass function betweentwo expansion factors.
Left panel:
Predicted ratios for the outputexpansion factors corresponding to different box sizes of our sim-ulation set. The reference redshift is z ref = 0 . Right panel:
Predicted and measured ratios for the expansion factors listed inthe legend. The reference redshift is z ref = 1 . deviations around the mean are much stronger, as we can seein Figure 8. APPENDIX B: REDSHIFT CORRECTION
For different box sizes the redshifts of the snapshots varyslightly. This effect is more pronounced at high redshift,where the difference of the output redshifts of different boxescan reach ∆ z ∼ .
01. In this time lapse the mass functionsmay have evolved, therefore, when combining different boxsizes we may be introducing some bias in our data set. Theleft panel of Figure B1 displays the expected ratio of the dif-ferential mass functions between the output redshifts of thedifferent boxes around z = 1. At M ∼ h − M (cid:12) , the dif-ferences can reach 10%.In order to test whether the predicted evolution of the massfunction is accurate, we make use of a simulation presentedin Section 6, for which we have many snapshots. In the rightpanel of Figure B1 we display the predicted and measuredratios for the redshifts listed in the legend. We can see thatthe predictions are a good description of the data. For othercosmologies the results are similar.Thus, we proceed to correct the mass functions of the bigbox sizes in the following way, dnd ln M | corrected ( z ref ) = dnd ln M | measured ( z ) × f ( M, z ref ) f ( M, z ) , (B1)where f ( M, z ) is some model for the differential mass func-tion.
APPENDIX C: EXTENSION TO OTHER MASSDEFINITIONS
In this appendix we present the main results of our modelingwith other mass definitions. These results are analogous towhat already described for M mass functions.In general, as seen in Figure 8 all mass functions showclear correlations with n eff and α eff for a given ν . Therefore,we keep the quadratic and linear functional forms for n eff and α eff in our model (Eq. 7 and 8), but slightly changethe functional form of the peak-height dependence for thedifferent mass definition. For M ∆ mass functions we use the Figure C1.
Ratio of the measured mass functions respect to ourmodel. The ratios are displayed against ν , n eff and α eff values inthe columns. In the rows we show the results for M FoF , M and M vir mass functions. functional form used in Despali et al. (2015) (Eq. 6). How-ever, the functional form used in Angulo et al. (2012) is moresuited to describe M FoF mass functions. Hence, for this massdefinition we replace f with˜ f ( ν ) = A ( b ν c + 1) exp( − d ν ) (C1)where { A, b, c, d } are the free parameters of the ν dependenceof our model. By applying the methodology explained in Sec-tion 4, we have obtained the best-fit parameters listed in Ta-bles C1, C2 and C3 for M FoF , M and M vir mass functionsrespectively.In Figure C1 we show the performance of the model foreach mass definition. For M FoF mass functions the residualscatter is consistent with the intrinsic uncertainties of our fit.However, even if the scatter is reduced significantly, for M vir and M the description is not as good as in the other cases.We think that, for both the cases, this may be a consequenceof the pseudo-evolution of the boundary definition. This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS , 1–14 (2021) on-universality of the mass function Table C1.
A table listing the best fit parameters of our M FoF fitting functions.
A b c d n n n a a f ( ν ) 0.231 1.6871 1.7239 1.1092 – – – – – f ( ν ) f ( n eff ) 0.2297 1.6824 1.6437 1.0975 -0.1565 -0.4757 0.6947 – – f ( ν ) f ( α eff ) 0.2218 1.8171 1.6643 1.1009 – – – -0.2908 1.2186 f ( ν, n eff , α eff ) 0.2276 1.7692 1.6249 1.09 -0.1398 -0.473 0.3671 -0.3715 1.6164 Table C2.
A table listing the best fit parameters of our M fitting functions. a p A mp n n n a a f ( ν ) 0.833 0.1753 0.263 – – – – – f ( ν ) f ( n eff ) 0.769 0.2936 0.2314 -0.979 -3.7864 -2.4854 – – f ( ν ) f ( α eff ) 0.8186 0.1796 0.2541 – – – f ( ν, n eff , α eff ) 0.7957 0.3069 0.2549 -0.502 -1.79 -0.8305 1.8695 -0.0937 Table C3.
A table listing the best fit parameters of our M vir fitting functions. a p A mp n n n a a f ( ν ) 0.7814 0.0854 0.3001 – – – – – f ( ν ) f ( n eff ) 0.7632 0.1867 0.2939 -0.4999 -1.761 -0.4741 – – f ( ν ) f ( α eff ) 0.7793 0.0981 0.2951 – – – 0.3581 0.7179 f ( ν, n eff , α eff ) 0.7693 0.2074 0.2861 -1.0439 -3.4809 -0.3037 0.1721 0.2899MNRAS000