© The Authors 2026

1. Introduction

The use of the halo mass function (HMF) to describe the abundance of dark matter halos as a function of mass is a cornerstone of modern cosmology. It encodes the growth of structure from the early Universe to the present day. It provides a fundamental link between theory, simulations, and observations, making it a key probe of cosmological models and parameters (Cooray & Sheth 2002; Mo et al. 2010; Castro et al. 2016). Accurate predictions of its mass and redshift dependence are therefore essential for interpreting data from current and upcoming large-scale structure surveys.

The theoretical foundation of the HMF was laid by the Press–Schechter (PS) formalism (Press & Schechter 1974), later developed through excursion-set theory (Bond et al. 1991; Lacey & Cole 1993) and refined with ellipsoidal-collapse models (Sheth & Tormen 1999; Sheth et al. 2001, hereafter ST). Empirical calibrations based on N-body simulations have since improved its accuracy over a range of cosmologies, halo masses, and redshifts (e.g. Jenkins et al. 2001; Tinker et al. 2008; Despali et al. 2016; Shirasaki et al. 2021; Seppi et al. 2021). However, all existing fits remain limited in scope: early models can be applied mainly to cluster-scale halos at z ≲ 3, while later extensions reach z ≲ 10, but typically only for halo masses above 10¹⁰ h⁻¹ M_⊙. As a result, no current model provides a self-consistent description of the HMF across the full halo mass hierarchy and cosmic history.

Recent advances in numerical simulations, particularly the UCHUUN-body simulation suite (Ishiyama et al. 2021; Ishiyama et al., in prep.), make it possible to overcome many of these limitations by combining extremely large cosmological volumes with a high mass resolution. This enables precise measurements of the HMF over a vastly extended range of halo masses and redshifts, offering a unique opportunity to revisit and improve theoretical models.

In this work, we build on the generalised Press–Schechter formalism of Betancort-Rijo & Montero-Dorta (2006a,b) and present a new model, named GPS+, which delivers a consistent description of the HMF from dwarf-scale halos to the most massive clusters, spanning 6.5 < log(M/M_⊙) < 16 and redshifts 0 ≤ z ≤ 20. While it remains rooted in the generalised PS framework, triaxial collapse plays a central role in the construction of the model. With minimal calibration to the UCHUU simulations, the resulting HMF depends primarily on the variance of the linear density field, σ(M, z), exhibiting only weak explicit mass dependence and with no explicit dependence on redshift. This unprecedented dynamical range and per cent-level accuracy are enabled by the scope and quality of the UCHUU suite, providing a robust framework for interpreting data from ongoing and future surveys such as DESI, Euclid, Subaru PFS, LSST, eROSITA, and JWST.

While our approach is still closely connected to the classical PS picture, it replaces the original collapse assignment with a more physically motivated prescription, while preserving the statistical properties of Gaussian random fields and the conceptual structure of Press–Schechter theory.

This paper is organised as follows. Section 2 describes the UCHUU simulations and HMF measurements. The GPS+ theoretical framework is presented in Section 3. In Section 4, we compare the model predictions with simulation results and Section 5 summarises our conclusions.

2. Uchuu simulations

In this work, we tested the theoretical framework using the UCHUU simulation suite spanning a wide range of mass resolutions, including PHI-4096, SHIN-UCHUU, UCHUU, MUCHO-UCHUU-140M, MUCHO-UCHUU-1G, and MUCHO-UCHUU-6G¹, as described in Ishiyama et al. (2021); Ishiyama et al. (in prep.). Their main properties are summarised in Table A.1. The first three simulations are based on Planck15 cosmology, while the remaining three follow Planck18. Despite these slight cosmological differences, their impact on the HMF is negligible compared to both the intrinsic statistical uncertainties of the simulations and the theoretical modelling uncertainties.

We used 50, 300, and 100 independent realisations of MUCHO-UCHUU-140M, MUCHO-UCHUU-1G, and MUCHO-Uchuu-6G, respectively. The HMF is obtained by averaging over all realisations, with uncertainties estimated from the diagonal of the covariance matrix, an ensemble approach that significantly reduces cosmic variance. For the purposes of illustration, within the total volume of 21 600 h⁻³ Gpc³ covered by the 100 MUCHO-UCHUU-6G realisations, the most massive halo is characterised by M = 8.3 × 10¹⁵ M_⊙, about 13 times the mass of the Virgo cluster (Kashibadze et al. 2020).

All simulations were run with the TreePM code GreeM (Ishiyama et al. 2009, 2012). PHI-4096, SHIN-UCHUU, and UCHUU were executed on the ATERUI II supercomputer at CfCA, National Astronomical Observatory of Japan, while the MUCHO-UCHUU simulations used the FUGAKU supercomputer at RIKEN, Kobe. Halos were identified with the RockStar halo/subhalo finder (Behroozi et al. 2013a) and its distributed version MPI-Rockstar² (Tokuue et al. 2024), and merger trees were constructed using Consistent Trees (Behroozi et al. 2013b).

Figure 1 shows the combined HMF derived from all the simulations used in this work. We express it in terms of the the multiplicity function (e.g. Sheth & Tormen 1999), defined as

Fig. 1. Multiplicity function f(M, z) (Eq. (1)) from UCHUU simulations (markers) vs. GPS+ predictions (solid lines) at z = 0–20 (right to left; see Table B.2). The y-axis shows the normalised differential halo abundance including dlnM/dlnσ. Only bins with > 200 halos are shown.

$$\begin{array}{c}\hfill f(M,z)=\frac{M_{200\mathrm{m}}}{{\rho }_{\mathrm{m}}}\frac{\mathit{dn}}{dln{M}_{200\mathrm{m}}}\left|\frac{dln{M}_{200\mathrm{m}}}{dln\sigma }\right|.\end{array}$$(1)

We adopted the halo mass M_200m, defined within an overdensity of 200ρ_b, where ρ_b is the mean background density. Unlike the virial mass, this fixed-threshold definition is redshift-independent. We refer to Appendix E for more details.

To construct the combined HMF at each redshift, we used halos with at least 407 particles to ensure low-mass completeness. The high-mass end of each simulation is limited to allow for a smooth, continuous transition across overlapping mass ranges. The resulting HMFs are available at the CDS (see Appendix B).

3. Theoretical framework

We modelled the HMF using the GPS+ framework, which describes halo abundances over an exceptionally wide mass range (6.5 ≲ log[M/(h⁻¹ M_⊙)] ≲ 16) and up to very high redshifts (z ∼ 20). The formalism builds on the PS idea that a mass element belongs to a halo of mass larger than M if the smoothed density field satisfies a collapse condition; however, it replaces the original centre-based criterion with a more physically motivated prescription. The model preserves the Gaussian-field statistics underlying PS and excursion-set theory while incorporating the physics of triaxial collapse. In this way, GPS+ retains the conceptual structure of PS but modifies the criterion for assigning mass elements to collapsed objects of mass larger than M.

In this framework, the HMF is given by Montero-Dorta & Betancort-Rijo (2006):

$$\begin{array}{c}\hfill \frac{dn(M,z)}{\mathit{dM}}=\frac{{\rho }_m}{M}\frac{dF(M,z)}{\mathit{dM}},\end{array}$$(2)

where ρ_m is the matter density, and F(M, z) is the mass fraction.

The mass fraction is defined as

$$\begin{array}{c}\hfill F(M)=\frac{\mathrm{erfc}[⟨{\delta }_c(\sigma (M,z),M)⟩/(\sqrt{2}\sigma (M,z))]}{V(\mathrm{\Sigma }(M,z),M)},\end{array}$$(3)

with V(Σ(M, z),M) given by

$$\begin{array}{cc}\hfill V(\mathrm{\Sigma }(M,z),M)& =3{\displaystyle {\int }_0^1}\mathrm{erfc}[\frac{⟨{\delta }_c⟩(\sigma (M,z),M)}{\sqrt{2}\mathrm{\Sigma }(M,z)}\hfill \\ \hfill & \times \sqrt{\frac{1-exp(-c(M){\xi }^2)}{1+exp(-c(M){\xi }^2)}}]{\xi }^{2}d\xi ,\hfill \end{array}$$(4)

where δ_l encodes the linear growth corrections at z = 0, Σ(M, z) is given by

$$\begin{array}{c}\hfill \mathrm{\Sigma }(M,z)=\sqrt{{\sigma }^2(M,z)+U^2(\sigma /{\delta }_c)},\end{array}$$(5)

where σ(M, z) is the standard linear-theory variance, and the correction function U(x) is defined as

$$\begin{array}{c}\hfill U(x)=-0.01507+0.17810·x+0.03835·x^2-0.00221·x^3,\end{array}$$(6)

where x = σ(M, z)/δ_c, and δ_c is the critical overdensity, equal to the linear-theory density contrast, δ_l, in the spherical-collapse approximation at collapse. In our model, the critical overdensity for collapse is further modified as a function of mass and redshift, expressed as

$$\begin{array}{cc}\hfill ⟨{\delta }_c⟩& ={\delta }_c{(1+0.845x-0.04x^2+0.0025x^3)}^B\hfill \\ \hfill & \times A{(1+0.17b(M)-0.087b^2(M))}^D\hfill \end{array}$$(7)

where b(M) is defined by

$$\begin{array}{cc}\hfill logb(x)& =-1.28+0.05781x-0.005622x^2\hfill \\ \hfill & +5.884\times {10}^{-4}{x}^3-1.365\times {10}^{-5}{x}^4,\hfill \end{array}$$(8)

where x = log(M/[M_⊙/h]), while A and B are free parameters fitted to the simulation data (see values below) and D is a parameter whose theoretical value is 1 (see Betancort-Rijo et al., in prep.). However, we did allow it to to test the robustness of this prediction. Equation (7) essentially comes from the physics of the triaxial collapse (see Appendix C for the details).

The function c(M) from Equation (4) is a polynomial defined as

$$\begin{array}{cc}\hfill logc(x)& =-1.124+0.01756·x+0.002539·x^2\hfill \\ \hfill & -6.438\times {10}^{-5}·x^3+4.726\times {10}^{-6}·x^4,\hfill \end{array}$$(9)

where x = log(M/[M_⊙/h]). The quantities entering Eqs. (6)–(7) follow from the physics of triaxial collapse: U² can be obtained from first-principles calculations in Montero-Dorta & Betancort-Rijo (2006) and Eq. (6) is a fit to those results. The functions b(M) and c(M) are fully determined by the linear matter power spectrum (Sánchez-Conde et al. 2007; Betancort-Rijo & Montero-Dorta 2006b). The polynomials given by Eqs. (8) and (9) correspond to fits done in this work for the Planck power spectrum.

This formulation enables the computation of the HMF across a broad halo mass range and out to very high redshifts. A concise overview of the theoretical framework is provided in Appendix C, while a full derivation of the equations will be presented in Betancort-Rijo et al. (in prep.).

The coefficients A and B in Eq. (7) were obtained by fitting Eq. (2) to the simulation data using the statistic,

$$\begin{array}{c}\hfill {\mho }^2={\displaystyle \sum _{i,j}}\frac{|n_{\mathrm{sim}}(M_i,z_j)-n_{\mathrm{th}}(M_i,z_j)|}{n_{\mathrm{sim}}{(M_i,z_j)}^2}[1-\mathrm{\Theta }(\frac{\mathrm{\Delta }{n}_{\mathrm{sim}}(M_i,z_j)}{0.05n_{\mathrm{sim}}(M_i,z_j)}-1)],\end{array}$$(10)

where n_sim and n_th are the simulated and theoretical HMFs, and Δn_sim is the statistical uncertainty. The Heaviside function Θ restricts the fit to bins with relative uncertainties below 5%.

Uncertainties were assumed to be Poissonian for simulations with a single realisation (PHI-4096, SHIN-UCHUU, UCHUU), while for MUCHO-UCHUU-140M, MUCHO-UCHUU-1G, and MUCHO-UCHUU-6G, they were derived from the diagonal of the covariance matrix over multiple realisations.

This fit gives A = 1.089 and B = 0.652. Varying D indicates the best agreement at D ≃ 1, consistent with theory. We note that only A and B were fitted to the simulations and their values are close to the expected theoretical value of 1.

A key feature of this formalism is that f(σ) depends explicitly on mass through the functions b(m) and c(m) (Eqs. (8) and (9)), with no explicit dependence on redshift. If this mass dependence is (incorrectly) neglected and we assume that f depends only on σ, we then obtain f(m, z) = f(σ(m, z),m) = f(σ, m(σ, z)) = f(σ, z). The first equality is correct: at fixed mass, any redshift dependence enters only through σ(m,z). The second equality instead removes the explicit mass dependence of f at fixed σ (a central aspect of our framework) by using the m–σ relation at fixed redshift, thereby introducing spurious σ and redshift dependences.

In our formalism, f depends explicitly on both m and σ; the artificial redshift dependence appears only when the mass dependence is forced to act exclusively through σ. This leads to the ∼ 50% variation shown in Fig. D.1, about four times larger than the typical modelling uncertainty of the (slightly simplified) framework used here. Since this behaviour follows directly from the explicit mass dependence predicted by the theory, we interpret the trend in Fig. D.1 as lending strong support to the theoretical basis of our formalism.

4. Performance of the theoretical framework and comparison with the Sheth and Tormen HMF

In this section, we present the HMF predicted by the theoretical framework developed and calibrated in this work, comparing it with the UCHUU simulation results described above; for reference, we also include the widely used Sheth–Tormen (ST) model³, while further comparisons will be presented in Ishiyama et al. (in prep.). Figure 2 shows the ratio between the simulated HMF and the GPS+ prediction across a broad redshift range. At z = 0, the ratio remains close to unity over most of the mass range, with deviations reaching ∼20% only at the high-mass end (log₁₀(M_200m/[h⁻¹ M_⊙]) > 14), while at z = 4.27 similar deviations appear at the low-mass end (log₁₀(M_200m/[h⁻¹ M_⊙]) < 11). At higher redshifts, the agreement improves further, with deviations typically below 10% across nearly the full mass range, demonstrating the robustness and predictive power of the GPS+ model over an exceptionally wide range in mass and redshift.

Fig. 2. Ratios of the HMF predictions from GPS+ and ST at different redshifts. Only mass bins with > 100 halos are shown.

The ST model, which extends PS by incorporating ellipsoidal collapse and introduces an explicit redshift dependence through the linear growth factor (noting that ST adopts a friends-of-friends mass definition, whereas we use M_200m), performs similarly to GPS+ at z = 0. At higher redshifts (z ≳ 4.27), however, GPS+ provides a markedly improved description, particularly for the most massive halos. This discrepancy increases toward earlier times: at z ≃ 18.37, ST deviates by ∼60–80%, while GPS+ remains within ∼10–20%. These results highlight both the accuracy of our framework across the full mass and redshift domains explored, along with the limitations of existing models when they are extrapolated to the high-redshift regime.

5. Summary

In this work, we present a physically motivated theoretical model for the halo mass function based on a generalised Press–Schechter framework with triaxial collapse (GPS+), which accurately reproduces the HMF measured in the UCHUU cosmological N-body simulation suite, spanning a wide range of volumes and mass resolutions within the Planck cosmology. The model is rooted in the PS formalism, but it includes a modified collapse criterion for assigning mass elements to halos above a given mass, as detailed in the appendices. Unlike many previous approaches, GPS+ introduces no explicit redshift dependence, with evolution entering solely through the variance of the linear density field, σ(M, z), enabling predictions up to z ∼ 20 over a halo mass range 6.5 ≲ log(M_h)≲16. Across this full mass and redshift domain, the model reproduces the simulated HMF with typical deviations below 10–20%, with the largest discrepancies confined to the high-mass end. Compared to the Sheth–Tormen model, GPS+ shows a similar accuracy at z = 0.0. However, it provides a substantially improved description at higher redshifts and for the most massive halos, where at z = 18.37 Sheth–Tormen deviates by 70–80%, while GPS+ remains within 20%. The uniqueness of this work lies both in the measurement of the HMF over an exceptionally wide mass and redshift range using the UCHUU simulations, paired with the ability of the GPS+ framework to predict the HMF across this range more reliably than any existing model.

We emphasise that our formalism is fully derived from first principles. Only A and B have been allowed to differ from their theoretical values, reflecting the details of our stabilisation model; thus, it is unsurprising that the best fit yields slightly different values.

Data availability

The GPS+ model is publicly available on the UCHUU GitHub⁴ and will be included in COLOSSUS⁵ (Diemer 2018). The combined HMF from the UCHUU suite is provided in electronic tables (Appendix B). The MUCHO-UCHUU-140M, 1G, and 6G simulation products will be released on the SKIES & UNIVERSES website.

Table B.1 is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/707/L4

Acknowledgments

This work was supported by project AST22_00001_27, funded by the European Union – NextGenerationEU, the Spanish Ministry of Science, Innovation and Universities, the Recovery, Transformation and Resilience Plan, the Andalusian Regional Government (Junta de Andalucía), and the Spanish National Research Council (CSIC). EFG acknowledges financial support from the Severo Ochoa programme (CEX2021-001131-S, MCIN/AEI/10.13039/501100011033). EFG and FP are supported by the Spanish MICINN grant PGC2018-101931-B-I00. EFG and FP thank the Instituto de Astrofísica de Andalucía (IAA-CSIC), the Centro de Supercomputación de Galicia (CESGA), and RedIRIS for hosting the Uchuu DR1, DR2, and DR3 data releases on the SKIES & UNIVERSES platform. TI acknowledges support from the IAAR Research Support Program at Chiba University, MEXT/JSPS KAKENHI (JP19KK0344, JP25H00662), the MEXT program “Promoting Researches on the Supercomputer Fugaku” (JPMXP1020230406), and JICFuS. The Uchuu simulations were performed on the Aterui II supercomputer at CfCA (NAOJ) and on the K computer at RIKEN, while the Mucho-Uchuu simulations were carried out on Fugaku at RIKEN (Project IDs: hp240184, hp250149). The Uchuu Data Releases made use of the IAA-CSIC computing facilities in Spain (MICINN EU-FEDER grant EQC2018-004366-P).

References

Appendix A: Simulation Suite Summary

Table A.1 provides a summary of the cosmological N-body simulations used in this work. The simulation suite spans a wide range of volumes and mass resolutions, allowing us to probe dark matter halo statistics across several decades in mass and redshift while controlling for numerical and cosmological systematics.

The Phi-4096, Shin-Uchuu, and Uchuu simulations form a hierarchical set with increasing box size and decreasing mass resolution, all assuming a Planck15 cosmology. These simulations are particularly well suited for studying halo formation and clustering from galactic to large-scale structure scales within a consistent cosmological framework.

Complementing these runs, the Mucho-Uchuu simulations consist of multiple realizations at fixed volumes, adopting a Planck18 cosmology. The Mucho-Uchuu-140M, Mucho-Uchuu-1G, and Mucho-Uchuu-6G sets are designed to improve statistical precision through ensembles of realizations, enabling robust estimates of cosmic variance and covariance in halo statistics. Their differing box sizes and particle masses provide additional leverage to test resolution effects and volume dependence.

For each simulation, Table A.1 lists the total number of particles, box size, dark matter particle mass, Plummer-equivalent gravitational softening length, number of realizations, adopted cosmological parameters, and relevant references. This diverse simulation suite allows for a comprehensive and systematic comparison between numerical results and theoretical predictions across a broad range of scales.

Appendix B: Halo mass function from simulations

Table B.1 provides an example of the release HMF data, corresponding to redshift z = 0. This table, along with 19 additional ones (20 in total, with redshifts listed in Table A.2), is available in electronic form at the CDS. Each table includes the following columns: M200m, f(M), Δf(M), #halos, Simulation, where M200m is the mass enclosed within overdensity 200ρ_b; f(M) is the multiplicity function (Equation 1); Δf(m) is its associated uncertainty; the number of halos corresponds to the count in each mass bin, #; and the simulation name indicates which simulation contributed the HMF measurement for that bin.

Appendix C: Basic considerations about the theoretical framework

The PS formalism yields a halo mass fraction F(M) that is not correctly normalised ($F(M=0)=\frac{1}{2}$, instead of one). Furthermore, it predicts a high-mass tail that is too steep. This is rather puzzling because, while in other mass ranges there is ample room for improving the theoretical treatment, in this limit the spherical collapse approximation used in the PS formalism should become asymptotically exact. This is due to a selection effect: only proto-halos undergoing an almost perfectly spherical collapse have any chance of having formed by the time being considered. ST combined ellipsoidal collapse with the PS approach and obtained a substantial improvement in the low-mass regime, but they did not solve either the high-mass problem or the normalisation issue.

In Betancort-Rijo & Montero-Dorta (2006a), it was pointed out that the origin of these two problems lies in the PS assumptions themselves. The criterion used by those authors to decide whether a mass element belongs to an object with mass larger than M was that a sphere with Lagrangian radius r(M), centred on that element, satisfies the condition for spherical collapse (i.e. the linear fractional mass density fluctuation within the sphere, δ_l, obeys ≥δ_c). This criterion cannot be exact, since not all mass elements can be located at the centres of the objects that contain them.

We adopted a more realistic and less restrictive criterion: a mass element belongs to a halo of mass M if it lies within a distance r(M) of a point that satisfies the collapse condition described above. This formulation allows an analytical treatment, leading to the mass fraction F(m) expression in Equation 3, where V(Σ, m)/2 represents the ratio between the probability that δ_l ≥ 1.686 within a sphere with radius r(M) centered on a randomly selected mass element and the greater probability that such an element belongs to an object with mass larger than M. In that work, we used the spherical collapse approximation, so that < δ_c > (σ, M) = δ_c is constant (taken as either 1.686, or, 1.676). Under this formulation, both the high-mass limit and the normalisation behave correctly.

In Montero-Dorta & Betancort-Rijo (2006), we used the Complete Zel’dovich Approximation (CZA) (Betancort-Rijo & López-Corredoira 2000), which effectively corresponds to ellipsoidal collapse, to compute the mean value of δ_l within the region (generally ellipsoidal) that collapses into an object of mass M. We denoted this mean value by < δ_l>, and its variance by U². Both quantities were found to depend on M and z only through σ(M, z), to a very good approximation). In the present work, the variance obtained in that study has been fitted using the expression given in Equation 6. For < δ_l> , the result obtained in that work was

$$\begin{array}{c}\hfill <{\delta }_l>={\delta }_c[0.814+0.688\left(\frac{\sigma }{{\delta }_c}\right)-0.033{\left(\frac{\sigma }{{\delta }_c}\right)}^2+0.02{\left(\frac{\sigma }{{\delta }_c}\right)}^3].\end{array}$$(C.1)

This expression corresponds to the first factor in Equation 7. The second factor introduces a small explicit dependence on M whose origin and derivation will be presented in Betancort-Rijo et al. in prep..

Equation C.1 was obtained using a model for shell-crossing and stabilisation in which the protohalo, upon collapsing along one direction in comoving coordinates to $\frac{1}{13}$ of its original size, becomes stabilised along that axis while continuing to collapse along the remaining directions. This stabilisation model is physically well motivated but not exact, so we introduced the parameters A and B, which are expected to remain close to unity, to account for its approximations.

It is important to note that in this formalism F(M) depends on mass not only through σ(M, z), but also through an explicit dependence on M that reflects the shape of the power spectrum at that scale M. This dependence appears in Equation 7 through b(M) (8) and in V(M) (Equation 4) through c(M) (Equation 9).

It should also be emphasised that there is no explicit dependence on redshift, and that F(M) is correctly normalised without requiring parameters fitting. The only fitted parameters are A and B in < δ_c > (σ, M), which, based on theoretical considerations, are expected to remain close to unity. In principle, their values could be determined directly from simulations by analysing the behaviour of this function, since analytically we cannot extend the derivation beyond C.1.

Appendix D: Multiplicity function as a function of σ(M)

In this appendix, we show the multiplicity function, f(σ, z), as a function of σ at different redshifts, as predicted by the GPS+ model presented in this work (Figure D.1).

Fig. D.1. Top panel: Multiplicity function f(σ, z) as a function of σ at different redshifts, as predicted by the GPS+ model presented in this work. Bottom panel: Ratio of the model at each redshift to z = 0 prediction.

Appendix E: Dependence of the HMF with the halo mass definition

As part of this work, we assessed the impact of different halo mass definitions in our theoretical framework. The HMF results presented in Section 4 are based on the M_200m mass definition; however, the virial mass, M_vir, can also be adopted. Figure E.1 compares the performance of the model under both definitions. The largest differences arise at low redshift. When using the virial mass, the agreement at the high-mass end worsens significantly, with ratios dropping to ∼0.5, corresponding to discrepancies of ∼50% between the model predictions and the measured HMFs from the simulations. At higher redshifts, the performance of the model using the virial mass becomes comparable to that obtained using the M_200m definition.

Fig. E.1. Comparison of the ratio between the simulated HMF and the GPS+ model at different redshifts, for two halo mass definitions in the data: the mass enclosed within overdensity 200ρb, M200m, and the mass within the virial radius, Mvir.

This analysis shows that the choice of halo mass definition has a strong impact on the performance of the model, particulary at low redshifts, where the difference between M_vir and M_200m becomes significant. The virial mass definition does not appear to be suitable in this context, as no choise of fitted parameters yields a good agreement at the high-mass end of the HMF at z = 0.0. This behaviour is related to the fact that the virial overdensity, Δ_vir, is not constant. In the UCHUU simulation suite, Δ_vir evolves following Bryan & Norman (1998), giving Δ_vir(z = 0)≈102, Δ_vir(z = 1.03)≈159, and approaching the Einstein–de Sitter value of 178 at higher redshift. That work argues that the density contrast at virialisation should depend on Ω_m (and thus on redshift for ΛCDM cosmologies). However, there are strong counterarguments suggesting that Δ_vir should be nearly independent of redshift. A detailed discussion of this point and its theoretical implications will be presented in Betancort-Rijo et al. in prep.; here, we simply highlight that the available evidence points toward a weak or negligible redshift dependence of Δ_vir. In Sánchez-Conde et al. (2007), we investigated the spherical collapse model with shell-crossing to study the stabilisation of proto-halos and no difference was found amongst the Einstein-de Sitter and ΛCDM cosmologies.

If the results of Bryan & Norman (1998) were correct and the stabilisation process truly depended on redshift, then the HMF constructed using M_200m should exhibit explicit redshift dependence. In that case, the formalism presented here - which contains no explicit dependece on z - would be unable to reproduce the simulation results accurately. Conversely, the HMF constructed using M_vir, together with its redshift-dependent Δ_vir, would account for the assumed redshift dependence of the stabilisation process and thus eliminate this source of spurious explicit redshift dependence. Under such circumstances, the HMF based on M_vir would be expected to appear more universal.

However, Figure E.1 shows the opposite behavior: the HMF constructed using M_200m provides a significantly better match to the formalism (which is explicitly independent of z), and therefore yields a more universal HMF.

It should be noted that in Figure E.1, the M_vir data from the simulations are compared with the predictions of the formalism using parameter values slightly adjusted relative to those employed for M_200m, in order to obtain the best possible agreement across the full redshift range. Even with this optimisation, the model fails to accurately reproduce the HMF at low redshift when using M_vir at z = 0, there is a pronounced discrepancy for high-mass end, while for z = 1.03 and larger z’s the difference residues for both mass definitions are within the accuracy of the formalism. Since for these redshifts Δ_vir changes little according to Bryan & Norman (1998) formalism, this support the conclusion that Δ_vir is in fact almost independent of redshift.

Appendix F: Comparison with the Dark Emulator

Fig. F.1. Comparison of halo mass function ratios using two different theoretical predictions. Each panel shows the ratio of the data to the GPS+model prediction and to the Dark Emulator prediction, highlighting relative differences between the two approaches. The panels correspond to redshifts z = 0.0 (left), z = 0.49 (middle), and z = 1.03 (right).

We compare the GPS+ model with the predictions from the widely used Dark Emulator (Nishimichi et al. 2019), which is calibrated only up to z = 1.48 and for M_200m ≳ 10¹²h⁻¹ M_⊙; therefore, the comparison is restricted to redshifts below this limit.

Figure F.1 presents the ratios of the data relative to the HMF predicted by the GPS+ model and by the Dark Emulator at z = 0.0 (left panel), z = 0.49 (middle panel), and z = 1.03 (right panel). The Dark Emulator fails to reproduce the UCHUU HMF at low halo masses and systematically overpredicts the HMF at the high-mass end. Overall, the GPS+ model shows significantly better agreement with the simulation data across the full mass range considered.

All Tables

All Figures

	Fig. 1. Multiplicity function f(M, z) (Eq. (1)) from UCHUU simulations (markers) vs. GPS+ predictions (solid lines) at z = 0–20 (right to left; see Table B.2). The y-axis shows the normalised differential halo abundance including dlnM/dlnσ. Only bins with > 200 halos are shown.
In the text

	Fig. 2. Ratios of the HMF predictions from GPS+ and ST at different redshifts. Only mass bins with > 100 halos are shown.
In the text

	Fig. D.1. Top panel: Multiplicity function f(σ, z) as a function of σ at different redshifts, as predicted by the GPS+ model presented in this work. Bottom panel: Ratio of the model at each redshift to z = 0 prediction.
In the text

	Fig. E.1. Comparison of the ratio between the simulated HMF and the GPS+ model at different redshifts, for two halo mass definitions in the data: the mass enclosed within overdensity 200ρb, M200m, and the mass within the virial radius, Mvir.
In the text

	Fig. F.1. Comparison of halo mass function ratios using two different theoretical predictions. Each panel shows the ratio of the data to the GPS+model prediction and to the Dark Emulator prediction, highlighting relative differences between the two approaches. The panels correspond to redshifts z = 0.0 (left), z = 0.49 (middle), and z = 1.03 (right).
In the text