Simulating the complexity of the dark matter sheet – II. Halo and subhalo mass functions for non-cold dark matter models Free

ABSTRACT

1 INTRODUCTION

While the cosmological parameters of our Universe have been measured to an astonishing accuracy, some key questions about our Universe remain unanswered. One of those key questions is: What is the nature of dark matter? So far, dark matter seems observationally consistent with a perfectly cold, collisionless fluid. However, most particle dark matter candidates are expected to deviate from this fiducial scenario at some scale. Such deviations can for example be caused by the warmth¹ of dark matter, its self-interaction or its quantum nature.

Dark matter cannot only be investigated by possible interactions with known species in a particle physics lab, but also by the imprints that its properties leave on the formation of structures. Different properties will in general cause different behaviour for the detailed small-scale physics. However, at first order, the main effect of these properties on to the formation of structures is the suppression of the power spectrum on small scales (Murgia et al. 2017). The usually considered canonical case for a dark matter candidate with such a small-scale suppression is thermal relic warm dark matter. However, for example, also fuzzy dark matter (e.g. Niemeyer 2020) and many ETHOS models (Cyr-Racine et al. 2016) exhibit a small-scale suppression, just with a different precise shape than canonical warm dark matter. This suppression will already go a long way for explaining a lot of the difference to cold dark matter (CDM) scenarios – for example in the halo mass function. Here, we will refer to any dark matter model with a small-scale cutoff as ‘non-cold’ dark matter, and we aim to understand the general impact of a small-scale suppression in the power spectrum.

There exist several promising probes that can be used to constrain (or possibly detect) non-cold dark matter (NCDM) candidates observationally. These include the absorption pattern in the Lyman-α forest (Viel et al. 2013; Murgia, Iršič & Viel 2018), perturbations in arcs of strongly lensed galaxies (Vegetti et al. 2014, 2018), flux-ratio anomalies of quadruply lensed quasars (Gilman et al. 2020), number counts of Milky-Way satellites (Lovell et al. 2014; Newton et al. 2020), and gaps in tidal streams from satellite galaxies (Yoon, Johnston & Hogg 2011; Banik et al. 2018, 2019). In all of these cases the aim is to distinguish cold from NCDM by either detecting small structures, as expected in CDM scenarios, or by failing to detect them in expected numbers, if the primordial power-spectrum is suppressed on small scales relative to CDM (see e.g. Enzi et al. 2020 for an overview).

All of these observational constraints depend crucially on the ability to predict how the properties and the abundance of small-scale structure depend on the dark matter candidate. For example, the modelling of flux-ratio anomalies simultaneously requires models for the halo mass function, for the subhalo mass function and for the concentration–mass relation in NCDM cosmologies (Gilman et al. 2020) and may even require models for non-halo structures (Richardson et al. 2021). While it is commonly assumed that N-body simulations can produce converged predictions for these quantities, this has never been proven reliably. Wang & White (2007) found that N-body simulations of warm dark matter universes fragment even well below the free-streaming scale into large numbers of artifical haloes. While such artificial haloes can be at least partially removed in post-processing (Lovell et al. 2014; Agarwal & Corasaniti 2015), it is not clear whether this approach leads to the same result as would have been found in the absence of the artefact. For example, the small artificial haloes may alter the accretion history of larger haloes and ultimately modify their concentrations. Since virtually all quantitative constraints on the warmth of dark matter depend crucially on the predictions from N-body simulations, it is important to test N-body results using independent methods.

Recently, such an alternative approach to simulating collisionless cosmological evolution has been developed, so-called ‘sheet-based’ simulation methods (Abel, Hahn & Kaehler 2012; Shandarin, Habib & Heitmann 2012; Hahn, Abel & Kaehler 2013; Hahn & Angulo 2016; Sousbie & Colombi 2016; Stücker et al. 2020). In such simulations an accurate density estimate without discreteness noise is obtained by reconstructing the dark matter sheet in phase space. This allows simulations of cosmologies with a small-scale cutoff in the power spectrum without any artificial fragmentation (Angulo, Hahn & Abel 2013). However, it comes at the cost of intractably complex dynamics inside of haloes – leading either to biased density estimates at their centers (Angulo et al. 2013) or to simulations that cannot be extended until a = 1 because of exploding refinement costs (Sousbie & Colombi 2016). In Stücker et al. (2020) we presented a solution to this dilemma; our ‘sheet + release’ (S+R) scheme combines the benefits of the N-body and sheet approaches, by using sheet reconstructions wherever possible, but switching to an N-body representation where the dynamics becomes too complex (in the inner regions of haloes).

Here, we will present the first large set of simulations carried out with the S+R method. These simulations span a wide parameter space of dark matter models with a primordial power spectrum cutoff. Our goal here is two-fold. First, we aim to understand the physical implications of the cutoff. Of course, these effects have been investigated in previous studies (e.g. Angulo et al. 2013; Lovell et al. 2014; Ludlow et al. 2016), but never with a scheme that is reliable from scales at and below the cutoff up to the dense centers of haloes. Here we can present the first analysis of simulations that are free of numerical artifacts. We use this to test the reliability of previous results based on N-body simulations, and we can point out the conditions under which the application of N-body results is problematic. We note that Colombi (2021) has already addressed a similar question by comparing ColDICE (Sousbie & Colombi 2016) sheet-based simulations with N-body simulations that use a large softening length (of at least one mean particle separation) in a detailed phase space analysis. However, here we are instead interested whether classical N-body simulations with a small softening length – as employed in the majority of previous investigations – can give robust results on the level of coarse-grained summary statistics, such as the halo mass function and the mass–concentration relation.

Further, we develop comprehensive quantitative descriptions of the halo mass function (HMF) and the subhalo mass function (SHMF) in cosmologies with a cutoff. Inspired by the study of Murgia et al. (2017), we capture the different cutoffs that power spectra can exhibit due to the particle nature of dark matter in a small set of parameters. Here, we use two, one indicating the scale of the cutoff and the other its sharpness. From our simulations we infer a simple functional form for the suppression in the (S)HMF as a function of our two parameters. This description can be used to evaluate the (S)HMF for almost any kind of dark matter model which exhibits small-scale suppression. We hope that these results will facilitate observational studies aimed at inferring constraints on warm dark matter as well as on the parameter space of other NCDM models.

This article is structured as follows: in Section 2 we present the simulations that we have run for this article and show how most NCDM cutoffs can be captured in a two-parameter description. In Section 3 we investigate the halo mass function on an absolute scale and quantify the uncertainties that arise due to issues about which objects should be considered haloes. Further, we compare the results of our S+R scheme with those from N-body simulations. In Section 4 we quantify the relative suppression of the (S)HMF and provide simple descriptions that can be rescaled to many different dark matter models. Finally, in Section 5 we have a brief look at the mass–concentration relations in our S+R simulations, comparing them to the relations that have previously been inferred from N-body simulations.

2 SIMULATIONS

We set up a suite of ‘sheet + release’ (S+R) simulations which test the effect of various power-spectrum cutoffs of the kind obtained generically in NCDM models. The goal is to use our novel S+R simulation scheme to create reliable fragmentation-free simulations of such universes and to measure how the halo- and subhalo-mass functions depend on the parameters that define the power spectrum cutoff. We will explain our parameterization of the space of initial power spectra in Section 2.1 and in Section 2.2 we review the main points of the S+R simulation scheme presented previously by Stücker et al. (2020).

We show visualizations of part of our simulations in Fig. 1. Clearly our S+R technique results in a smoother and better defined density field which does not suffer from artificial fragmentation – unlike the corresponding N-body simulation. Note, however, that in high-density regions the simulations are indistinguishable in these images

Figure 1.

Visualization of a region of the different sheet + release simulations presented in this article (top four panels) and of an N-body simulation for the ‘normal’ cutoff case (the bottom panel). Three panels show cases with differing slope β but the same scale for the cutoff, while one shows results for the ‘normal’ slope but a smaller cutoff in mass (larger k_hm). The sheet + release simulations produce a high-quality density field with no sign of the artificial small-scale structures which are quite prominent in the N-body case. A video (showing the 1.75 and 1 kev cases with β = 6) can be found here: https://bacco.dipc.org/ncdm.html.

Open in new tabDownload slide

2.1 Transfer Functions

Figure 2.

Examples of different dark matter models that can be described through our two-dimensional parameter space (β, k_hm). See Fig. A1 for a more detailed version of this Figure that also displays the corresponding best-fitting parameters. Almost all dark matter models with a small-scale cutoff are well approximated by a two-parameter fit of the cutoff scale k_hm and the steepness parameter β.

Open in new tabDownload slide

Therefore, we set up a suite of simulations which explores cutoff functions with fixed γ = 5, but varying parameters k_hm and β. In this article we will mostly focus on four sheet + release simulations (S+R) that explore this space. However, we will also validate our results and extrapolations through a larger set of N-body simulations which include additional different values of k_hm and β and further explore the difference found when using γ = 2 and γ = 10. See Table 1 for a full overview over our set of simulations. Note that we label the simulations e.g. as 1 keV ‘sharp’/‘normal’/‘flat’, indicating that these simulations have a similar cutoff scale to a thermal relic with m_X = 1 keV, but with a sharper or flatter cutoff with respect to the thermal relic case. This approximate correspondence is just pointed out to guide the intuition, but it is not relevant for the results presented in this article.

In Fig. 3 we show the transfer functions of these simulations relative to k_hm.

Figure 3.

A subset of the transfer functions used in the simulations. The solid lines show the transfer-functions used by the S+R simulations. The dotted lines show simulations with other values β that have been run as N-body simulations for validation purposes. The dashed lines vary γ. It seems clear that varying γ at a fixed β and k_hm induces only minor differences in the transfer-function. See also Table 1.

Open in new tabDownload slide

2.2 The sheet + release simulation method

Here we briefly review the main aspects of the sheet + release (S+R) simulation scheme which was first introduced in Stücker et al. (2020). The main purpose of this scheme is to reliably simulate universes with a small-scale cutoff in the power-spectrum – like warm dark matter universes or other NCDM scenarios.

As discussed in the introduction, N-body simulations of warm dark matter universes suffer from severe discreteness effects (Wang & White 2007). These discreteness effects are avoided by ‘sheet’ simulation schemes (Hahn et al. 2013; Hahn & Angulo 2016; Sousbie & Colombi 2016), which use superior density estimates as can be seen in the top panels of Fig. 1 in comparison to an N-body simulation in the bottom panel. However, pure ‘sheet’ simulations come at the cost of intractable dynamics inside of haloes.

The S+R simulation scheme addresses these issues by combining the benefits of N-body and sheet simulation techniques (Stücker et al. 2020). In the S+R simulation scheme, initially all mass is followed through a sheet interpolation approach. However, in regions where the dynamics becomes too complex – that is in haloes – mass elements are ‘released’, which means that their mass is subsequently followed through a set of N-body particles. Further, we have shown in Stücker et al. (2020) that one can discretize the density field in this scheme well through an oct-tree of cubes where the density distribution in each cube can be well inferred through a joint assignment of sheet and N-body mass elements.

The S+R is the first scheme that can reliably simulate warm dark matter universes – all the way from large scale structures without any artificial haloes down to the very centers of haloes. In this paper we present the first set of simulations run with this scheme and can therefore test many of the predictions that have been previously made through N-body simulations (e.g. Lovell et al. 2014; Ludlow et al. 2016; Corasaniti et al. 2017; Lovell 2020; Bohr et al. 2021) through an independent scheme.

2.3 Simulation parameters

For the simulations we choose the following parameters: The boxsize is fixed to |$L = 20 \rm \ {Mpc} / h$| and the softening scale is set to |$\epsilon = 0.5\rm \ {kpc}/h$| for both the N-body and the S+R simulations. For the N-body simulations we use N = 512³ particles and for the sheet simulations, we use N_T = 256³ sheet-tracing particles, which are used to reconstruct the interpolated sheet, and N = 512³ ‘silent’ particles which are turned into N-body particles for released mass elements. Therefore, the S+R simulations have ‘infinite’ mass resolution in regions where the sheet can be reconstructed and have the same mass resolution as the N-body simulations (M_part = 5.0 × 10⁶ M_⊙/h) in released regions. For the structure finding algorithms we only use the 512³ silent particles for the S+R simulations, so that the results can be directly compared with the N-body simulations.

Since the S+R simulations only sample the initial conditions grid at a resolution of 256³ they have a smaller Nyquist frequency (⁠|$k_{\rm {Ny}} = 40.2 \ h^{-1}\, {\rm Mpc}$|⁠) than the reference N-body simulations. To make sure that the cutoff is well resolved we required that the total variance of the density field when only considering the power spectrum up to half of the Nyquist frequency deviates less than |$1 \ \mathrm{ per \, cent}$| from the total variance of the unsmoothed density field. This is why all our S+R simulations have k_Ny/k_h ≳ 4 and why we only use N-body simulations for colder test cases as shown in Table 1. Further, with the mass resolution of M_part = 5.0 × 10⁶ M_⊙/h the half-mode mass is very well resolved in our S+R simulations, by M_hm/M_part ∼ 10³ − 10⁴ silent particles.

3 UNCERTAINTIES OF THE HALO MASS FUNCTION

In this section, we will explore the shape of the halo mass function in our ‘sheet + release’ (S+R) simulations. We show that the mass function in the strongly suppressed regime (⁠|$n_{\rm {X}} / n_{\rm {CDM}} \lt 5 \ \mathrm{ per \, cent}$|⁠) depends sensitively on the definition of what is considered a halo, in particular, on the numerical details of the halo-finding algorithm. Large fractions of the objects detected by Friends-of-Friends (FoF) halo finders and also of those considered bound by the subfind algorithm (Springel et al. 2001) appear unbound when the tidal field is considered in the binding criterion. Further, many of these are unvirialized. We demonstrate these problems in Section 3.1; we introduce the new binding check from Stücker, Angulo & Busch (2021) that also considers the effects of the tidal field in Section 3.2; and we investigate the virialization of haloes in Section 3.3. Finally, in Section 3.4 we consider the implied uncertainties in the abundance of very low-mass haloes, and we demonstrate that these primarily affect halo masses where the suppression exceeds a factor of 20 so that the abundance relative to CDM is well determined throughout the intermediate regime where it is of practical relevance. This suppression factor forms the basis of our quantitative analysis in Section 4. In the current section we will focus on a single simulation, the ‘normal’ 1 keV simulation.

3.1 Unsuitability of traditional mass definitions

Angulo et al. (2013) showed that, when applied to warm dark matter (WDM) simulations, the popular FoF algorithm detects many small-scale over-densities which do not resemble conventional haloes. (This is a short-coming of the halo-finding algorithm and should not be confused with the issue of artificial haloes that reflects a short-coming of the simulation method.) We demonstrate the problem visually in Fig. 4 where we show random examples of FoF-groups in the mass range |$5\times 10^8\, h^{-1}\, {\rm M}_{\odot } \lt M \lt 6\times 10^8\, h^{-1}\, {\rm M}_{\odot }$| identified with a linking length equal to 20 per cent of the mean interparticle separation. For the N-body simulations most (but not all) of the groups in this mass-range are artificial haloes spaced regularly along filaments. However, for the S+R simulations the FoF-groups in this regime look very different. In the bottom panels of Fig. 4 we show the sheet-density estimate for these same regions and it becomes clear that the overdensities detected by the FoF-algorithm do not correspond to haloes, but rather to overdense shell or caustic regions in the outskirts of larger haloes.

Figure 4.

The density field around typical FoF groups. The different columns represent four randomly selected groups in the mass range |$5\times 10^8\, h^{-1}\, {{\rm M}_{\odot }} \lt M \lt 6\times 10^8\, h^{-1}\, {\rm M}_{\odot }$| (∼ 100 particles). The top row displays the 1 keV N-body simulation – typical N-body FoF groups in this mass range correspond to artifical haloes. The central row shows groups in the 1 keV sheet simulations, but just displaying the 512³ silent particles used in the group-finding in order to compare fairly to the N-body case (also 512³). Typical FoF groups in this mass range correspond to caustic structures in the outskirts of massive haloes. The bottom row visualizes the same regions using the sheet-density field, showing that there are, in fact, no haloes at these locations. The colour bar has been clipped at 4e2. Black lines delineate the convex hulls of the FoF groups. The projection depth is 62.5 kpc in all cases.

Open in new tabDownload slide

The problem with the FoF algorithm is its assumption that all connected regions with over-densities above a certain threshold (e.g. ρ/ρ₀ > 80 for l = 0.2; More et al. (2011)) should correspond to virialized haloes. However, in the continuum limit N → ∞ the density becomes arbitrarily high ρ → ∞ near caustics (e.g. White & Vogelsberger 2009; Vogelsberger & White 2011) – which are also abundant within and outside of haloes (e.g. Feldbrugge et al. 2018). This effect can be neglected in most CDM simulations because the density field fragments on all scales and increasing the resolution decreases the density of the smooth component near caustics. N-body simulations with a power spectrum cut-off fragment artificially on small scales, reducing here also the likelihood of linking through caustics. However, from a theoretical perspective, all algorithms based on detecting haloes purely as overdensities are ill-defined and this becomes problematic in simulations with a small-scale cutoff, particularly if these are, like our S+R simulations, devoid of haloes formed by artificial fragmentation.

We make a simple experiment to investigate this further. Using the sheet interpolation method on the 512³ silent particles of our S+R simulation, we resample its density field with an increasingly larger number of particles and then run a FoF group-finder on these new particle sets. The resulting mass functions can be seen in Fig. 5. No matter how high the resolution, these seem to diverge at small masses. The low mass objects in this tail cannot be haloes of any kind, since the mass-resolution of the parent simulation is limited to 5 × 10⁶ M_⊙/h. At much lower masses, the FoF mass function counts apparently disjoint peaks in overdensity which are not virialized quasi-equilibrium haloes. Such overdensities, like caustics, are not persistent material objects (particles move in and out of them over time) and appear to connect smoothly to the cosmic web. An investigation of their statistics thus does not seem particularly useful.

Figure 5.

FoF mass functions for different resampling levels. Each mass function is plotted down to objects with 20 resampled particles. The FoF algorithm identifies small-scale structures that are indeed physical features of the density field, e.g. high-density regions associated with caustics (see Fig. 4) but they are not at all similar to what we typically think of as haloes. The standard procedure of identifying haloes as compact, highly overdense regions does not work on small scale in WDM universes.

Open in new tabDownload slide

A straight-forward idea for cleaning up the mass function is to adopt a structure-finder that also performs a binding check, the subfind algorithm, for example. We do indeed find that many of the more curious looking FoF groups are unbound according to this binding check, which reduces the enhancement of the mass function at small masses. However, a visual inspection of the remaining seemingly bound objects shows that many still appear very unlike ‘normal’ haloes. An example is shown in the first two panels of the second row of Fig. 6. We find this to be a short-coming of the particular binding-criterion used by subfind which incorrectly tags these objects as self-bound, as we will discuss next. Note that similar issues apply to all other structure-finders with a binding check that we know of.

Figure 6.

A few hand-picked examples to show how typical results from the subfind binding criterion compare to those from our boosted-potential binding check. The first column shows the sheet-density around three structures identified as self-bound subhaloes by subfind. The second and last columns show particles that are considered to be members of these structures by subfind and by our boosted-potential algorithm, repectively. The third shows the comoving potential ϕ while the fourth shows a boosted potential obtained by subtracting a uniform gradient term (equation 5). The red point marks the subfind position of the halo and the black point the saddle-point of the boosted potential. The examples, from top to bottom, are: a massive halo which is bound according to both criteria; a caustic structure that is considered bound by subfind, but is actually unbound due to the large external tidal field; a structure which is bound according to both criteria, but which is not virialized, −G/T = 0.52.

Open in new tabDownload slide

To identify which particles are bound to a candidate halo or subhalo, subfind considers the self-potential of the proposed members and discards those that have enough kinetic energy to escape from the group. This binding check neglects any effects due to tidal forces from surrounding material, even though these are expected to be substantial in many situations. Imagine that the self-gravity of a proposed group of particles creates a nearly spherical potential valley. Now, an external tidal field can deform this valley anisotropically – bending the potential landscape downwards in one direction and upwards in another. If the tidal field is strong enough, it can even make the net potential curvature in one direction negative, thereby removing the minimum and the valley. An example of such a potential field is shown in the fourth panel of the second row of 6. Here, the equipotential lines open up along the filament. The tidal field is preventing the selected group of particles from further collapsing along the filament axis. We refer to Stücker et al. (2021) for an in depth discussion of this topic. We will discuss in the next section how to define such potential landscapes and how to use them to define a binding check that properly takes the tidal field into account.

3.2 The boosted potential binding check

In Stücker et al. (2021) we proposed a new binding check that accounts properly for the effects of tidal fields and is no more complex than traditional checks based on the self-potential of the proposed (sub)halo. We review the idea briefly here, but refer the reader to Stücker et al. (2021) for a more detailed discussion.

We visualize this boosted potential in the fourth column of Fig. 6. A comparison with the third column clearly illustrates the advantages of the boosted potential. For well defined haloes, like the one in the first row, a halo boundary is naturally defined as the equipotential passing through the saddle-point of the boosted potential (marked as a black point). This defines what would be commonly be called the ‘tidal-radius’ of the halo, but can more usefully be thought of as bounding surface. Particles that have enough energy to cross this surface can exit the halo and never come back. More examples of boosted potential landscapes around haloes and subhaloes can be found in (Stücker et al. 2021).

The second row of Fig. 6 demonstrates that in regions where the tidal field is strong (like filaments), it is possible to have a large overdensity which appears bound according to its self-potential, yet does not actually correspond to a bound structure. By construction, the gradient of the boosted potential vanishes at |$\boldsymbol{x}_0$|⁠, but this will only correspond to a local minimum if all three eigenvalues of the tidal tensor (which is invariant to boosting) are positive. Here we use the ‘boosted potential binding check’ (BPBC) as introduced by Stücker et al. (2021) to filter the halo-/subhalo-catalogue produced by FoF+subfind. If suitably refined, the boosted potential could be used to define a structure finder based on the comoving potential alone, requiring neither the density field nor any arbitrary threshold, but we reserve this development for later work.

Our current version of the BPBC operates in three steps:

• Determine the reference acceleration |$\boldsymbol{a}_0$| by averaging over a group of particles around a candidate center |$\boldsymbol{x}_0$|⁠.

• Find the critical contour of ϕ_boost where the valley around |$\boldsymbol{x}_0$| merges with a deeper valley.

• Consider all particles inside the critical contour as candidate members of the (sub-)halo, but unbind all those that have enough kinetic energy to escape this surface.

A more lengthy description with additional numerical details can be found in the appendix of Stücker et al. (2021).

In the fourth column of Fig. 6 the critical contours are marked and we have normalized the colour scale so that the limiting potential is white, lower potentials are blue and higher ones are red. The fifth column of the same figure shows the bound particles according to the BPBC – to be compared with the bound objects identified by subfind and shown in the second column. In the top row the object is clearly more extended than what is suggested by subfind. In the second row the subfind object has no corresponding bound object in column five, reflecting the large effect of the tidal field.

We present halo mass functions inferred using different binding criteria in Fig. 7. Since we do not want to be sensitive to differences in mass definition, different lines in each panel are all inferred using the same definition, but in the various cases, we reject objects as unbound if less than |$40 \ \mathrm{ per \, cent}$| of their reference mass is bound according to subfind or to the BPBC.² Differences between the mass functions in each panel are thus caused purely by differences in the halo sets that pass the binding check. For the FoF mass function in the 1 keV WDM simulation (top panel), the filtering causes large differences, while for a CDM simulation (central panel) the differences are much smaller. In the bottom panel we apply the same filtering schemes to the WDM simulation, but now using a different mass definition |$M_{200\rm {c}}$|⁠.

Figure 7.

Top and central panels: FoF mass functions for the 1 keV WDM simulation and for a CDM simulation after filtering by various binding and virialization checks. There is considerable uncertainty arising from the definition of a halo. The CDM mass function is relatively robust to these differences, but in WDM the mass function is both highly uncertain and definition-dependent at low mass, even when no artificial haloes are present. Bottom panel: M_{200, c} mass function in the WDM case filtered by the same set of criteria (for our simulations this is only defined for subfind bound objects). This mass definition is considerably more robust, but still allows definition-dependent variations up to a factor of 5 or so.

Open in new tabDownload slide

Many small-mass FoF haloes are considered unbound by subfind, but the resulting mass function still has an upturn in abundance at the small-mass end. When the BPBC check is applied, many of these objects are rejected as unbound, but a small upturn still remains, while at least the peak of the mass function becomes already clearly identifiable. As we will see in the next section, most of these objects are far from dynamical equilibrium. The bottom panel of Fig. 7 shows that the number of low-mass objects is reduced in all cases if M_200c is used instead of M_fof, and in addition the differences between the cases become smaller.

3.3 Virialization of WDM structures

Angulo et al. (2013) argued that many of the small-mass objects in their simulations, are ‘protohaloes’ – i.e. objects that are in the process of forming and will grow rapidly in their subsequent evolution. Motivated by this, we test here whether the low-mass objects that pass the BPBC correspond to virialized haloes.

We show an example of a structure that is bound according to all binding criteria, but not virialized (−G/T ≃ 0.52) in the last row of Fig. 6. This object might indeed be interpreted as a protohalo that is in the process of collapsing, in agreement with the arguments in Angulo et al. (2013).

3.4 Uncertainties in mass functions

Fig. 7 demonstrates that different halo definitions give very different results for the low-mass end of the halo mass function in WDM but agree reasonably well in the CDM case. Although this may, in part, be because most halo finders have been tested and tuned using CDM simulations, WDM objects in the strongly suppressed regime of the mass function are very different in nature from the clearly bound, nearly spherical quasi-equilibrium haloes that dominate the CDM mass function at all masses. As a result, in the WDM case the halo mass function at low mass depends strongly on the definition of a halo; Fig. 7 should be understood as an indication of the strength of this definition-dependence. In Appendix  B, we provide a similar figure for satellite subhaloes (Fig. B1) showing that a strong definition-dependence is also present in this case. We conclude that answering the question ‘What is the shape of the halo mass function far below the half mode mass?’ requires first finding a clear and precise answer to the question ‘What should we call a halo?’

It may seem disappointing that although our simulations have eliminated the artificial haloes produced by N-body artifacts, they still do not provide an unambiguous measurement of the halo mass function at small masses. However, it is important to realize that this is a result of a change in the physical nature of structure below the cut-off scale. In addition, we note that details of the mass function in this small-mass limit are quite irrelevant for most practical purposes. At the point where the objects found by subfind and the BPBC become qualitatively different, the mass function is already suppressed with respect to the CDM case by about a factor of 20.

This is illustrated in Fig. 8 which plots on a linear scale the suppression of the mass function of our 1 keV WDM case relative to the CDM case for all the combinations of mass definition and binding check discussed above. All cases appear to give equivalent answers in this representation, except for the FoF mass definition with no binding check. Thus our new simulations give robust results for the primary purpose for which they are needed, namely to provide quantitative input to observational and experimental set-ups trying to distinguish WDM from CDM. A similar figure for satellite subhaloes can be found in Appendix  B (Fig. B2) and is consistent with the same conclusion. We therefore concentrate on this relative suppression when discussing different dark matter models in Section 4, where we also use M_200c masses for haloes and M_sub masses for satellite subhaloes despite believing that the BPBC is preferable to subfind for WDM simulations. This is because the two binding checks give similar results for the relative suppression and using the more common definition facilitates reproducibility.

Figure 8.

The suppression of the halo mass function in the 1 keV WDM case with respect to that in the CDM case for the various definitions of halo mass and the various filtering schemes presented in Fig. 7. When plotted on a linear scale, the relative suppression appears quite robust to the details of the halo definition as long as a binding check is imposed. Only the case of FoF masses with no binding check gives significantly discrepant results. For the rest of this article we will investigate the relative suppression using M_{200, c} (for haloes) and M_subfind (for subhaloes), since these are common choices in the literature and are easily reproducible for other researchers. See Fig. B2 for a similar figure for subhaloes.

Open in new tabDownload slide

4 QUANTIFYING THE SUPPRESSION RELATIVE TO CDM FOR GENERIC NON-CDM MODELS

We have seen in the last section that the suppression of the halo and subhalo mass functions relative to CDM can be reliably and robustly inferred in the intermediate suppression regime, |$n_X(M) / n_{\rm {CDM}}(M) \gt 5 \ \mathrm{ per \, cent}$|⁠. In this section, we will introduce a generic expression for this relative suppression that it is suitable for use in future studies to predict the mass function for generic NCDM models.

4.1 Measurements

Figure 9.

Our estimates of halo (top) and satellite subhalo (bottom) abundance suppression ratios as a function of mass in units of the half-mode mass M_hm of equation (3). From left to right, the panels correspond to different sharpness parameters β of the power spectrum cutoff. Error bars are simulation-based jackknife estimates (see text). Solid black lines are fits using equation (16), while the other black lines in the top panel show previously published results.

Open in new tabDownload slide

Fig. 9 shows the abundance suppression ratios to agree well for simulations with different cutoff scales, k_hm, but similar sharpness parameter, β. However, they vary significantly as a function of β. An enhancement at very low masses is evident in the N-body simulations, reflecting the presence of artificial haloes, but above the characteristic scale of these haloes, the suppression ratios estimated from N-body simulations agree well with those from ‘sheet + release’ (S+R) simulations. Indeed the S+R simulations appear very similar to the artificial-fragmentation-free limit of the N-body simulations as estimated, for example, by Lovell (2020). Note that the agreement of simulations with different cutoff scales k_hm also shows that the measurements are relatively resolution and box-size independent, since these simulations have different particle numbers and abundance statistics at the same value of M/M_hm.

4.2 Fits and abundance suppression models

In general, we find that functions of this form do indeed give a very good fit to our estimates of suppression ratios, but, as was also the case when fitting the power spectrum suppression using equation (1), the parameters of this function are highly degenerate – i.e. one can easily find very different combinations {a, b, c} which create very similar functions. This space is thus not well suited to finding regularities as a function of β.

We fit these functions in three different cases using only the S+R simulations:

• β = 1.5 using the ‘1 kev-flat’ simulation,

• β = 2 using the ‘1 keV-normal’ simulation, and

• β = 6 using both the ‘1vkeV-sharp’ and the ‘1.75 keV-sharp’ simulations.

Figure 10.

Suppression mass-scales as a function of β. The dots and crosses indicate estimates based on the halo and satellite subhalo mass functions, respectively, while the lines indicate power-law fits as in equation (16).

Open in new tabDownload slide

With these functions it is straightforward to estimate the mass function for most non-CDM cosmologies. First, one fits the high wave-number suppression of the initial power spectrum using the model of equation (1) with γ = 5. Then one uses equation (16) with the values from Table 2 to infer the mass function parameters |$\lbrace M_{20{{\%}}}, M_{50{{\%}}}, M_{80{{\%}}}\rbrace$| and maps these to {a, b, c} by numerically inverting equations (14a)–(14c). To simplify this procedure we provide a simple python script to do these steps in a code-repository along with some examples.³

We show mass functions ‘emulated’ in this way as black solid lines in Fig. 9. By construction they provide a good fit to our S+R simulations. More significantly, they also fit well to colder N-body simulations, demonstrating that the mass function can indeed be scaled accurately in proportion to M_hm.

4.3 Interpretation and comparison with previous work

It is interesting to see that the mass function cutoff behaves regularly as a function of the sharpness β of the cutoff in the initial power spectrum, so that a two-parameter description is sufficient to represent the mass function shapes over the full range we have tested. We will test this further in the next subsection. It is also interesting to see that the satellite subhalo mass function (SMF) is very similar to the halo mass function (HMF) in the slightly suppressed regime |$n/n_{\rm {CDM}} \gt 50 \ \mathrm{ per \, cent}$|⁠, but is more weakly suppressed at smaller mass. This is probably because the low-mass SMF is built in part from larger mass objects which have been tidally stripped, and so originated on more weakly suppressed scales.

A number of other studies have studied similar issues. Lovell (2020) measured the halo and satellite mass functions from N-body simulations for the thermal relic case. The cutoff function they adopt (with β ≈ 2.2, γ ≈ 4.5) corresponds approximately to our β = 2 case and is therefore shown only in the central panels of Fig. 9. Our results agree very well with the HMF of Lovell (2020), but for the SMF we find good agreement only at low mass end, with a clear difference at weak suppressions where Lovell (2020) find the SMF to be less suppressed than the HMF. It is hard to judge where this difference comes from, but we suggest that it may not be statistically significant, given the relatively small number of high-mass satellite subhaloes on which the SMF is based in each of our two studies.

In addition, we compare with HMF predictions from two different excursion set (ES) approaches (see e.g Press & Schechter 1974; Bond et al. 1991; Sheth & Tormen 1999; Sheth, Mo & Tormen 2001; Hahn & Paranjape 2014). We first consider the ES approach presented by Schneider, Smith & Reed (2013) which uses a sharp k-space filter as a window function. We find that this approach gives good results for the β = 2 case. This makes sense, since this is a good representation of the suppression in the simulations that that the approach was fitted to. However, for the β = 6 case this ES formalism produces much too sharp a cutoff in the HMF. This particular scheme appears to work only in the case for which it was explicitly tuned.

4.4 Validation and high redshift

So far we have shown that the emulation procedure discussed in Section 4.2 provides a good description of the mass functions of simulations with β = 1.5, β = 2, and β = 6 and with substantial variations in half-mode scale. Fig. 11 explores how well the procedure works for other values of β and for |$z$| > 0. For these validation plots we use N-body simulations, since they are computationally cheaper and we have already shown them to agree well with the S+R scheme as long as results are used only above the scale of artificial fragmentation.

Figure 11.

Validation plots for our suppression ratio model for a variety of deviations from the fiducial case. Top left: Varying γ hardly changes the mass function, showing that the cutoff is well described by the proposed two-parameter model. Bottom left: The predicted functions for different values of β. The model works well throughout the interval 1.5 ≤ β ≤ 6. Other panels: HMF and SMF suppression ratios at |$z$| = 0.63 and 1.66. The model remains reasonably accurate for |$z$| = 0.66, but clear discrepancies are visible at 1.66.

Open in new tabDownload slide

The top left-hand panel of the figure demonstrates that if we keep β fixed while varying the parameter γ (which we have fixed to γ = 5 so far), the mass functions barely change. This is because γ has only a weak impact on the initial power-spectrum. It shows that our two-parameter description is sufficient for the problem at hand.

The lower left-hand panel of Fig. 11 shows suppression ratios estimated from simulations with a range of different β values and compares them with predictions from the emulation procedure described above. The predicted versus measured mass functions agree very well. Our emulation appears statistically compatible with all validation simulations for β ∈ [1.5, 6]. Given the validation plots from Fig. 11 and also the tests against N-body simulations with varying effective resolutions from Fig. 10 we estimate that our model of the suppression factors is accurate at the level |$\Delta f \lt 10 \ \mathrm{ per \, cent}$| in this β range.⁴

Fig. 11 also explores whether the model for the HMF and SMF suppression ratios can be applied at higher redshift, in particular at |$z$| = 0.63 and |$z$| = 1.66. We see that for both HMF and SMF the suppression ratios at |$z$| = 0.63 are well described by the model, but at by |$z$| = 1.66 some apparently significant discrepancies are visible in both cases.

5 THE MASS–CONCENTRATION RELATION

Finally, in this section we will briefly consider the mass–concentration relation (MCR) in our S+R simulations. The MCR is crucial for the inference of reliable constraints on the nature of dark matter. For example, Gilman et al. (2020) find that much of the constraining power on the WDM particle mass from flux-ratio anomalies in quadruply lensed quasars comes from its effect on the MCR, since concentrations are significantly reduced at halo masses well above the half-mode mass. There are two interesting questions that we would like to answer here.

• Do S+R simulations and N-body simulations produce the same mass–concentration relation? So far all published results for the MCR are based on N-body simulations. However, it is possible that accretion of artificial haloes could alter the inner profiles of much more massive objects and so affect their concentrations. We can use our S+R simulations as an independent validation of previous results based on N-body simulations.

• Is the model of the MCR from Ludlow et al. (2016) able to predict its variation as a function of the sharpness β of the cutoff in the initial power spectrum? The excursion set model proposed by Ludlow et al. (2016) is widely used to model the MCR in universes with thermal-relic dark matter. However, if it is to be used to constrain other dark matter candidates, it is important to verify that it predicts the MCR reliably as a function of the shape of the small-scale cutoff in the initial power spectrum, parametrized in this paper by β.

The top panel of Fig. 12 compares the MCR’s found from S+R and from N-body simulations. Overall it seems that the different schemes agree fairly well. For 1.75 keV and β = 6, the MCR of the S+R simulation seems to be slightly lower than the N-body case at small mass. The internal scatter in the relation (Δc ∼ 2) is much larger than this difference, but the uncertainties we estimate for the median values do indicate a systematic difference that is not seen in any of the other cases. This is unlikely to reflect the effects of artificial haloes, which we would expect to be more pronounced in the warmer cases. However, it is still rather small, and a more sophisticated analysis of the convergence of halo profiles in the S+R simulations would be required to exclude any numerical effects.

Figure 12.

Top: A comparison of mass–concentration relations between the S+R and N-body simulations. Bottom: The ratio of the median concentrations of the NCDM cases to those of the CDM case in comparison to the model from Ludlow et al. (2016). Error bars are the statistical uncertainty of the median as determined through a jackknife technique. The error bars of the N-body simulations have been slightly displaced horizontally to avoid overlap. Our S+R simulations independently confirm results from previous studies that were based on N-body simulations.

Open in new tabDownload slide

Overall, our S+R results are in good agreement with the results from N-body simulations. We find no evidence that artificial haloes modify the MCR, thus validating previous results obtained with N-body schemes. We note that Colombi (2021) has shown in a full phase-space analysis that N-body simulations give robust results when benchmarked against the ColDICE sheet-based simulation scheme (Sousbie & Colombi 2016), at least when a large softening (exceeding one mean inter-particle separation) is used. While such a large softening would be generally desirable, it is hard to implement in practice, since a gigantic particle number would be required to resolve the centres of haloes. Our results suggest that N-body simulations with a small softening may give correct results for coarse-grained summary statistics like density profiles, mass–concentration relations, and halo mass functions provided that the applicable range is carefully tested. Our S+R simulations establish that these statistics are not biased by artificial haloes, which are present in N-body simulations with a small softening but not in S+R simulations, but artifacts caused by discreteness effects inside larger haloes cannot be excluded since these would be present in both types of simulation.

In the bottom panel of Fig. 12 we plot the ratio of the median-concentrations of the NCDM and the CDM universes for our S+R simulations and compare them to the predictions of the excursion set model of Ludlow et al. (2016). Note that the high mass end might be affected by finite-box-size effects, yet the measured ratio between NCDM and CDM seems to be well described. For the 1 keV cases we find excellent agreement. The model of Ludlow et al. (2016) seems to capture the differences between the models as a function of β very well, and it even reproduces the point where the concentrations cross at around 3 × 10¹¹M_⊙. This is quite remarkable, since Ludlow et al. (2016) developed their model specifically for thermal relic universes (β ≈ 2), and its use for other cases has never been tested (as far as we know). Only for the 1.75 keV β = 6 case we again find that the MCR of the S+R simulations is below the prediction. As noted above, it is hard to judge the significance of this deviation, since this it is much smaller than the scatter around the MCR.

We conclude that our S+R simulations independently confirm the validity of the mass–concentration relations that have been estimated. previously from N-body simulations. Further, we find that the model of Ludlow et al. (2016) can be used to predict the MCR in universes with a cutoff that has different shape than in the fiducial thermal relic case, β ≈ 2.

6 CONCLUSIONS

In this paper we have presented a first set of ‘sheet + release’ (S+R) simulations of NCDM universes, specifically models in which the initial linear power spectrum is strongly suppressed on small scales. These simulations are free from the artificial fragments that appear in classical N-body simulations of such models and are able to follow the dynamics even in the centres of haloes where previous ‘sheet’ simulation schemes break down.

With these simulations we investigated the objects that form in the strongly suppressed, low-mass tail of the halo mass function, finding them to be quite different from traditional haloes. The majority are elongated caustic overdensities that appear to be gravitationally unbound when the tidal field is considered in the binding check. In addition, we find objects which appear bound but not virialized, which may be protohaloes in an early stage of collapse, as suggested by Angulo et al. (2013). From this, we conclude that the halo mass function in the strongly suppressed regime |$n_X(M)/ n_{\rm {CDM}}(M) \lt 5 \ \mathrm{ per \, cent}$| cannot be inferred reliably, since it depends too strongly on the precise definition of a halo. Any potential observable sensitive to the mass function in this regime would need to be inferred directly from the simulations, without specific assumptions about density structure or dynamical state, since for at least some observational probes the effect of the large number of caustic overdensities may surpass that of the few traditional, quasi-equilibrum haloes at low mass (see e.g. Richardson et al. 2021 for the case of gravitational lensing).

However, it is rather unlikely that the strongly suppressed tail of the mass function will ever be relevant for distinguishing between non-CDM candidates, since more easily measurable differences already arise in the regime of intermediate suppression where |$n_X(M)/ n_{\rm {CDM}}(M) \gt 5 {{\%}}$|⁠. In the second part of this article we have focused on a quantitative description of halo and subhalo mass functions in this regime, finding that for most NCDM candidates the suppression of the primordial power spectrum can be adequately represented by just two parameters, namely k_hm, which describes the spatial scale of the cutoff, and β, which describes its sharpness. We provide fitting formulae for the halo and subhalo mass functions as a function of these two parameters. These are applicable to a wide variety of non-CDM models, provided that the small-scale suppression of the initial power spectrum is well described by our two-parameter model.

Finally, we have found that the mass–concentration relation for low-mass haloes is very similar in our S+R simulations and in standard N-body simulations, agreeing well with the model proposed by Ludlow et al. (2016). Further, we have tested this model for the first time for power spectrum cutoffs that are different in shape from that in the standard thermal relic case. We find that, in these cases also, it predicts the behaviour of the mass–concentration relation astonishingly well for the range in sharpness parameter, β, that we have tested.

Overall, we find good agreement between our S+R simulations and classical N-body simulations for the summary statistics we studied. We can confirm that the N-body scheme reliably predicts both the mass function and the mass–concentration relation in NCDM scenarios when evaluated above the mass scale of artificial fragmentation (see Wang & White 2007) and above the most heavily suppressed regime of the mass function.

ACKNOWLEDGEMENTS

We thank Thomas Richardson for helping with evaluating and testing the simulations. JS thanks Marcos Pellejero-Ibañez and the BACCO group for interesting and motivating discussions. JS and RA acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program with grant agreement No. 716151 (BACCO). OH acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 679145, project ‘COSMO-SIMS’). The authors thankfully acknowledge the computer resources at MareNostrumIV and technical support provided by the Barcelona Supercomputing Center (RES-AECT-2019-3-0015).

DATA AVAILABILITY

The data and algorithms underlying this article will be shared on reasonable request to the corresponding author. We provide some python scripts to facilitate the evaluation of the inferred mass functions in the ncdm-mass-functions repository under https://github.com/jstuecker/ncdm-mass-functions.

Footnotes

REFERENCES

APPENDIX A: TRANSFER FUNCTIONS FITS

In Fig. A1 we show a number of examples of dark matter models with a small-scale cutoff in the power spectrum that we have tried to represent with our two-parameter fitting function. These models are taken from Murgia et al. (2017), who show that their power spectum cutoffs are reasonably well approximated by the three-parameter (α, β and γ) model of equation (1). As explained in Section 2.1) and demonstrated explicitly in this figure, effectively equivalent behaviour is achieved by fixing γ = 5 so that only two parameters remain. Only in the mixed dark matter case is γ ≠ 5 clearly required. For simplicity, the two- and three-parameter fits are matched in these plots by requiring agreement on the scales where T² is 0.5 and 0.85. With one exception, the models considered by Murgia et al. (2017) are almost equally well fit with two parameters as with three.

Figure A1.

Transfer functions from Murgia et al. (2017) when parameters are reduced from 3 (α, β, γ) to 2 (α, β). Solid lines correspond to three-parameter models and dashed lines to two-parameter fits with γ = 5. The labels give α, β and γ with α in units of |$\rm {kpc}/h$|⁠. In all but the mixed dark matter case, the three-parameter model can be adequately reproduced varying just α and β. For more detailed discussion of these dark matter models see Murgia et al. (2017).

Open in new tabDownload slide

APPENDIX B: SUBHALOES

In this section, we show some complementary plots to the ones presented in Section 3, but for the case of satellite subhaloes instead of haloes.

Fig. B1 shows how the satellite subhalo mass function depends on mass definition and filtering criteria. It is directly analogous to the halo mass function plots of Fig. 7. Note that here we plot on the abcissa either the bound mass according to subfind, or the bound mass according to the BPBC. The counts of numbers in each bin are either unfiltered (in which case the subfind mass is on the x-axis) or are filtered to require non-zero BPBC mass and various levels of virialisation (in which case the BPBC mass is on the x-axis). In the CDM case there is an almost constant offset between the subfind and BPBC mass functions. For satellite subhaloes (which are typically subject to strong tidal fields) the BPBC typically finds fewer bound particles than subfind – on average the ratio is about 0.5. Note that this is in contrast to haloes where the BPBC often finds larger masses than subfind, since it does not depend on the FoF pre-selection.

Figure B1.

The satellite mass function for different filtering criteria / mass definitions in a WDM universe (top) and for CDM (bottom). In the CDM universe all of the mass definitions agree well, with the exception that there is an offset between the subfind bound mass and the tidally bound mass (inferred through the BPBC). For the WDM case different filtering choices lead to very different results at the low-mass end. For example the unfiltered subfind catalogue shows a steep rise at small masses whereas other filter choices show a flattening or even a decrease.

Open in new tabDownload slide

For satellite subhaloes the lower end of the mass function depends strongly on the mass definition. The subfind mass function turns up at small masses which is primarily due to unbound overdensities rather than actual subhaloes as shown by the difference between subfind and BPBC masses. When visually inspecting these objects, we found that most of them correspond to caustic-like structures in the outskirts of haloes. However, these effects are dominant in a mass range where the mass function is already suppressed by a factor of around 20 and are therefore probably irrelevant for quantitative studies which try to distinguish between different dark matter models.

In Fig. B2 we show the relative suppression of the satellite mass function for the different mass definitions and catalogue filterings, in direct analogy to Fig. 8. As for the haloes, it becomes clear that the differences between the various cases appear small when considered on a linear scale and are substantial only where the SMF is suppressed by more than a factor of 20. Note, however, that the subfind mass function suppression does not seem to approach zero at small masses (corresponding to the upturn in Fig. B1). We consider this a failure of the mass definition to distinguish caustic overdensities from bound substructures.

Figure B2.

The relative suppression of the satellite subhalo mass function with respect to CDM for the various mass definitions presented in Section 3. The relative suppression in the regime |$n/n_{\rm {CDM}} \gt 5 \ \mathrm{ per \, cent}$| is quite robust to the details of the mass definition. Compare with Fig. 8 for haloes.

Open in new tabDownload slide