The abundance of dark matter haloes down to Earth mass

We use the Voids-within-Voids-within-Voids (VVV) simulations, a suite of successive nested N-body simulations with extremely high resolution (denoted, from low to high resolution, by L0 to L7), to test the Press-Schechter (PS), Sheth-Tormen (ST), and extended Press-Schechter (EPS) formulae for the halo abundance over the entire mass range, from mini-haloes of $10^{-6}\ \mathrm{M_\odot}$, to cluster haloes of $10^{15}\ \mathrm{M_\odot}$, at different redshifts, from $z=30$ to the present. We find that at $z=0$ and $z=2$, ST best reproduces the results of L0, which has the mean cosmic density (overdensity $\delta=0$), at $10^{11-15} ~\mathrm{M_\odot}$. The higher resolution levels (L1-L7) are biased underdense regions ($\delta<-0.6$). The EPS formalism takes this into account since it gives the mass function of a region conditioned, in this case, on having a given underdensity. EPS provides good matches to these higher levels, with deviations $\lesssim 20\%$, at $10^{-6-12.5} ~\mathrm{M_\odot}$. At $z \sim 7-15$, the ST predictions for L0 and the EPS for L1-L7 show somewhat larger deviations from the simulation results. However, at even higher redshifts, $z \sim 30$, EPS fits the simulations well again. We confirm our results by picking more subvolumes from the L0 simulation, finding that our conclusions depend only weakly on the size and overdensity of the region. The good agreement of EPS with the higher-level simulations implies that PS (or ST) gives an accurate description of the total halo mass function in representative regions of the universe.


INTRODUCTION
In the Lambda Cold Dark Matter (ΛCDM) model structures grow hierarchically from primordial quantum fluctuations to the galactic haloes and large-scale structures in the universe observed today (Davis et al. 1985).The abundance of these structures as a function of their mass provides a fundamental basis for galaxy formation models (White & Frenk 1991) and a framework for constraining cosmological parameters (Frenk et al. 1990;White et al. 1993;Henry et al. 2009;Allen et al. 2011).It is important, therefore, to develop theoretical models that describe how initial density perturbations collapse into non-linear structures, and to use these as tools for interpreting the predictions of N-body simulations.
Using the statistics of Gaussian random fields and the spherical collapse model, Press & Schechter (1974) proposed the well-known analytical model for halo mass functions, the so-called Press-Schechter formalism (PS model hereafter), which, at least qualitatively, predicts halo abundances that are comparable to numerical simulations.However, the PS model overpredicts the halo mass function in the low-mass regime (by even up to 60% at  = 0) and exhibits too sharp ★ Email: hnzheng@nao.cas.cn† Email: lgao@bao.ac.cn a decrease at the high-mass end (White et al. 1993;Gross et al. 1998;Governato et al. 1999;Jenkins et al. 2001;Lukić et al. 2007;Pillepich et al. 2010).To remedy this, Sheth & Tormen (1999, 2002) replaced the spherical collapse ansatz in the PS model with the ellipsoidal collapse model and provided another formula (ST model hereafter) whose predictions provide a somewhat closer match to those of numerical simulations.For example, according to Reed et al. (2003), the difference between the ST halo mass function and their simulated one is ≲ 10% in well-sampled mass bins.Bond et al. (1991), Bower (1991) and Lacey & Cole (1993) extended the PS model (hereafter the EPS model), to include predictions for the halo abundance and assembly history in different environments (Gao et al. 2005;Faltenbacher et al. 2010).Thereafter, several studies have attempted to provide accurate and universal fitting formulae for halo mass functions from numerical simulations (e.g.Jenkins et al. 2001;Warren et al. 2006;Reed et al. 2007;Lukić et al. 2007;Tinker et al. 2008;Crocce et al. 2010;Manera et al. 2010;Watson et al. 2013;Despali et al. 2016), improving upon earlier fits by extending to higher redshifts and wider mass ranges, and also considering different definitions of halo mass.
It is challenging to perform cosmological simulations with high mass resolution over a wide mass range down to redshift  = 0 because of computational cost.Previous studies, therefore, have fo-cused primarily on either high-mass haloes (10 10−15 M ⊙ ) down to low redshift (e.g.Jenkins et al. 2001;Tinker et al. 2008) or lowmass haloes (10 5−10 M ⊙ ) at high redshifts (e.g. = 10, 30, see Lukić et al. 2007;Reed et al. 2007).As a result, there are comparatively few studies comparing the theoretical halo mass functions (e.g. the PS, ST, and EPS models) to simulations at the mass range of mini-haloes (≲ 10 5 M ⊙ ) down to low redshift (Angulo et al. 2017).Furthermore, most comparisons of EPS predictions with simulations focus on overdense or mean density regions.In this work, we use the VVV simulations (Wang et al. 2020), a series of successive nested zoom cosmological simulations of void regions with extremely high mass resolution, to test the accuracy of theoretical halo mass functions over the full mass range of CDM haloes as a function of time.By construction, this work focuses on the abundance of haloes in preferentially underdense regions.
The paper is organized as follows.Section 2 briefly describes the PS, ST, and EPS theoretical halo mass functions and Section 3 the details of our simulations.Our results, discussions and conclusions are presented in Sections 4, 5, and 6, respectively.An examination of the conversion between linear and non-linear overdensities are illustrated in Appendix A. Further tests with the EAGLE simulations (Schaye et al. 2015) are presented in Appendix B.

THEORETICAL HALO MASS FUNCTIONS
Press-Schechter theory (Press & Schechter 1974) assumes that the initial density field follows a Gaussian random distribution, and that gravitational collapse occurs when the smoothed density field, , exceeds the critical overdensity for collapse,  c , during its random walk in the space of () − , where () represents the variance of the density field smoothed with a filter of mass, , by the conventional real-space top-hat filter (e.g.Lacey & Cole 1993;Mo & White 2002).These crossing events correspond to structure formation on a certain scale, yielding a halo mass function given by: where ρ0 = Ω m  crit,=0 is the mean matter density, (, ) =  c /[ ()()],  c = 1.68647, and  () is the growth factor given by  () = ()/[(0) (1 + )], with (assuming a flat universe, i.e.Ω m + Ω Λ = 1): and  is the scale factor defined as (1 + ) −1 .While PS theory considers a spherical model for the collapse of perturbations, the Sheth & Tormen (Sheth & Tormen 1999, 2002) model takes ellipsoidal collapse 1 into consideration, leading to a modified formula: 1 The spherical collapse model considers the evolution of a spherically symmetric density perturbation; the ellipsoidal model makes the more realistic assumption that the overdense region is an ellipsoid, which leads to different collapse times in different directions. (4) where  ′ = √  ST , with free parameters  ST = 0.707,  ST = 0.322 and  ST = 0.3.
Extended PS theory (Bond et al. 1991;Lacey & Cole 1993;Mo & White 1996) is a useful tool to quantify halo assembly history and assembly bias.In particular, it predicts the probability that matter in a spherical region of mass,  0 , at redshift,  0 , and linear overdensity,  0 , is contained within dark matter haloes of mass in the interval ( 1 ,  1 + d 1 ) at redshift  1 : where, for a virialized structure,  1 =  c / (), and  0 is given by2 : as an approximate solution of Eqs.(16-17) of Mo & White (1996).Equation 6is just an approximation to the exact solution that relates the nonlinear density to the linear density during the collapse of a homogeneous sphere.Here,  nl denotes the non-linear overdensity in Eulerian space, which can be calculated directly from an N-body simulation.We modified the original formula from Mo & White (1996); Sheth & Tormen (2002) by including an additional factor,  ( nl ), to improve the accuracy of the fit when the non-linear overdensity is close to -1.We carried out our own fit to the spherical collapse model and obtained a correction factor of  ( nl ) = 1 − 0.0053977 + 0.00184835 2 + 0.00011834 3 with  = min(0, ln (1 +  nl )).This makes a difference of about 5% for  nl = −0.99.Once these are calculated, we may obtain the corresponding halo mass function as: where a factor of (1+ nl ) is introduced to accommodate the difference of the sphere sizes in Eulerian and Lagrangian space.It is worth noting that, in the limit where the density of a large region is equal to the cosmic mean matter density (i.e.,  nl = 0), the EPS formula reverts to the standard PS formula.

DETAILS OF THE SIMULATIONS
We use the nested zoom N-body simulations of the VVV project (Wang et al. 2020), which consists of 8 levels of resimulation, covering a wide halo mass range spanning around 20 orders of magnitude (∼ 10 −6 − 10 15 M ⊙ ).These simulations were performed with the Gadget-4 code (Springel et al. 2021), adopting cosmological parameters derived from Planck (Planck Collaboration et al. 2014): Ω m = 0.307, Ω Λ = 0.693, ℎ = 0.6777,  s = 0.961 and  8 = 0.829.On large scales ( ≤ 7 Mpc −1 ), the initial linear power spectrum is computed with the Camb code (Lewis et al. 2000), while the BBKS fitting formula (Bardeen et al. 1986) with Γ = 0.1673 dn/d log 10 M/(Mpc) 10 −8 10 −6 10 −4 10 −2 10 0 10 2 10 4 10 6 10 8 10 10 10 12 10 14 10 16 Predicted and simulated halo mass functions in different resolution levels at  = 0.In the top panel, the colour solid lines show the halo mass functions in the VVV (error bars represent Poisson errors); black, gray and colour dashed lines show predictions of the PS, ST, and EPS models for each resolution level respectively.In the bottom panel, the solid lines show the ratio VVV/EPS.Thick lines indicate mass bins containing at least 20 haloes.The black dotted line represents a ratio of unity, while the gray dashed line represents the ratio ST/PS.At  = 0, ST gives the best prediction for L0; EPS overpredicts the halo abundance in higher resolution levels by ∼ 20%.
Table 1.Properties of the different resolution levels in the VVV simulation used in this work.Column 1: name of the level; column 2: size ( box or  sphere ) of the region selected at each resolution level at  = 0 -we use the entire cube in L0, while in L1-L7, we use the central sphere of diameter ∼ 0.8 times the diameter of entire high-resolution region; column 3: mass of highresolution particles; column 4: softening length of high-resolution particles; column 5: overdensity (  level = ρlevel / m − 1) of the selected region at  = 0. and  8 = 0.8811 is adopted to extrapolate the power spectrum to small scales ( ≥ 70 Mpc −1 ), with a smooth transition between 7 − 70 Mpc −1 .
In the VVV, we select candidate resimulation regions to be nearly spherical in shape and underdense (relative to the cosmic mean density); initial conditions for these candidate regions are then generated at higher resolution for each subsequent VVV resolution level.This nested zoom technique is ideal to study extremely small structures embedded within cosmologically representative environments (Jenkins 2010(Jenkins , 2013;;Jenkins & Booth 2013).We refer the reader to Wang et al. (2020) for a detailed description of the zoom-in strategy.Wang et al. (2020) present 8 levels of resolution with uncut initial power spectra 3 , labelled L0-L7.L0 refers to the periodic, 'parent' simulation cube with  box = 738 Mpc.Following Wang et al. (2020), in 3 There are two additional levels (L7c and L8c) where the initial power spectrum is cut off on small scales to reflect the free streaming of a 100 GeV neutralino.Simulations where this free streaming cutoff is resolved are known L1-L7 we only consider the halo population contained within 0.8 times the radius of the high-resolution region, so as to avoid any potential contamination from low-resolution particles in the boundaries of the high-resolution regions.Details of the simulations at each resolution level are given in Table 1.
Haloes were identified using a friends-of-friends (FOF) algorithm (Davis et al. 1985, assuming a linking length,  = 0.2 times the mean interparticle separation).Subhaloes within haloes were identified using the SUBFIND algorithm (Springel et al. 2001;Dolag et al. 2009); both methods are built into Gadget-4.There are many different definitions of halo mass: a basic approach is to adopt the mass of the FOF group ( FOF ) directly, while others may use  Δ , defined as the mass scale within which the average density is equal to Δ times the critical or the mean density of the universe.The overdensity, Δ, may be set in various ways, e.g.200, 200Ω m or Δ  from the spherical top-hat collapse model (Eke et al. 1996;Bryan & Norman 1998).There is no consensus as to which definition is 'best'.Warren et al. (2006) argued that the FOF mass could suffer from a systematic bias for haloes with small particle numbers, while Tinker et al. (2008) calibrated the parameters in a fitting formula with different mass definitions.In this paper, we adopt  Δ=200Ω m ( 200 hereafter), as the definition of halo mass.We consider only central haloes (excluding subhaloes) in the following analysis.

Halo mass function at 𝑧 = 0
We begin by considering the halo mass function at  = 0, shown in Fig. 1, for each of the different resolution levels (L0-L7, solid lines of different colours).We only include haloes containing at least 50 particles.We compare with the corresponding EPS prediction (dashed colour lines) based on the local overdensity and total mass to produce spurious structures (e.g.Wang & White 2007;Lovell et al. 2014).Thus, for simplicity, we exclude these two levels from our analysis.It can be seen that at  = 0, EPS gives relatively more precise predictions for the results of the simulations (within ∼ 20%) at all resolution levels compared to PS (black dashed line).This is especially true at the higher levels (i.e., lower underdensity regions).ST (gray dashed line) gives the best prediction at L0 (the only at cosmic mean density), resulting in the similarity between VVV/EPS (i.e. the ratio between the simulation halo mass function and the EPS prediction) and ST/PS (i.e. the ratio between the ST and PS halo mass functions, gray dashed line truncated at 50 p, L0 in the bottom panel) lines for L0.
Our results from PS and ST for L0 at  = 0 agree with previous studies (e.g.Jenkins et al. 2001;Reed et al. 2003Reed et al. , 2007;;Yahagi et al. 2004;Lukić et al. 2007;Pillepich et al. 2010), suggesting that the ST model is a better approximation than PS for volumes simulated at mean density and low redshift.This is true even when we use  200 rather than  FOF as the definition of halo mass.For example, Reed et al. (2003) studied the halo mass function in the mass range ∼ 10 10.2 − 10 14.2 M ⊙4 and found that PS overpredicts the abundance by around 20% when  < 10 13.7 M ⊙ at  = 0; ST, on the other hand, performs comparatively better in the same mass range -this is consistent with our findings that the halo mass function of L0 at the low mass end ( ≲ 10 14 M ⊙ ) aligns with the ST prediction and is ∼ 20% lower than PS prediction.On the other hand, it is worth noting that while the simulated halo mass function tracks the ST prediction well at the high mass end, the PS model underpredicts the halo mass function of L0 at  ≳ 5 × 10 14 M ⊙ .This is also consistent with the results of White et al. (1993) In the higher resolution level simulations (L1-L7, corresponding to increasingly underdense regions), the EPS model performs significantly better than either PS or ST.A simple comparison can be made with the halo mass function at the high mass end of the L1 volume in Fig. 1 -even with the most moderate underdensity ( = −0.607),both PS and ST significantly overpredict the halo abundance, while EPS provides a relatively accurate prediction.This is to be expected because these simulations focus on regions far below the average density of the universe.Indeed, it is known that the halo mass function depends strongly on local overdensity (Gao et al. 2005;Rubiño-Martín et al. 2008;Crain et al. 2009;Faltenbacher et al. 2010;Tramonte et al. 2017).For example, Gao et al. (2005) find that at  = 49, the halo mass function in a region with  nl = 4.3 has a larger amplitude than that in a region with  nl = 2.8 by a factor of ∼ 4.Only the EPS model takes local environment into consideration when predicting the abundance of haloes.The mass function of L6 (blue line) is less aligned with other levels.This may be due to cosmic variance (i.e. the scatter in halo abundance across volumes of the same size and matter density; we refer readers to Section 4.3 for a detailed illustration).We examine this by comparing the halo mass function in L6 and the corresponding volume in L5 (ensuring the same non-linear overdensity) and find good agreement, which excludes the possibility of problems related to numerical convergence in L6.In general, EPS provides reasonably accurate predictions (par-ticularly at the low mass end) which overestimate our simulated halo mass functions by ∼ 20%.This might due to various reasons, for example, the ambiguous definition of a "virialized" halo in the original PS theory.

The redshift evolution of the halo mass function
After examining the accuracy of theoretical predictions of the halo mass function at  = 0, we now consider how well these models perform at higher redshift.In Fig. 2, we show the halo mass functions at redshifts ∼ 2, 7.8 and 30.For L0, at low redshift, the simulation result is in good agreement with the ST prediction.With increasing redshift, however, better agreement is gradually found with the PS and EPS models (the PS and EPS predictions are rather similar when  = 0 and the enclosed mass is large enough).
Cohn & White ( 2008) discussed the differences arising from using different definitions of halo mass at  = 10: the halo mass function with  FOF in the mass bin 10 8.2−10.7 M ⊙ agrees with the ST prediction at  = 10, while its counterpart with  180 is almost half of the  FOF measurement.Klypin et al. (2011) (using  Δ  , with Δ  the value from the spherical top-hat collapse model) found that the deviation between ST and the simulations becomes larger at higher redshift -up to a factor of 10 in mass bin 10 10.2−11.2M ⊙ at  = 10.When using a similar definition of halo mass, our results agree with the last two studies, suggesting that the halo mass function with mass defined as  200 is lower than the ST prediction and approaches the PS or EPS prediction at high redshift5 .We note that, when defined using  FOF , the halo mass function (not shown) of L0 deviates less from the ST prediction at high redshifts, which agrees with the conclusion of Reed et al. (2003) that the halo mass function obtained using spherical overdensity masses is lower than the mass function based on the FOF mass.
For the higher VVV resolution levels, the results at  ∼ 2 are quite similar to the results at  = 0: the halo mass functions in the simulations are still in fairly good agreement with the EPS prediction.On the other hand, at  ∼ 7.84 the difference between the simulations and the EPS prediction becomes larger, with the ratio of VVV/EPS dropping to ∼ 50% for L2-L7.At  = 30 results from L3-L7 show that the low-mass end of the halo mass function (computed so that haloes contain at least 50 particles and each bin contains at least 20 halo samples) goes back to being aligned with EPS.
In Fig. 3, we plot the evolution of the geometric mean of the ratio VVV/EPS6 in each resolution level.Error bars have been propagated from the Poisson errors of selected mass bins using the following relations: where  and   represent the value and the error of the geometric mean ratio respectively;   and   represent the value of the ratio and the Poisson error in the -th mass bin; and  is the number of mass bins.Fig. 3 suggests that the ratios in all resolution levels follow similar tracks with time: the average VVV/EPS ratios are close to one at the present, drop at intermediate redshifts, and increase again at earlier times.Generally, the higher the level is, the earlier their track intersects the dashed horizontal line at unity, which is likely related to the different characteristic mass scales of haloes in each level.
We find that the average VVV/EPS ratios are close to unity at low redshifts for large masses and, at high redshifts, for low masses.This is consistent with the results of Gao et al. (2005), who studied the halo mass function in denser environments at  = 0 and  = 49, and showed that EPS gives accurate predictions in the mass ranges 10 8.2−10.7 M ⊙ ( = 0), and 10 2.7−5.2M ⊙ ( = 49), confirming that EPS is a good fit at low redshift and extremely high redshift.

Halo mass functions with different realizations
The analyses in the preceding section use only one realization for each resolution level and are thus subject to "cosmic variance".To assess the effect of this variance and test the accuracy of the theoretical models of the halo mass function in general, we picked spherical regions of different sizes and overdensities from the L0 cube, and measured the average halo mass function in them.The L0 volume is large enough to sample a range of overdensities with statistical fidelity.We generated 10 5 randomly located points and measured their local overdensities inside spheres of different radii, and then selected 100 samples with the closest overdensity to a set of preselected values.We then computed the mean halo mass function across these 100 samples; this is displayed in Fig. 4. Black (gray or colored) dashed lines represent the predictions of PS (ST or EPS), and the colour solid lines show the mean halo mass function of samples in the VVV simulation.Different colors represent different overdensities, the values of which are adapted to the redshift and sphere size of each panel (marked at the upper right corner), in order to obtain a large enough sample.
It can be seen that the ratio VVV/EPS (for  200 ) depends only weakly on sphere size and overdensity, but evolves with redshift -it agrees with the theoretical ST/PS line at  = 0 and 2, but deviates at  = 7.84.The mean halo mass function averaged over the random samples is quite similar to that of the full volume, corroborating our conclusion from the full simulation that the ST model provides a good approximation to the halo mass function at low redshifts, but overpredicts it at high redshifts.However, the size of the error bars suggests that even for regions with the same overdensity and total mass, there is considerable variance in the halo mass function.This is especially true in smaller regions and might be part of the reason why the predicted halo abundance in some levels (e.g.L6) are less accurate compared to other levels.In the upper right panel of Fig. 4, we also show the halo mass function of L1 at  = 0 ( = 41 Mpc,  = −0.61)with an orange dash-dotted line.Comparing to its counterpart in L0 (purple solid line), we find that the L1 mass function is a little larger, by ≲ 10%, at  200 ∼ 10 12.5 M ⊙ , showing reasonable convergence in the halo mass function between the two different levels.

DISCUSSION
In the previous section, we have shown that for different subvolumes of the L0 simulation (corresponding to regions with a range of overdensities), the halo abundance in the EPS model differs from the halo abundance in the simulation at a similar level (20 -40% at  = 0, higher at high ) as in the L0 full volume, which has cosmic mean density.We expect this conclusion to hold not just for a relatively small range in halo mass, but also for smaller masses, down to the cutoff in the power spectrum; we are unable to perform similar tests for higher levels due to the limited sample size in zoomed regions.Consequently, the deviations between the EPS predictions and the (L1 -L7) VVV simulations (Fig. 1 and Fig. 2) should approximately reflect the difference at other overdensities, rather than just the underdense regions probed by the VVV simulations.Combining these two arguments, given that the PS formalism is equivalent to integrating the EPS formula over all linear overdensities,  0 , we can use the PS formula to predict the abundance of haloes across the entire range of halo masses at the mean density.
Based on this, in Fig. 5, we present the evolution of halo abundance in the whole universe (at mean density) as predicted by the PS (black) and ST (coral) formalisms.The predictions are similar to those of Mo & White (2002), but we show the differential mass function, and display it over a much more extended mass range than Mo & White (2002) and with updated cosmological parameters.The comparison with the halo mass function in the L0 simulation (red lines) corroborates our previous results that, at low redshift ( ≲ 2) and high masses (10 11−15 M ⊙ ), the ST model provides a better prediction, while the PS model underpredicts the abundance of large haloes (≳ 10 15 M ⊙ ) and overpredicts the abundance of smaller haloes (∼ 10 11−14 M ⊙ ).At higher redshifts, the simulation results lie between the PS and ST predictions.
For the two lowest mass bins (i.e. 10 −6 M ⊙ and 10 −4 M ⊙ ), we show the effect of the free-streaming cut-off for the 100 GeV neutralino with dotted lines, using the cut-off initial power spectrum of Wang et al. (2020).Note that here we have used the sharp -space window function to evaluate () in the theoretical mass functions, as suggested by Benson et al. (2013).This modification provides a

P S S T L0
Figure 5.The evolution of halo abundance within the mass range 10 −6 -10 15 M ⊙ in the whole universe (or mean density regions) predicted by PS (black lines) and ST (coral lines), with the numbers beside each line representing the corresponding log( /M ⊙ ).The dashed lines show the prediction with the uncut initial power spectrum, while the dotted lines take account of the free-streaming effect at 10 −6 -10 −4 M ⊙ .The red lines are the halo mass function in the L0 simulation, with the error bars representing Poisson errors.Within the mass range 10 11 -10 15 M ⊙ , ST is a good approximation at low redshift  ≲ 2, but overpredicts the halo number density at earlier times.The free-streaming effect suppresses the abundance of 10 −6 M ⊙ haloes by approximately one order of magnitude, while barely affecting larger haloes.more accurate prediction of the halo abundance near the cut-off scale compared to a real-space top-hat filter, which instead overpredicts the mass function near the cut-off by re-weighting long wavelength modes.We find that the effect of free-streaming is to suppress the abundance of 10 −6 M ⊙ haloes by one order-of-magnitude, while larger haloes are barely affected7 .This is consistent with previous studies of warm dark matter models (e.g.Schneider et al. 2012;Benson et al. 2013;Bose et al. 2016), which suggest that free-streaming only suppresses the halo mass function near and below the cut-off scale.

SUMMARY AND CONCLUSIONS
In this paper, we have made use of the VVV simulations, a suite of multi-zoom nested simulations at very high resolution, to test the accuracy of the Press-Schechter (PS), Sheth-Tormen (ST) and extended Press-Schechter (EPS) models for the abundance of haloes as a function of their mass and redshift.In particular, this work focuses on the halo mass function at extremely small scales (down to ∼ 10 −6 M ⊙ ) and its evolution to high redshift ( = 30).The resolution levels are labelled L0-L7, corresponding to smaller, higher density regions of progressively higher mean underdensity (below the cosmic mean).The non-linear underdensity of each level is mapped into a linear underdensity using fitting formulas from the EPS formalism.
We find that at  = 0, the ST model provides the most accurate fit to the halo mass function in the L0 volume (overdensity,  = 0; 10 11−15 M ⊙ ) but the EPS model provides the best fit to the higher resolution levels (L1-L7;  < −0.6 at 10 −6−12.5 M ⊙ ), to better than 20% accuracy (see Fig. 1).The results at  = 2 are similar, but there are larger deviations at  ∼ 7 − 15 from the ST prediction in L0 and from the EPS predictions in the higher resolution levels.However, at even higher redshift,  ∼ 30 (10 −6−2.5 M ⊙ ; −0.55 <  < −0.30), the EPS model provides, once more, a good fit to the halo mass functions in the simulations (see Fig. 2 and Fig. 3).We validated our results by selecting regions of different volume and overdensity from the full L0 volume, and find that the VVV/EPS ratio (for  200 ) depends only weakly on region size and overdensity, but increasingly deviates from unity at higher redshifts (Fig. 4).Finally, we tested convergence by comparing the halo mass function in L1 with those in selected subvolumes of L0 spanning a range of sizes and overdensity.
Having demonstrated that the EPS formalism gives a good description of the halo mass function in the biased regions of the VVV high resolution levels, given the arguments in Section 5, it is reasonable to assume that the PS formalism (or the ST formula) also gives a good description of the halo mass function in representative, mean density regions.Thus, we are able to present the actual halo mass function (number of haloes per unit volume) over the entire mass range in ΛCDM, from 10 −6 to 10 15 M ⊙ (Fig. 5).While the mass function at the high mass end is, of course, well known, our study is the first to explore the mini-halo regime (≲ 10 5 M ⊙ ) where the EPS model provides a prediction for the halo mass function, with deviations at the ∼ 20 − 50% level, which could be reduced with further theoretical work.Finally, we note that since we have analyzed dark matter only simulations, the impact of baryons on the haloes is ignored.A more complete study based on hydrodynamics simulations with full baryon physics will be presented in a forthcoming paper.

Figure 2 .
Figure 2. As Fig.1, but for  ∼ 2 (top panel),  ∼ 7.8 (middle panel), and  ∼ 30 (bottom panel).For  ∼ 30, we only show results from L3-L7 as very few haloes have formed at this time in the lower resolution levels.The PS and ST predictions are calculated at the corresponding redshifts, while the simulation results and the EPS predictions are shown at the closest available redshift output in the simulation.At higher redshifts, the EPS prediction deviates more from the simulations, by ∼ 20 − 50%, especially at  ∼ 7.8.

Figure 3 .
Figure 3. Evolution of the geometric mean of the ratio of the mass functions, VVV/EPS, in different resolution levels represented by the different colours.The error bars are the error of the mean value obtained by error propagation; see the main text for details.The deviations peak at  ∼ 2 − 9 (depending on scale), indicating that the EPS predictions fit the simulations best at low and extremely high redshifts.

Figure A1 .
Figure A1.Conversion between non-linear and linear overdensities.The colour of each cross represents the fraction of traced particles in the volume at  = 31.39,indicating the extent to which the measurement of overdensity in that volume is contaminated by background particles; see the text for details.The purple solid line shows the median values for all samples and the red and orange lines the median values of samples with  traced > 0.7 and 0.8 respectively.The black dashed line shows the prediction of Eq. 6.The bottom panel shows the ratio of the difference relative to the predicted value.Eq. 6 provides a reasonably accurate fit to the simulations, especially for samples with high values of  traced .