{'fraction_non_alphanumeric': 0.029589280291527482, 'fraction_numerical': 0.017233525660491952, 'mean_word_length': 4.504492269118262, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-universality of the mass function: dependence on the growth rate and power spectrum shape

The abundance of dark matter haloes is one of the key probes of the growth of structure and expansion history of the Universe. Theoretical predictions for this quantity usually assume that, when expressed in a certain form, it depends only on the mass variance of the linear density field. However, cosmological simulations have revealed that this assumption breaks, leading to 10-20% systematic effects. In this paper we employ a specially-designed suite of simulations to further investigate this problem. Specifically, we carry out cosmological N-body simulations where we systematically vary growth history at a fixed linear density field, or vary the power spectrum shape at a fixed growth history. We show that the halo mass function generically depends on these quantities, thus showing a clear signal of non-universality. Most of this effect can be traced back to the way in which the same linear fluctuation grows differently into the nonlinear regime depending on details of its assembly history. With these results, we propose a parameterization with explicit dependence on the linear growth rate and power spectrum shape. Using an independent suite of simulations, we show that this fitting function accurately captures the mass function of haloes over 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.

optically-detected galaxies. Despite these differences, a necessary ingredient for all such analyses is accurate predictions for the abundance of haloes of a given mass as a function of cosmological 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 peaks of different types in a Gaussian random field. Specifically, the halo mass function reads: where ρ b is the background matter density of the Universe; f (ν) = 2/π exp(−0.5ν 2 ); ν is the so-called "peak height" 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 the redshift of interest, z.
In this approach, cosmological parameters and the shape of the power spectrum of fluctuations, P (k), are considered only through modifications to σ: where D(z) is the linear growth factor, andW (k, R) is the Fourier transform of the top-hat window function and M = 4π 3 ρ b R 3 . On top of this, n(M ) is affected by the cosmologcal parameters 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 property because it allows for accurate predictions even if PS itself is inaccurate. For instance, if the mass function is universal, a single simulation is needed to measure f (ν), and then use Eq. 1 to make predictions for any cosmological model. Therefore, computational resources can be focused on accurately measuring f (ν) using high force and high mass resolution simulations of large cosmic volumes, rather than requiring large ensembles of simulations spanning the full range of cosmological parameters of interest.
In fact, several early studies found that the PS halo mass function describes only qualitatively the abundance of dark matter haloes in N -body simulations. Motivated by the universality of the mass function, these works have provided much more precise fitting functions for f (ν), usually employing functional forms inspired by ellipsoidal collapse, but still assuming that all cosmology dependence can be captured through σ(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 quantified the "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 (ν) does depend on redshift and cosmology in a complicated manner, which depends on the halo definition, and can modify by up to 10% the expected abundance of haloes of a given mass. This can be easily the leading theory systematic error in the cosmological analysis of future cluster catalogues (Artis et al. 2021).
One of the main goals of this paper is to explore the nonuniversality of the halo mass function. That is, the dependence of the abundance of dark matter haloes on cosmology and/or redshift in addition to that on the linear rms variance of fluctuations, σ(M ). For this, we will consider cosmologies with identical values for σ(R) at z = 0, but with very different growth histories. In this way, any signal of non-universality can be attributed to the way in which haloes grow, since the statistics of the initial fluctuation field are identical. This can shed light on the origin of the mass function non-universality and allow for a more accurate modelling. In addition, we will consider simulations with fixed growth history but varying the power spectrum of primordial fluctuations.
Indeed we will show that by explicitly accounting for the dependence of f (ν) on the growth rate and power spectrum slope, we are able to predict the halo mass function with a 2-3% accuracy over essentially the whole currently viable cosmological parameter space, including dynamical dark energy. Moreover, this modelling allows to improve the accuracy with which cosmology-rescaling algorithms predict the abundance of haloes. The outline of this paper is as follows. In §2 we describe the cosmological models and the respective N -body simulation we carry out. In §3 we measure the non-universality of our simulated halo catalogues and illustrate its physical origin by comparing haloes across simulations. In §4 we model the departures from universality as a function of an effective growth rate and power spectrum slope in each cosmological model. In §5 we validate these predictions against the halo mass function as measured in a suite of simulations spanning a broad range of cosmological parameters. We further show in §6 that our proposed model can be employed to improve the accuracy of cosmology-rescaling techniques. Finally, we conclude and summarise our findings in §7.

NUMERICAL SIMULATIONS
In this section we will describe our set of cosmological simulations and our measurements of the halo mass function. Specifically, in §2.1 we describe the cosmological models we consider and in §2.2 the numerical setup of the respective simulations. In §2.3 we discuss our measurements of the halo mass function, and how we account for numerical and discretization errors.

Cosmological Models
We will consider 9 cosmological models given by a combination of 3 different growth histories and 3 linear power spectra. In this way we can explore the effect of the growth history at a fixed linear mass variance, and of the power spectrum shape at a fixed growth history. We note that, in practice, we obtain varying growth histories by defining them with different values of the matter density parameter, Ωm, and vary the power spectrum shape by considering different values of the primordial spectral index ns (see Table 1).
In the left panel of Figure 1 we show the linear growth factor, D(a), and growth rate, f ≡ d log D d log a , as a function of expansion factor a of the models we will consider. By construction, at z = 0 all models have the same linear amplitude, however, they show very different values for the linear growth rate. At one extreme (green lines) we have a cosmology where structure initially grew very quickly and then stalled, where we expect very little mass accretion today. At the other extreme (blue lines) is a cosmology where structure has been growing at the same pace through the history of the Universe, and in particular, we expect it to yield the highest present-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 given by linear predictions for a cosmology consistent with recent observational constraints (c.f. Table 1), for three different values of the primordial spectral index, ns = {0.75, 0.96, 1.25}. Although these values are clearly inconsistent with current data, they will allow us to clearly identify the role of the shape of fluctuations at a fixed growth history. Specifically, the cosmology with ns = 0.75 displays a very flat power spectrum, thus the density field has more similar fluctuations on all scales. We expect this to yield to similar collapse redshifts among different halo masses. On the other hand, the case with ns = 1.25 features stronger small scales fluctuations, thus we expect small haloes collapsing at high redshifts and large haloes forming at progressively later time.

N -body simulations
For each cosmological model described in the previous subsection, we have carried out a suite of cosmological simulations with N = 1024 3 particles and 4 different box sizes, L = {200, 600, 1200, 2400} h −1 Mpc. This allows us to compute the halo mass function over a broad range of halo masses with a sufficient statistical accuracy at a moderate computational cost. Therefore, in total we have a suite of 36 simulations. The details of the simulations are listed in Table 2.
Each of our simulations is initialized at z = 49 using second-order Lagrangian perturbation theory. As recently pointed out by Michaux et al. (2020), this configuration is expected to be accurate at the 2% level for the abundance haloes resolved with more than 100 particles.
We carry out our simulations with an updated version of the L-Gadget3 code , employing 48 MPI Tasks. In all cases, we set the Plummer-equivalent softening length to a 2% of the mean interparticle separation. Each of our simulations took approximately 1 to 3 thousand CPU hours, depending on the mass resolution of the simulation.

Halo Catalogues and discreteness correction
We construct halo catalogues employing a Friends-of-Friends (FoF) algorithm with a linking length parameter b = 0.2 at the z = {0, 0.5, 1} simulation outputs. Additionally, for each FoF halo we compute the spherical-overdensity masses  The middle and bottom panels display the measurements relative to the expectations of the fitting function with dependence on n eff and α eff developed in this work. The shaded regions correspond to the Poisson uncertainty of the measurements. Moreover, for comparison, dotted lines display the fitting functions developed in Despali et al. (2015) and Tinker et al. (2008), as indicated by the legend.
mean matter density of the Universe, and ∆vir ≡ ρc{18π 2 − 82[1 − Ωm(z)] − 39[1 − Ωm(z)] 2 } is the virial overdensity expected at each cosmological model. For each mass definition, we compute the halo mass function by considering haloes with more than 32 particles in equally-spaced logarithmic bins, ∆ log M = 0.155. It is known that the FoF algorithm suffers from effects related to particle discreteness, which leads to an overestimation of the mass function (see e.g. Leroy et al. 2021, and references therein). In Warren et al. (2006) an empirical formula was derived to correct for these effects. In agreement with Lukić et al. (2009) and More et al. (2011), we have, however, found that the performance of this correction varied greatly with cosmology and redshift. In addition, finite numerical precision in the force calculation and time integration, as well as the effect of softening length, also affect the abundance of haloes detected by FoF (Ludlow et al. 2019). Thus, we have followed a conservative approach and impose a cut of 200 particles per halo without any additional correction. This limit, as shown by Ludlow et al. (2019) is enough to keep all the numerical effects below 5% for all mass definitions. Consequently, we will add in quadrature to Poisson errors this 5% to account for possible systematic errors in the measurement of our mass functions.

Finite volume and output redshift corrections
In order to span a broad mass range we have combined simulations of many box sizes. For different box sizes, however, the output times can vary slightly since in L-Gagdet3 we choose them to coincide with a global timesteps which, in turn, can vary from simulation to simulation. In Appendix B we describe and validate a simple model with which we account for this effect in our measurements.
In addition, the lack of modes larger than the simulated box will induce systematic differences among different box sizes (e.g. Power & Knebe 2006;Lukić et al. 2007;Reed et al. 2007). However, we have checked that for all the boxes these effects are sub-percent at the relevant masses.
In Figure 2 we display our measured M 200b halo mass function at z = 0 for the cosmological model with Ωm = 0.307 and ns = 0.9611. Results from simulations of various box sizes (after the corrections described above) are denoted by different line colours, as indicated by the legend.
The top panel displays the differential mass function whereas the middle and bottom panels display the ratio with respect to the predictions of a fitting function we will develop later in this work. Note we only display bins with more than 400 objects resolved with at least 200 particles. We can see how our suite of various box sizes complement each other to cover a very large range of halo masses, from 5 × 10 11 h −1 M up to 10 15 h −1 M . The agreement in the overlapping regions is always better than 5%, consistent with our systematic error estimate. Although not shown here, we have checked that this also holds for the other 8 cosmological models.
For comparison, in the bottom panels we also display the fitting functions of Despali et al. (2015) and Tinker et al. (2008). Although some differences among our data and these models are expected due to differences in the group finder, the comparison readily highlights the impact of non-universality of the mass function. At z = 0 our model and that of Despali et al. (2015) and Tinker et al. (2008) are in reasonable agreement. However, at z = 1, these fits overestimate by more than 10-15% the abundance of haloes in our simulations. In subsequent sections we will explore this issue in greater detail.

THE DEPENDENCE OF THE MASS FUNCTION ON GROWTH HISTORY
In this section we will compare how the same linear fluctuation turns into collapsed objects of different mass for different cosmologies. We will then explore the dependence of the mass function on both growth rate and the slope of the power spectrum. . The projected simulated density field normalised by the mean background density at z = 0 for two cosmological models that share the same linear density field but differ significantly in their current growth rate. Top panels show the full simulated box, L=200h −1 Mpc, whereas the middle panels and bottom panel show zooms into regions of 30h −1 Mpc and 7h −1 Mpc a side centered in a halo of normalised mass M 200b /ρ b ∼ 8.5 × 10 3 hMpc 3 . In the left column we plot the cosmological simulation with the lowest matter density Ωm = 0.148, whereas in the right panel we display that with the highest matter density, Ωm = 1. The solid, dashed and dotted circles show r 200 b , r 200c and r vir radii of the halo.

Examples of haloes matched across simulations
The universality of the mass functions assumes that the mass function is completely described by the linear density field. In order to test this assumption, we have run simulations with very different growth histories that share the same linear density field at z = 0. In Figure 3 we show the simulated density field at z = 0 for our L = 200h −1 Mpc simulations with the most dissimilar growth histories for the ns = 0.9611 cosmology. In the top panel we show a region of 200h −1 Mpc wide, whereas in the middle and bottom panels we zoom on a massive dark matter halo of normalised mass M 200b /ρ b ∼ 8.5 × 10 3 hMpc 3 . We can see that, although both cases corresponding to identical z = 0 linear density peaks, their nonlinear counterparts are different. Specifically, in the case with the highest Ωm value, and thus, highest current growth rate (in the right), haloes are significantly less dense in its center, which is consistent with its expected lower formation redshift and thus lower concentration parameters. Various definitions of halo radii are displayed by white circles in each case. By comparing r200c and rvir radii in the two cosmologies, we see that they identify very different regions of the halo, unlike r 200b which is similar in both cases.
To explore this further, we will compare haloes of the same peak height in different simulations. In particular, we have cross-matched halo catalogues among simulations that share the same linear power spectrum at z = 0. For this, we associate two haloes based on their position and peak height. From Figure 3 we expect that the same fluctuation in the linear density field will end up having a different mass depending on its nonlinear evolution.
We will characterise each halo by an "effective growth rate" and an effective "local power spectrum slope", which we define respectively as: where aev is defined implicitly via D(aev) = γD(a) with γ = 4/5, and where κ = 1, and RL is the Lagrangian radius of a halo of mass M . Physically, these two parameters will be capturing how quickly haloes have recently grown and the density profile of the collapsing region, which can be considered as a proxy for the full mass accretion and merger history of a given halo. Note that the effective growth rate is not evaluated at the redshift in which we identify a halo, but it is evaluated in the past, i.e. γ < 1. By this we seek to capture not the rate of current mass accretion, but instead the amount of mass that has been accreted recently. We have tried different definitions of α eff and found that this distinction was particularly important for models with dynamical dark energy. We chose the numerical value for γ as that which provided the most accurate and simplest model for the halo mass functions, as we will show in Section 6.
In Figure 4 we display the spherically averaged mass distribution around crossmatched haloes in 2 bins of the peak height, ν ∼ 1 and ∼ 2. We display the average density profile, the cumulative mass profile, and the logarithmic slope of the halo density profiles. Vertical lines indicate the radius at which the average enclosed density reaches a value equal to 200 times the background, virial and critical density, as indicated in the legend. Coloured lines indicate three growth histories for the simulations with ns = 0.9611.
We can see that generically the mass profiles differ systematically with α eff , at all values of ν. The higher the growth rate the lower the enclosed mass respect to the background density at a given physical radius. We emphasise that all these objects share the same shape and amplitude of their linear overdensity field at z = 0. Thus all changes necessarily are caused by the different growth history.
The different growth histories are expected to influence the internal structure of haloes. In particular, lower growth rates are expected to cause lower current accretion rates onto haloes, which implies haloes formed earlier and thus are expected to have higher concentrations. In the first row of panels we see that this is indeed the case. Inner regions of haloes appear more concentrated. However, the changes are not limited to the concentration, as external parts are also modified increasing their density the higher the current growth rate. In fact, there seems to be an inflection point located at around r200c radius.
Despite the systematic dependence on α eff , the profiles are very similar when expressed in physical units. However, as a consequence of the pseudo-evolution of the halo boundaries (Diemer et al. 2013), when expressed in r∆ units the profiles become very different. The pseudo-evolution of the boundaries is clear in Figure 4. In the α eff = 1 cosmology r200c is almost three times larger than in the α eff = 0.343 case, while r 200b radii remain roughly constant. Thus, depending on how we define the boundary of our halo, the mass differences will be enhanced or suppressed.
In the lower panels of Figure 3 this can be appreciated visually. 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 r200c, rvir and r 200b radii of the halo. While r 200b defines a halo boundary roughly at the same physical location, r200c compares very different regions of the density field. This effect is less important for rvir, which might explain why Despali et al. (2015) found that Mvir is the mass definition that leads to the most universal behavior.
In order to explore this effect more systematically, in Fig-ure 5 we show how the masses of the crossmatched haloes differ depending on the growth history and the mass definition. We display the ratio of the masses of the crossmatched haloes respect to the α eff = 1 cosmology as a function of ν. We see that in the cosmology with the lowest current growth rate value M200c masses are around 30% smaller than in our reference cosmology. However, for the same cosmologies and haloes, M 200b masses are around 20% more massive. This effect has two contributions. On the one hand, r200c radii lie in the inner parts where the effect of the growth history on the mass profile is larger. On the other hand, because of the pseudo evolution of r200c, we compare the masses enclosed in different physical radii. As a consequence, even if at a given physical radius the enclosed mass is always larger for haloes in low growth rate cosmologies, when comparing M200c masses it seems that haloes in high growth rate cosmologies are more massive. Note that this is solely because we compare masses enclosed in different physical regions. Thus, the non-universality of M200c mass function is in a big part due to the evolution of the boundary itself.
Finally, we want to explore the effect of the local slope of the power spectrum in the crossmatched haloes. For a given power spectrum, redshift and ν, n eff is completely determined. Therefore, in order to see the effect of this variable on the mass of the haloes, we crossmatch the cosmologies with α eff = 0.52 and α eff = 1 for the three power spectra defined with the three ns values we have considered in this work. Note that we only crossmatch cosmologies with the same linear power spectrum. However, if the local power spectrum slope affects the mass of the haloes, we expect the departures of the masses of α eff = 0.52 from α eff = 1 cosmology to be different in the three linear density fields. In Figure 6 we show the results following Figure 5, coloured by the n eff values of the crossmatched haloes. Indeed, we see that for a given ν, the departures of α eff = 0.52 haloes from α eff = 1 haloes are different depending on the n eff value of the halo. Nevertheless, these differences are much smaller than the differences that haloes with different α eff show.
In summary, the whole density profile of the halo is affected by the growth history in a non trivial way. This effect will be reflected in the mass function in a different fashion depending on how masses are defined. Specifically, we expect the nonuniversality of the mass function to change with the mass definition.

The non-universality of the mass function
From the numerical simulations described in the previous section, we have obtained measurements of the halo mass function in 9 cosmologies -3 growth histories and 3 different power spectrum slopes -covering a broad range of masses at various redshifts.
In order to compare these measurements and estimate the impact of the non-universality in the mass function, we have computed νf (ν), where ν is δc/σ(M ). Operationally, we first measure the mass function n(M ) and then estimate νf (ν) inverting Equation 1.
The critical density for collapse, δc(z), can be estimated as 3/5(1.5π) 2/3 Ωm(z) 0.005 (Kitayama & Suto 1996), with an explicit dependence on Ωm. In this work we approximate it with the value that corresponds to a universe where there is only matter, δc = 1.686. We do this to simplify the redshift The bottom panel shows the ratio relative to the mean value in each ν bin. We show in gray the z > 0 mass functions, while we colour the z = 0 mass functions according to their α eff value. The shaded area corresponds to measurements with ν < 0.7, which will be excluded when developing a fitting function for f (ν). and cosmology dependence and make it easier to model. We have checked that the deviations from universal behaviour of the mass function are stronger (up to ±30%) than the effect of taking into account the redshift evolution of the critical overdensity for collapse, which changes the results around ±10% (see Appendix A).
In Figure 7 we display the measurements of νf (ν) from M 200b mass function in the cosmologies listed in the Table 1 at z = 0, z = 0.5 and z = 1. The mass functions with z > 0 are displayed as gray lines, while z = 0 mass functions are coloured according to their α eff value. The average value in each ν bin is displayed as a black solid line, and it is used as a reference for the ratio displayed in the bottom panel. To avoid possible biases due to differential coverage of our models 1 , we will restrict our subsequent analysis to the range 0.7 < ν < 5. For the cosmology most consistent with observational constraint, this implies a mass range of 10 11 < M/[h −1 M ] < 10 15 at z = 0.
In this figure we can clearly see deviations from an universal behavior. For ν values above unity, haloes of a given peak height in one cosmology can be up to 70% more abundant than in others. By construction, the origin of this nonuniversal behavior must be in a combination of the different statistics of the initial Gaussian random fields and the different growth histories. Indeed, we can already see that this is 1 At a fixed volume and number of particles, the differences in the power 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 density field at z = 0. At a given ν, haloes seem to be more abundant the lower the growth rate value. Recall that, as we saw in the previous section, this is a consequence of the same fluctuation being more massive for low growth rate values.
We now explore how these deviations correlate the value of the effective growth rate and power spectrum slope at any redshift. In Figure 8 we display the deviations from the mean νf (ν), at a fixed ν as a function of n eff and α eff (horizontal and vertical axis). Each panel shows the result for a different mass definition.
For all mass definitions we can see that the non-universality clearly correlates with these properties -despite them being only a proxy of very different merger and assembly histories. In the second panel we see that for M 200b , deviations around the mean can reach ±20%. Cosmologies that have higherthan-the-mean νf (ν) typically have lower growth rate values, whereas those with higher growth rate values lead to lower abundances. At fixed α eff , deviations from universality are much smaller, about 15%, and they correlate with n eff . Note that here we are plotting measurements of many redshifts, therefore, we expect that the redshift evolution of the mass function could be described through the dependence on these physically-motivated variables.
However, the non-universality of the mass function depends on the mass definition. Among those considered in this work, M200c (MFoF) mass functions are the most (least) nonuniversal with deviations up to ±30% (±10%) around the mean. It is interesting to note that, as Despali et al. (2015) found, Mvir mass functions are the most universal among the SO mass functions. Furthermore, the dependence on the growth rate is inverted in M200c and Mvir cases with respect to M 200b case (in agreement with Diemer 2020).
In summary, haloes of a given peak height, ν 200b , are more abundant the lower the growth rates and the shallower the power spectrum slope. In other words, a halo that forms early and has grown mostly trough minor mergers, will be more massive than another that has recently formed and has experienced a lot of major mergers, even if both have an identical peak height in the linearly extrapolated initial field.
There could be different paths to follow at this point. One could be to find the halo boundary definition that minimises the non-universality of the mass function. In fact, we have seen that mass definitions based on the critical density induce strong pseudo-evolution in the mass function driven by the change of the boundary of the halo. One quantity that has been argued is more physical is the turnaround radius, which by definition encloses the outermost shell that has collapsed. In the same direction, the first explorations of the splashback mass functions have been done (Diemer 2020). Other alternatives have been recently proposed, which are claimed separate better the linear and nonlinear regimes (e.g. Garcia et al. 2020;Fong & Han 2020).
Any universal mass definition, as we saw earlier, would correspond to very large scales, which although perhaps better suited for describing the mass distribution, might not describe equally well, and thus it might display less correlation with the properties of collapsed gas and of the galaxies hosted by the halo. In addition, many of the proposed halo definitions are ambiguous to implement numerically.
Another option would be to develop a model for the changes of the full density profile as a function of the mass ac- cretion history. For instance, in an analogous manner to the models developed for the relationship between the concentration and the expected mass accretion history in Extended Press Schechter (Ludlow et al. 2016(Ludlow et al. , 2019, it is perhaps possible to develop a model for the outer regions of a halo, which would then predict the changes in halo mass at any radius.
The option we will follow here is to adopt a standard halo definition but calibrate the predictions for the halo abundance to be a function of the peak height but also of the properties of the cosmological model. We will show that with a simple parameterization in terms of α eff and n eff , we can accurately describe the halo mass function for a large region in cosmological parameters space.

MODELLING THE DEPENDENCE ON GROWTH FUNCTION AND POWER SPECTRUM SLOPE
In the previous sections we showed how the halo mass function varies systematically with growth rate and slope of the power spectrum. In this section, we will model this dependence explicitly. We will focus on the M 200b mass function, because as it has been discussed in the previous Section, it is the most physically motivated choice and presents the least pseudoevolution among the overdensity mass definitions. FoF mass function could also be a good candidate, but the fact that it has no clear observational counterpart and that the masses are very subject to numerical effects make it a less interesting candidate. However, in Appendix C we show that our approach is valid to model the mass function of any of the other mass definitions considered in this work.
In each case, we find the best fitting parameters by minimizing the χ 2 of the quantity νf (ν). The uncertainty in each measurement is given by the Poisson statistics and a systematic uncertainty of 5%, as discussed in Section 2. We impose a limit of 200 particles per halo and 400 haloes per mass bin. The minimization is done with the optimize.minimize package of scipy, imposing bounds on the possible values that the parameters may take. Other than our main model (Eq. 5), we have found the best fit parameters for the functional forms that depend only on ν (Eq. 6), ν and n eff (Eq. 6 × 7) and ν and α eff (Eq. 6 × 8). We list all the best-fit parameters in Table 3.
We now asses how well this model is able to describe our calibrating data. In Figure 9 we display the ratio of measured νf (ν), in all our simulations at all three redshifts, to their corresponding predictions of Equation 5 (red). In the left, middle, and right panels we display the residuals as a function of ν, n eff , and α eff , respectively. In all cases, solid and dashed lines display the mean and the median, whereas the shaded areas denote the regions enclosing 90% of the measurements. Thus, this plot quantifies the overall accuracy of each model in 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 of values explored, with no noticeable remaining dependence with either parameter. For comparison, we display also residuals with respect to a version of Equation 6 where we have measured their parameters to our whole dataset but only adopting dependence with respect to ν (blue). As expected, in this case, the residuals are significantly larger, reaching variations of ±20%.
Although the lack of residual dependence with α eff and n eff is achieved by construction, it is in principle not guaranteed that the amplitude of these residuals decrease significantly. For instance, the mass function could have shown dependence on many more details of the assembly history of haloes and the statistics of peaks than simply on the effective growth rate and power spectrum slope. It is, therefore, remarkable that the residuals in the mass function are all contained within a region of ±5%, consistent with our statistical uncertainties for the mass function.  . Deviations between νf (ν) measured in our N -body simulations and the predictions of the fitting functions developed in this work. Left, middle, and right panels display these deviations as a function of the peak height, ν; the effective power spectrum slope, n eff ; and the effective growth rate, α eff , respectively. In each panel, shaded regions indicate the region that contains 90% of our simulated results when employing a fitting function calibrated only as a function of ν (blue) or additionally including dependence with respect to n eff (green), α eff (orange), or both of them (red).
In the next section we will explore whether our approach is actually able to describe accurately the mass function in multiple cosmologies currently allowed by observational data.

VALIDATION: HALO ABUNDANCES AS A FUNCTION OF COSMOLOGY
To asses the accuracy of our description for the halo mass function, we will compare its predictions against a suite of simulations with 30 different cosmologies. Each of our simulations evolved 1536 3 particles inside a box of approximately L = 512h −1 Mpc. The initial conditions where created using 2nd-order Lagrangian Perturbation theory at zstart = 49 and fixing the amplitude of Fourier modes (Angulo & Pontzen 2016). The cosmologies were chosen so that they cover a region of approximately 10σ around Planck's best fit values. Specifically, they cover the following parameters ranges:  . The vertical lines display the limits of the bins with at least 400 haloes resolved with more than 200 particles. Each row, from top to bottom, display variations in ns, Ωm, Ω b and σ 8 for the first column and h, Ων , w 0 and wa for the second column. The gray lines display the result for the fiducial cosmology. The shaded areas denote regions of ±5% and ±2%.
(2013); and dynamical dark energy with an equation of state w(z) = w0 + (1 + z) wa. The cosmology of each simulation is obtained by changing one cosmological parameter of a fiducial cosmology while keeping the rest fixed. The fiducial cosmology assumes flat geometry, massless neutrinos (Mν = 0), a dark energy equation of state with w0 = − 1 and wa = 0, an amplitude of matter fluctuations σ8 = 0.9, cold dark matter density Ω cdm = 0.265, baryon density Ω b = 0.05, and normalised Hubble constant h = 0.6.
In Figure 10 and 11 the solid lines represent the measured mass functions of the 30 simulations relative to the predictions of the main model developed in this work (Eq. 5). We show the limits of the bins with at least 400 haloes with resolved with more than 200 particles as vertical lines.
In each row we show the mass functions of the cosmologies where we vary one cosmological parameter keeping the rest fixed. In Figure 10 we display the results at z = 0, and in Figure 11 the results at z = 1. For , we also show the residuals respect the model developed in Despali et al. (2015) as dotted lines, which assumes universality of the mass function. We recall that the functional dependence on f (ν) is the same in both models, while in our model we have added extra dependences on n eff and α eff in order to capture the effect of growth history on the mass function.
At z = 0, our model describes the low mass end of the mass function at an accuracy of 3%, while Despali et al. (2015) predicts that haloes are 10% more abundant. This may be a consequence of using different group finders, SO and FoF respectively. For haloes with masses above M > 10 14 h −1 M , there seems to be an underprediction of our fitting function. To investigate this, we have compared our predictions against the simulations of Angulo et al. (2020), which feature the same mass resolution as our test suite but on a volume 27 times larger. Although not shown here, in such case we find an agreement to better than 5% up to 10 15 h −1 M . Combined with the good agreement of our predictions with those of Despali et al. (2015), we speculate that there is a systematic over prediction of the abundance of haloes in our test sims for M > 10 14 h −1 M , which could be caused by finite-volume effects.
At z = 1, the redshift evolution of the mass function is evident. Even if at z = 0 Despali et al. (2015) is a good description to the mass function, at z = 1 it overpredicts the abundances for more than 10%. Our model captures this and yields results that are accurate within 5% at all masses considered.
Compared to the redshift evolution, the cosmology dependence of the mass function seems to be weak. However, the mass functions of w0 cosmologies present strong deviations from universality. The scatter of the ratio with respect to Despali et al. (2015) is of ∼ 5% among cosmologies with different w0 values at both redshifts. After taking into account the dependences on α eff and n eff , this scatter vanishes. We emphasise that we have only used ΛCDM cosmolgies to calibrate the fit, and so α eff and n eff are physically meaningful proxies of the non-universality of the mass function.
It is interesting to notice that the largest part of the improvement is obtained when adding α eff to the universal description. This is expected, because as discussed in Section 3, the deviations from universality correlate much stronger with α eff than with n eff .

APPLICATION: IMPROVING THE ACCURACY OF COSMOLOGY-RESCALING METHODS
To approach an optimal exploitation of current and future observations of the abundance of dark matter haloes and the clustering of galaxies, very accurate theoretical predictions for these quantities are required. Although fitting functions and calibrated recipes are indeed extremely valuable, they fall short in providing correlations among different observables or the full three dimensional distribution of clusters of galaxies. One option to obtain those predictions is to employ cosmological N -body simulations together with cosmology rescaling algorithms.
The basic idea of cosmology rescaling is to employ a few simulations carried out adopting specific cosmological parameters, and then manipulate their outputs so that represent nonlinear structure in any other set of cosmologies. These algorithms have been extensively discussed and tested in Angulo & White (2010); Angulo & Hilbert (2015); Ruiz et al. (2011);Mead & Peacock (2014b,a); Renneby et al. (2018). In particular, Zennaro et al. (2019) showed these are applicable to cases of massive neutrinos, and Contreras et al. (2020) showed that the clustering of dark matter and dark matter haloes and subhaloes can be obtained to better than 3% accuracy from large to very small scales (0.01 < k/h Mpc −1 < 5). This technique has been recently employed by Angulo et al. (2020) to predict the nonlinear power spectrum as a function of cosmology, by Aricò et al. (2019Aricò et al. ( , 2020 to model the effect of baryonic physics and predict the suppression of the power spectrum due to baryons, and by Zennaro et al. (2021) to model the clustering of biased tracers.
In Figure 12 we compare the performance of cosmology rescaling in predicting the halo mass functions. Specifically, we compare measurements in the 30 simulations described in the previous section to the halo mass function after rescaling one simulation. We refer the reader to Contreras et al. (2020) for details on how the cosmologies of these simulations were chosen.
In the original cosmology-rescaling, the simulation volume and the particle mass are rescaled by a single factor, found .The vertical lines display the limits of the bins with at least 400 haloes resolved with more than 200 particles. Each row, from top to bottom, display variations in ns, Ωm, Ω b and σ 8 for the first column and h, Ων , w 0 and wa for the second column. In each panel we show results before (dotted lines) and after (solid lines) applying our additional correction accounting for dependence on growth history, as indicated by the legend (see text for details). The shaded regions denote regions of ±5% and ±2%. by minimizing the difference in the linear mass variance in the target and rescaled cosmologies. In this operation the number of particles in each halo is left invariant. Using this recipe the scaling algorithm sets the same linear density field in the rescaled and target simulations, which is equivalent to assuming the universality of the mass function. However, in this work we have shown that the mass function depends not only on the linear density field but also on the entire growth history. In fact, we can see this effect in the dotted lines of Figure 12, where we compare the mass functions from the rescaled simulations with the target ones. By assuming universality of the mass function, the rescaled M 200b mass functions differ from the target mass functions up to 10% in some cosmologies at z = 0.
Our model for the dependence of the halo mass function on growth history gives us the possibility to construct an additional correction for cosmology rescaling by taking into account the different growth histories the target and rescaled cosmologies have gone through. Specifically, we first predict the rescaled and target mass functions. In the target cosmology case, the prediction is straightforward. In the rescaled original cosmology case, we compute the expected mass function of the original cosmology once we have applied the corresponding mass and length scalings, i.e, once we have set the linear density field equal to the target cosmology's linear density field. As we have seen in Section 3.1, the growth-history affects the masses of the haloes, not the abundance of them. Thus, the difference between target and rescaled mass functions is given by a change in mass, which we find by mapping the rescaled masses to the target masses where the abundances are the same. Therefore, for a given pair of originaltarget cosmologies, we can predict a halo-by-halo mass correction that describes the effect of the non-universality of the halo mass function.
We show the results after applying this correction as solid lines in Figure 12. We can see that in all cosmologies, the accuracy of the predictions improve in a clear manner. At z = 0 the differences are in most of the cases smaller than 2%. Note that our model, calibrated on simulations where we only vary Ωm and ns, is able to capture the non-universality of general cosmologies, even in beyond ΛCDM cosmologies with 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 ±1 − 2% over the full range of masses where we can measure the mass function accurately.

SUMMARY AND CONCLUSIONS
In this paper we have studied the non-universality of the halo mass function. We have run simulations with very extreme cosmologies to maximise the deviations from universal behaviour and we have shown that the halo masses are affected by the entire growth history. As a consequence, given the same linear density field in two different cosmologies, the halo mass functions are different.
In order to shed light in the origin of the non-universality of the mass functions, we have crossmatched haloes of different cosmologies that share the same linear density field and we have compared their density profiles. Generally, we have observed that all the density profiles up to very large radii are affected by the growth history of the haloes. Furthermore, the physical boundaries of haloes selected with density criteria are subject to pseudo evolution, and correspond to different physical radii for different cosmologies. This effect is more pronounced for overdensities defined with respect to the critical density of the universe. Therefore, different mass definitions yield mass functions with different dependences on redshift and cosmology.
We have modelled the non-universality of the mass function adding two additional parameters other than the peakheight ν: the effective growth rate, α eff , and the local slope of the power spectrum, n eff . Using a total of 8 free parameters, our model captures the non-universality and can lower the scatter on the halo mass functions in all the cosmologies considered from ±20% to about ±5% up to z = 1. In the literature, the redshift evolution of the mass function is typically parametrised explicitly within a fiducial cosmology (see e.g. Tinker et al. 2008). On the contrary, here we have modelled simultaneously the cosmology and redshift non-universality of the mass function by using physically motivated parameters.
We have tested our model on an independent set of simulations of 30 different cosmologies, including massive neutrinos and dynamical dark energy. By considering the α eff and n eff dependences, we have been able to reproduce the halo mass functions within a 5% accuracy in all the cosmologies up to M ∼ 5 × 10 14 h −1 M until z = 1. We emphasise that the simulations that we have used to calibrate the models have been run with ΛCDM cosmologies. Thus, it is not a trivial result that our model is able to describe the halo mass functions within, for instance, cosmologies that include massive neutrinos or dynamical dark energy.
As an application of our model, we have applied it together with the cosmology rescaling method presented in Angulo & White (2010). We have found that the accuracy in the scaling of the halo mass function improves from 10% to 2% in all 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 future works. It is well known that baryonic processes alter in non-trivial way the halo mass function. In particular, astrophysical feedback ejects a large amount of gas outside the haloes boundaries, and therefore haloes result less massive, even by when including a baryonic modelling (see e.g. Castro et al. 2021;Debackere et al. 2020). We plan to extend our formalism to include the effect of baryons on the halo mass function, by using the so-called baryonification technique (Schneider & Teyssier 2015; Aricò et al. 2019). By combining it with cosmological-rescaling algorithms, we will construct an emulator of the halo mass functions, as a function of cosmological and astrophysical parameters.
In the near future, more precise and accurate predictions of the halo mass function will be necessary in order to fully exploit the data of the future surveys. As an example, Artis et al. (2021) estimated that, only considering the precision of the parameters of the fitting functions in the analysis (i.e assuming universality of the mass function), an improvement from 30% to 70% is required. This framework provides us with a very accurate fit of the halo mass function, which can be eventually exploited to directly compare against observed clusters count, from optical, X-ray or Sunyaev-Zel'dovich surveys. Thus, we anticipate that this model will be of great value value in constraining the cosmological parameters of the Universe.

DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author. Figure A1. The difference of the deviations of the mass function respect to the mean computed in each ν bin between the mass functions obtained with the critical density for collapse presented in (Kitayama & Suto 1996) and with the critical density corresponding to a universe with only matter. ∆ Kitayama,EdS = (f (ν) Kitayama  In this appendix we show the effect of taking into account the redshift dependence of the critical density for collapse on the non-universality of the mass function. Specifically, we compute the relative difference of the deviations in each ν bin between the mass functions with δc(z) and δc = 1.686. In Figures A1 and A2 we show these relative differences for two different approaches of computing δc(z). We see that, for most of the cases they do not exceed the ±10%, while the deviations around the mean are much stronger, as we can see in Figure 8.

APPENDIX B: REDSHIFT CORRECTION
For different box sizes the redshifts of the snapshots vary slightly. This effect is more pronounced at high redshift, where the difference of the output redshifts of different boxes can reach ∆z ∼ 0.01. In this time lapse the mass functions may have evolved, therefore, when combining different box sizes we may be introducing some bias in our data set. The left panel of Figure B1 displays the expected ratio of the differential mass functions between the output redshifts of the different boxes around z = 1. At M ∼ 10 15 h −1 M , the differences can reach 10%.
In order to test whether the predicted evolution of the mass function is accurate, we make use of a simulation presented in Section 6, for which we have many snapshots. In the right panel of Figure B1 we display the predicted and measured ratios for the redshifts listed in the legend. We can see that the predictions are a good description of the data. For other cosmologies the results are similar.
Thus, we proceed to correct the mass functions of the big box sizes in the following way, where f (M, z) is some model for the differential mass function.

APPENDIX C: EXTENSION TO OTHER MASS DEFINITIONS
In this appendix we present the main results of our modeling with other mass definitions. These results are analogous to what already described for M 200b mass functions.
In general, as seen in Figure 8 all mass functions show clear 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 change the functional form of the peak-height dependence for the different mass definition. For M∆ mass functions we use the where {A, b, c, d} are the free parameters of the ν dependence of our model. By applying the methodology explained in Section 4, we have obtained the best-fit parameters listed in Tables C1, C2 and C3 for MFoF, M200c and Mvir mass functions respectively. In Figure C1 we show the performance of the model for each mass definition. For MFoF mass functions the residual scatter is consistent with the intrinsic uncertainties of our fit. However, even if the scatter is reduced significantly, for Mvir and M200c the description is not as good as in the other cases. We think that, for both the cases, this may be a consequence of the pseudo-evolution of the boundary definition. This paper has been typeset from a T E X/L A T E X file prepared by the author.