The locations of halo formation and the peaks formalism

Abstract

INTRODUCTION

Where and when do dark matter haloes form? The problem of identifying the locations where gravitational collapse leads to bound haloes of dark matter, and predicting the cosmic time at which this will occur, is among the oldest problems of cosmic structure formation theory. The idea that small perturbations in the primordial matter density formed the seeds of the large-scale structure we observe at present is among the cornerstones of our current picture of the evolution of the Universe. An understanding of the relevant processes and a robust theoretical model enables us to map properties such as the abundance and clustering of dark matter haloes – which are directly tied to the corresponding observed properties of galaxies – to well-understood statistical properties of the initial dark matter density.

Although this problem can now be tackled directly using numerical simulations of large, cosmological volumes, it is still important to explore analytical approximations and identify the key physical features that decide the sites of halo formation. The main motivation behind this exercise is to gain a better understanding of the physical processes that affect structure formation in the Universe. From a practical viewpoint, however, this can also lead to useful, fast approximations to the halo mass function, clustering and predictions of collapse time for a given patch in the initial conditions. The latter especially could be useful from the point of view of ‘semi-analytic’ mock catalogue algorithms such as pthalos (Scoccimarro & Sheth 2002), pinocchio (Monaco, Theuns & Taffoni 2002; Monaco et al. 2013), cola (Tassev, Zaldarriaga & Eisenstein 2013), alpt (Kitaura & Heß 2013), etc., which are becoming increasingly popular in the construction of covariance matrices in current and upcoming surveys such as Baryon Oscillation Spectroscopic Survey (BOSS; Manera et al. 2013), WiggleZ (Marín et al. 2013), Euclid (Laureijs et al. 2011) and others.

Our focus in this paper is on the mass function of dark matter haloes, which is the most basic diagnostic of the fully non-linear density field. Analytical descriptions of the halo mass function have traditionally used two parallel approaches: the excursion set approach (Press & Schechter 1974; Epstein 1983; Peacock & Heavens 1990; Bond et al. 1991; Lacey & Cole 1993; Sheth 1998; Sheth, Mo & Tormen 2001; Maggiore & Riotto 2010; Achitouv et al. 2012; Musso & Sheth 2012, 2013; Paranjape, Lam & Sheth 2012) and the peaks formalism (Bardeen et al. 1986; Bond 1989; Appel & Jones 1990; Manrique et al. 1998; Hanami 2001), both of which aim to characterize the locations of collapse in the initial conditions using some criteria. The former relies on counting sufficiently overdense regions in the initial conditions, which it maps to collapsed haloes in the final, gravitationally evolved density field, while the latter associates haloes specifically to peaks in the initial matter density. In other words, while both approaches rely on the statistical properties of the initial conditions to predict final halo abundances, the excursion set approach does this by treating all locations in the initial conditions on the same footing, while the peaks formalism treats density peaks as being special.

The key aspect of the excursion set approach (Bond et al. 1991), which is missing in the traditional peaks approach (Bardeen et al. 1986), is that it explicitly accounts for the so-called ‘cloud-in-cloud’ problem which avoids overcounting overdense regions embedded in larger overdense regions as individual objects. The ‘peak-patch’ approach of Bond & Myers (1996) is a numerical prescription for unifying the two approaches to solve the cloud-in-cloud problem for peaks, or, equivalently, to study excursion sets for a special subset of initial locations, namely peaks. Recent work (Musso & Sheth 2012; Paranjape & Sheth 2012) has shown that this can also be achieved analytically by making some simple but accurate approximations (see also Bond 1989). There are several motivations for doing so (Sheth, Mo & Tormen 2001; Paranjape & Sheth 2012), not least the fact that N-body simulations of cold dark matter (CDM) show that a large fraction of haloes do, in fact, originate from initial density peaks (Ludlow & Porciani 2011a). Further, Paranjape, Sheth & Desjacques (2013) showed that this unified analytical formalism of excursion set peaks (ESP) gives a self-consistent description of the CDM halo mass function as well as clustering which is accurate at the ∼10 per cent level.

It is worth asking whether this formalism has correctly captured all the relevant aspects of structure formation that affect the mass function. One way of addressing this issue is to apply the same formalism in an ‘extreme’ situation which it was not explicitly built to describe. Structure formation from an initial matter power spectrum with highly suppressed small-scale power, as found in warm dark matter (WDM) cosmologies, offers the perfect playground. The reason is that, apart from having a truncated initial power spectrum, simulations of WDM in fact solve exactly the same problem as those of CDM: the evolution of a cold, collisionless, self-gravitating fluid.

Analytically, one then expects that the same ESP expressions, which correctly describe the CDM mass function and clustering, should work for the WDM case as well, with the simple replacement of the CDM initial power spectrum with that of WDM. In this regard, as we describe in detail below, the ‘out-of-the-box’ ESP calculation does considerably better than traditional top-hat-filtered excursion sets: it correctly predicts a turnover in dn/d ln m at the correct scale whereas the latter predicts a monotonic rise at low masses. We will see, however, that ESP predicts a power law decrease at low masses dn/d ln m ∼ m^2/3 which is incompatible with the results of simulations. This analytical prediction is very robust and hints at a missing physical ingredient in the excursion set logic.¹

Our goal in this paper is to characterize the collapsed objects identified in a WDM simulation in terms of the properties of the initial density field. This will allow us to understand the reasons behind the mismatch of the measured mass function and the ESP prediction. The paper is structured as follows.

In Section 2, we describe the numerical simulation and halo finding algorithm, which are the same as presented by Angulo, Hahn & Abel (2013). We then compare the resulting halo mass function with theoretical expectations based on the ESP formalism, and discuss possible reasons for the differences we see between the theory and numerics. To better understand where haloes form, in Sections 3 and 4 we turn to an in-depth analysis of the initial conditions of the simulation. In Section 3 we analyse the initial density field at the ‘Lagrangian patches’ of the haloes (i.e. the initial locations of groups of particles that will eventually be identified as haloes) and demonstrate that all haloes in the simulation are consistent with forming near peaks of the initial density. We also explore correlations between the initial overdensity and shape of the Lagrangian patches, and use these results to motivate the construction of an empirical catalogue of ‘ESPeaks’, which we describe in Section 4. These ESPeaks are a numerical realization of what the ESP calculation aims to accomplish, and we compare their Lagrangian properties with those of the haloes.

Our main conclusion from this exercise is that, while the ESP calculation on average correctly identifies the locations of halo formation, it systematically underpredicts the mass of the resulting object, and that this effect is especially enhanced at halo masses that are small compared to the characteristic mass scale where the WDM mass function turns around. In Section 5 we argue that this mass mismatch is related to a systematic overprediction of the time of collapse of a given perturbation, and propose a modification to the ESP calculation to re-assign masses by correcting for this effect. We show that the resulting mass function not only agrees very well with the WDM result, but also describes the CDM mass function accurately with the simple replacement of the WDM power spectrum with that of CDM.

We close with a summary and discussion in Section 6. The appendices collect technical details and arguments used to motivate some of the results in the main text.

THE HALO MASS FUNCTION: CONFRONTING SIMULATIONS AND THEORY

Matter power spectra with an initial small-scale truncation arise naturally in warm and hot dark matter cosmologies where density fluctuations on small scales are suppressed due to the late transition to the non-relativistic regime of the respective dark matter particle. Such a power spectrum leads to a corresponding turnover in the late time halo mass function. The numerical determination of such mass functions has, however, proven extremely challenging due to the presence of low-mass objects that arise – completely unphysically – from the fragmentation of filaments (see e.g. Avila-Reese et al. 2001; Bode, Ostriker & Turok 2001; Melott 2007; Wang & White 2007; Hahn, Abel & Kaehler 2013). In the presence of artificial fragmentation, mass functions can only be measured indirectly after filtering or correcting for the spurious haloes (see e.g. Lovell et al. 2012, 2013; Schneider et al. 2012). Only more recently has the behaviour of the halo mass function around and below the turnover scale been explicitly demonstrated by Angulo et al. (2013, hereafter AHA13).

While such WDM cosmologies are of course of genuine physical interest in their own right, we are mainly concerned with a different aspect here: the suppression of low mass haloes provides powerful additional leverage to test models of structure formation in such cosmologies. The exponential fall in the halo mass function at large masses – whose sensitivity to cosmological parameters has been exploited for decades – is replicated here at the small-mass end. Any theoretical model must now describe both of these strong features in the mass function.

We begin by briefly discussing the numerical simulations of AHA13 and the WDM halo mass function they measure. We will see that these numerical results do not meet theoretical expectations based on the ESP formalism. We discuss possible reasons for this, which will motivate our subsequent analysis.

Numerical simulation

The numerical simulation discussed by AHA13 employs the novel t4pm method (Hahn et al. 2013), which completely suppresses artificial fragmentation and allows the determination of the halo mass function at and below the turnover scale in the absence of numerical artefacts.

Specifically, this simulation resolves a 80 h⁻¹ Mpc cosmological volume with 1024³ particles with cosmological parameters Ω_m = 0.276, |$\Omega _\Lambda =0.724$|⁠, Ω_b = 0.045, h = 0.703, σ₈ = 0.811 and n_s = 0.96, consistent with the 7-year Wilkinson Microwave Anisotropy Probe (WMAP7) data release (Komatsu et al. 2011). The normalization of the power spectrum using σ₈ was set using a CDM spectrum, so that the amplitude of fluctuations on large scales is independent of the truncation scale of the power spectrum.

Adopting the fit of Eisenstein & Hu (1999) as the fiducial CDM transfer function T_CDM, initial conditions were generated using the music code (Hahn & Abel 2011) at an initial redshift of z = 63 using the Zel'dovich approximation. We note that a simulation initialized at such a rather low redshift using first-order Lagrangian Perturbation Theory is to some degree affected by transients from the initial conditions (e.g. Crocce, Pueblas & Scoccimarro 2006). The high-mass end of the halo mass function is thus expected to deviate from the true one. The small volume of the simulation adds further to a systematic deviation. Furthermore, it is possible that the detailed behaviour of the mass function around the half-mode mass is also affected by transients. This possibility needs to be considered for precision determinations of the halo mass function but we do not expect it to alter the qualitative behaviour with which we are mostly concerned here (the results of AHA13 are roughly consistent with e.g. the predictions of Schneider et al. 2013 who use second-order Lagrangian Perturbation Theory). We note that the half-mode mass is resolved with almost 100 000 particles in the simulation of AHA13.

Halo identification

AHA13 found that the suppression of artificial fragmentation leads to a failure of the Friends-of-Friends (FoF) halo finder. The dense cores of filaments (in the absence of artificial fragmentation) lead to a percolation of large regions of several haloes when the standard linking parameter b = 0.2 is used. Instead, they first adopted a linking parameter of b = 0.05 times the mean interparticle separation and then determined the spherical-overdensity (SO) mass centred on the centre-of-mass of the parent FoF group. A halo was defined as the sphere of radius R₂₀₀, which has a mean density of 200 times the critical density, ρ_crit. This corresponds to a halo mass |$M_{200} = (4\pi /3)R_{200}^3(200\rho _{\rm crit})$|⁠.

Further, by analysing all haloes individually, AHA13 found that the halo sample could be divided into various subsamples or ‘types’. ‘Type-1’ objects are virialized haloes, while ‘type-2’ include haloes in late stages of formation; the latter do not show an isotropic density structure and instead contain larger scale caustics that are remnants of their formation. The remaining objects were haloes in early stages of formation that have e.g. just started collapsing along the third axis.

In what follows, we only consider the ‘type-1’ objects clearly identified as haloes. The red histogram in Fig. 1 is identical to the line labelled ‘haloes’ in fig. 7 of AHA13 and shows the mass function of these ‘type-1’ haloes, which has a sharp cut-off between 10¹¹ and 10¹² h⁻¹ M_⊙. (The other two histograms will be discussed in Section 4.) It is important to note here that the classification was performed visually and we thus expect that, on an object-by-object basis, the distinction between types 1 and 2 is likely not perfectly robust.

Figure 1.

Halo mass functions: the abundance of virialized objects in the simulation (red) as well as the mass functions of ‘ESPeaks’ determined using empirical excursion set walk starting at peak locations from the same initial conditions as the simulation using the barrier in equation (28) (orange; see Section 4.1 for details). The dashed blue histogram shows the result of the same algorithm but now using the deterministic barrier B = δ_c + 0.5σ₀. The solid black line shows the ‘standard’ analytic ESP calculation from Paranjape et al. (2013) which used a stochastic barrier adjusted to match CDM simulations. The dashed black line shows the ESP calculation using the deterministic barrier mentioned above (see Section 5 for details). Although peaks theory correctly predicts a turnover in the mass function at the correct scale, it gets the asymptotic behaviour below this scale wrong. For comparison, the dotted black line shows the result of using the Tinker et al. (2008) fitting formula.

Open in new tabDownload slide

Theoretical expectations

As discussed earlier, as far as late time structure formation is concerned, the only difference between the CDM and WDM cosmologies is the lack of initial small-scale power in the latter. So one might expect that a physically motivated description of CDM structure should apply equally well to WDM – at least at scales much larger than the free-streaming scale – with the simple replacement of the CDM transfer function with the one in equation (1). One can already anticipate that the standard hierarchical excursion set calculation would have difficulty in describing the low-mass end of the WDM mass function where the effective logarithmic slope of the power spectrum n_eff ≡ d ln P(k)/d ln k becomes steeper than −3. This is the boundary beyond which hierarchical prescriptions are known to fail, with predicted mass accretion rates becoming ill-defined (Lacey & Cole 1993, 1994). The ESP formalism, however, introduces a new ingredient into the picture – the peaks constraint – and since we know that it works well for CDM, we can ask what it predicts for the WDM case.

It is easy to see that the spectral integrals σ_j, and consequently also ν, γ and V_*, will each approach a constant value for WDM as m → 0. The behaviour of ν and V_* in this respect is very different from that in CDM in the same limit, where ν, V_* → 0 as m → 0. This reflects the fact that peaks can only form on scales large enough to be inhomogeneous; reducing the smoothing scale cannot wipe out existing peaks for any power spectrum and, for a truncated spectrum, cannot introduce new peaks. This can be seen in the top panel of Fig. 2: for each choice of m_dm, V_* at large m is close to the CDM value (and hence approximately proportional to |$V=m/\bar{\rho }$|⁠, where |$\bar{\rho }$| is the comoving mean density), but deviates as m approaches the half-mode mass M_hm, finally approaching a constant value at m ≪ M_hm, this value being close to |$M_{\rm hm}/\bar{\rho }$|⁠. This last aspect gives us an interesting physical interpretation of the half-mode mass as being essentially the same as the asymptotic peaks scale.

Figure 2.

Top panel: characteristic peak volume V_* (solid curves) defined in equation (8), for CDM (black) and three choices of WDM particle mass m_dm; from bottom to top: 1 (blue), 0.5 (yellow), 0.25 keV (red). Dotted black line shows the Lagrangian volume |$V=m/\bar{\rho }$|⁠. For CDM, V_* closely tracks V. For WDM, as m decreases, V_* tracks V until m ∼ M_hm (vertical dotted lines) and then approaches a constant value close to |$M_{\rm hm}/\bar{\rho }$| (horizontal dashed lines). All the V_* curves used Gaussian filtering as described in Appendix A. Bottom panel: the relation between σ₀ and mass m with top-hat filtering for CDM (black) and when using equation (1) with the same particle masses (and colour coding) as in the top panel. For WDM we clearly see a ‘freezing-out’ of σ₀(m) at small masses.

Open in new tabDownload slide

The behaviour of V_* discussed above shows that the turnover occurs at around m ∼ M_hm. The solid curve in Fig. 1 shows the ESP calculation of Paranjape et al. (2013) using the WDM transfer function (1). The dashed curve shows the ESP result with a somewhat different, convenient choice of barrier for the random walks, which we will discuss in detail in Section 5.2. The main point is that both of these curves show the turnover and the asymptotic scaling of ∼m^2/3. As discussed above, the latter is a very robust prediction that cannot be easily altered by technical modifications within the framework of a single barrier.

Theory versus simulation – what could be going wrong?

Although the physical requirement of being a density peak naturally accounts for a turnover in the mass function at the correct mass scale, the asymptotic scaling (very robustly) predicted by ESP is clearly wrong. There are several issues which could in principle affect this result.

1. Dynamics. A dramatic possibility is that, since small-mass haloes in WDM do not form hierarchically (e.g. at some point in time the first object forms, with no virialized progenitor), the peaks calculation might simply not be applicable. This would lead to the interesting question of just what it is that characterizes the locations and dynamics of the collapse of small-mass objects. The cut-off scale in the initial spectrum could in principle allow for higher order catastrophes (cf. Arnold, Shandarin & Zeldovich 1982) to become relevant, and these need not necessarily appear as peaks when filtered on the protohalo scale. In CDM, the situation is quite different in this respect, since fluctuations persist down to very small scales and so every protohalo has a progenitor at a smaller scale.

2. Barrier shape. It has been argued that the collapse barrier appropriate for WDM haloes is very different from the corresponding CDM one due to thermal effects in WDM, and that this can introduce a sharp cut-off in the mass function (Benson et al. 2013). However, since WDM simulations see a cut-off despite ignoring thermal effects (e.g. AHA13; Schneider et al. 2013), the origin of the cut-off must be rooted in the suppression of initial small-scale power, and must therefore be a generic feature of cold collisionless dynamics in such conditions. The failure of excursion set (peaks) models to reproduce the correct mass function hence indicates quite clearly that these models are still not accounting for some important physical processes. One of the primary goals of this work is to investigate the cause of this behaviour.

3. Patch shape. Traditional excursion sets, as well as the ESP calculation of Paranjape et al. (2013), use spherical filters when assigning masses to objects, and it could be that asphericity of the Lagrangian patches affects the mass assignment significantly at small masses. For example, recently Despali, Tormen & Sheth (2013) have demonstrated using CDM simulations that accounting for halo asphericity using an ellipsoidal halo finder can lead to small increases in mass for low-mass haloes (see also Ludlow & Porciani 2011b).

4. Stochasticity. Regardless of the importance of thermal effects, the specific details of the barrier, e.g. those related to stochasticity in the barrier height, are in fact somewhat uncertain (even in the CDM case). The ESP calculation for CDM is self-consistent but not fully predictive, and needs some inputs from simulations (Paranjape et al. 2013). In particular, the barrier used in that calculation was adjusted to match measurements by Robertson et al. (2009) of protohalo overdensity in CDM simulations, and the same results might not apply in the case of WDM.

5. Peak-in-peak. Another possible source of error is that the ESP framework treats the peak-in-peak problem approximately, by introducing the effects of the peaks constraint as an extra weight in the mass function, rather than by explicitly accounting for spatial correlations between walks centred at different locations in space (see e.g. Scannapieco & Barkana 2002), and this approximation needs testing.

We address these issues in the next two sections by exploring the properties of the initial conditions of the simulation in greater detail.⁴

LAGRANGIAN PROPERTIES OF HALOES

In this section, we turn to the initial conditions of the simulation and perform an in-depth study of the Lagrangian properties of regions that will eventually form haloes; we call such regions protohaloes and give a precise definition below. Several authors have performed such studies in CDM simulations (e.g. Bond & Myers 1996; White 1996; Sheth et al. 2001; Porciani, Dekel & Hoffman 2002; Robertson et al. 2009; Ludlow & Porciani 2011a,b; Elia, Ludlow & Porciani 2012; Despali et al. 2013). To our knowledge, the current work is the first to extend these studies to the case of WDM, and is interesting for the reasons discussed in Section 2.

An advantage of using a WDM model with m_dm = 0.25 keV is that the number of objects is reasonably small. A disadvantage is that the half-mode mass is close to being unit significance, ν(M_hm, z = 0) ≃ 0.9, which does not allow us to explore low-significance objects with sufficient statistical precision. This could also potentially confuse non-linear assembly-bias-like effects with the peculiarities of halo formation at and below the half-mode mass scale. Nevertheless, this simulation provides us with an invaluable testing ground for several ideas in the peaks framework.

We focus on the overdensity of the protohalo patch (which is indicative of the collapse threshold), its curvature, velocity shear (ellipticity and prolateness) and moment of inertia. We will demonstrate two important features of the protohalo patches: (a) that they are all consistent with forming at initial density peaks and (b) their overdensities are strongly correlated with their ellipticities but not their prolateness.

Lagrangian density and shear fields

The initial conditions code music allows us to output the density field that was used to generate the simulation initial conditions as three-dimensional grid data. We used this function to re-generate the density field directly on a 512³ mesh with the same Fourier modes as the original simulation. We refer to this as the unsmoothed field. Using the particle IDs that encode the three-dimensional Lagrangian coordinate |$\boldsymbol {q}$| on the unperturbed initial particle lattice (i.e. before applying the Zel'dovich approximation), we can directly evaluate the density at |$\boldsymbol {q}$| without interpolating the perturbed particle position back on a grid. We linearly scale the density field to z = 0.

Protohaloes and their properties

For each halo, we recorded the particle IDs and recovered their respective Lagrange coordinates |$\boldsymbol {q}$|⁠. We call the set of N_k Lagrangian particles comprising halo k its Lagrangian patch or protohalo, denoted |$L_k\equiv \left\lbrace \boldsymbol {q}_i\,\left| \,i=1\cdots \ N_k \right. \right\rbrace$|⁠.

Fig. 3 shows three examples of protohaloes with masses ∼10¹³ h⁻¹ M_⊙. The top panel shows a well behaved protohalo. We notice two disconnected shells surrounding the connected interior of this patch. This is a beautiful example of the mapping between Lagrangian and Eulerian space. The gaps between the shells appear because the outer caustics of the halo are not inside the virial radius and are thus cut off. The two shells correspond to material on first and second infall. The other two examples show evidence of mixing due to large-scale interactions.

Figure 3.

Volume renderings of the (unsmoothed) density field in the protohalo patches for three haloes of mass 10¹³ h⁻¹ M_⊙. Each panel shows a light grey sphere centred on the protohalo centre, containing the same mass. The darker shaded grey volume indicates the actual protohalo patch. Coloured contours indicate isodensity regions at the values given underneath each image. Top panel: a good example of an ellipsoidal protohalo patch, with clearly visible disconnected outer regions. This protohalo was also assigned a matching ESPeak by the algorithm described in Section 4. Middle and bottom panels: two examples of protohaloes with evidence for substantial substructure and mixing due to large-scale interactions. Our algorithm did not find any matching ESPeak for these two objects.

Open in new tabDownload slide

Haloes form at peaks

In Fig. 4, we show the distribution of the total curvature x (equation 15; top panel) and the magnitude of the gradient |$\eta \equiv \sqrt{3}|{\nabla }\delta |/\sigma _1$| (bottom panel) averaged over the halo patches. For computing σ₁ and σ₂ which make these quantities dimensionless, we used a top-hat filter at the Lagrangian scale R_L of each object.

Figure 4.

(Top panel: the distribution of peak curvatures x (equation 15) at protohalo locations, averaged over the Lagrangian patch (blue histogram) and using a spherical aperture of the same mass (green histogram). The dashed curve shows the theoretically expected distribution (equation 36), while the dotted curve shows the distribution at random locations (a Gaussian with zero mean and unit variance). Bottom panel: the distribution of density gradient |$\eta \equiv \sqrt{3}|{\nabla }\delta |/\sigma _1$| at protohalo locations, averaged over the Lagrangian patch (blue histogram) and using a spherical aperture of the same mass (green histogram). The dotted line is the distribution at random locations; in this case η² is distributed as chi-squared with 3 degrees of freedom. These plots show that all haloes in our sample are consistent with having formed near initial density peaks.

Open in new tabDownload slide

The distribution of x for a Gaussian random field would be a Gaussian with zero mean and unit variance (dotted black curve in the top panel). The measured distribution on the other hand has only positive values as mentioned above, and its shape is very similar to the analytical prediction using ESP with a deterministic barrier (dashed black curve, see Section 5, equation 36), although the measured mean value for x is lower than the predicted mean by about 0.4.

The distribution of η for a Gaussian random field would be |$p(\eta )=\sqrt{2/\pi }\eta ^2{\rm e}^{-\eta ^2/2}$| (because in this case η² is chi-squared distributed with 3 degrees of freedom.) This is shown as the dotted black curve in the bottom panel; the measured values clearly populate the low tail of this distribution. (Ideally all the values would be zero.) We also see that the patch-averaged values of η have a significantly lower scatter than the spherically averaged ones. This is not surprising since the requirement η = 0 is quite unstable to choices of filtering, and the spherical filter is known to introduce an additional randomization as compared with the actual Lagrangian patch (BBKS; Despali et al. 2013).

We therefore conclude that all haloes in our sample are consistent with having formed near initial density peaks. In principle, we should also have explicitly checked for the presence of local density maxima at or near the protohalo locations, e.g. along the lines discussed by Ludlow & Porciani (2011a). This, however, would involve making a specific choice regarding the smoothing scale. We defer such a calculation to Section 4.1, where we implement an algorithm that makes this choice while simultaneously centring the smoothing filter at locations that are most likely to collapse according to the excursion set formalism.

Overdensity of Lagrangian patches

Fig. 5 shows the patch-averaged (left-hand panel) and spherically averaged (middle panel) overdensities of the protohaloes as a function of their mass. We find that the spherical overdensities are strongly correlated with the corresponding spherically averaged values of Y (equation 17). This is evident in the middle panel where we have coloured the points using 〈 Y 〉^(s). (The patch-averaged overdensities show a similar strong correlation with the patch-averaged Y; we omitted the colouring in the left-hand panel for clarity.)

Figure 5.

Patch overdensity as a function of protohalo mass. Left panel: the overdensity δ (equation 16) at protohalo locations, averaged over the Lagrangian patch (points). The thick blue line and the blue band show, respectively, the mean and standard deviation of the spherically averaged protohalo overdensities measured by Robertson et al. (2009) in CDM simulations, extrapolated to our WDM initial power spectrum (i.e. we used our WDM transfer function in the fit in their equation D2 with the Δ = 200 parameters from their table 1). Middle panel: protohalo overdensity δ averaged over spherical apertures of the same mass, coloured by the spherically averaged ellipticity Y = e_vδ (equation 17). We see a strong trend of δ with Y, as expected from ellipsoidal collapse arguments (Sheth et al. 2001). Right panel: the quantity δ − Y, coloured by the spherically averaged prolateness Z = p_vδ (equation 18). The residual scatter shows no trend with Z, which is at odds with the predictions of ellipsoidal collapse models.

Open in new tabDownload slide

Robertson et al. (2009) performed similar spherically averaged measurements in the initial conditions of the CDM simulations presented by Tinker et al. (2008). The left-hand panel of Fig. 5 shows the mean and standard deviation of the distribution of overdensities reported by Robertson et al. (2009), but using the WDM transfer function. We see that their spherically averaged measurements, extrapolated to WDM, are in reasonable agreement with our patch-averaged overdensities. Our spherical overdensities, on the other hand, have a higher mean and scatter than theirs (see also Elia et al. 2012, who found similar results in their CDM simulations). There is no clear reason for this discrepancy.

These results for the spherically averaged overdensity and shear ellipticity will form the basis of the empirical walks that we describe in the next section. In general, however, we note that our measurements of spherically averaged quantities tend to have larger scatter than the corresponding patch-averaged ones. For completeness, in Appendix C we also show the distributions of ellipticity and prolateness defined using the tidal tensor and the Hessian of the density.

In summary, the results of this section show us that haloes form at peaks and have Lagrangian (spherically averaged) overdensities that are consistent with equation (27), with a residual scatter that is uncorrelated with other properties such as shear prolateness or peak shapes. One aspect we have not explored here is the relative (mis-)alignment between the velocity shear, density Hessian and moment of inertia tensors, which can be an important ingredient in any recipe for predicting collapse time based on dynamical arguments. This is especially interesting given previous results from CDM simulations (Porciani et al. 2002; Despali et al. 2013) which suggest that, contrary to expectations based on Gaussian statistics (e.g. van de Weygaert & Bertschinger 1996), the direction of maximum initial compression is on average well aligned with the longest geometrical axis of the protohalo, rather than the shortest. We will return to an analysis of tensor alignments and their dynamical consequences in future work.

EMPIRICAL EXCURSION SET PEAK WALKS

To test this, in this section, we explicitly perform such random walks in the actual initial density field that was used for the numerical simulation and identify spherical peak patches which are predicted to form haloes using this barrier. We will refer to these objects as ESPeaks below.

Methodology

Our algorithm is essentially a more accurate version of what the analytical ESP calculation tries to achieve. We note that it is less sophisticated than the original peak-patch algorithm implemented by Bond & Myers (1996), since we are not interested in the final locations, profiles and velocity dispersions of the haloes, but only in their mass.

We consider a hierarchy of N_s = 100 smoothing scales R_i logarithmically spaced in the top-hat spherical mass contained in R_i between M₀ = 10¹⁰ h⁻¹ M_⊙ and |$M_{N_{\rm s}}=10^{15}\,h^{-1}\,{\rm M}_{{\odot }}$|⁠. We then proceed as follows, starting with the largest smoothing scale. A small fraction of objects have partially overlapping Lagrangian volumes. We flag the smaller of such pairs as ‘subpeaks’ and, for the current analysis, do not include them in the sample of protohaloes. (These form about 10 per cent of the total sample.) At the end, we arrive at a catalogue of ESPeaks whose Lagrangian properties we can analyse in exactly the same way as for the actual protohalo patches.

• We determine the coordinates |$\boldsymbol {x}_k$| of all peaks in |$\delta _{R_i}$|⁠. Being a peak requires that |$\delta _{R_i}$| is larger at |$\boldsymbol {x}_k$| than in all 26 surrounding cells.

• We discard all peaks for which |$\delta _{R_i}$| is below the barrier, i.e. where |$\delta _{R_i}(\boldsymbol {x}_k)< B(\boldsymbol {x}_k;R_{i})$|⁠; B is given by equation (28).

• Additionally, we discard all peaks that are within the Lagrangian radius of a peak that has been identified before. This explicitly solves the cloud-in-cloud problem.

• Finally, we also discard all peaks where the density was above threshold on a larger scale, i.e. where |$\delta _{R_{i+1}}(\boldsymbol {x}_k)>B(\boldsymbol {x}_k;R_{i})$|⁠. This step improves numerical stability but is otherwise redundant.

• We proceed to the next smaller scale i → i − 1 and start over at step 1.

The mass function of ESPeaks

The solid orange histogram in Fig. 1 corresponds to the mass function of ESPeaks obtained using the algorithm described above. The dashed blue histogram shows the result of the same algorithm, but now using a deterministic barrier equation (33) which, as we argue in the next section, is a useful approximation to equation (28). Indeed, we see that these two histograms agree quite well, indicating that stochasticity in the barrier arising from statistical fluctuations in the initial conditions does not lead to dramatic effects in the mass function.

The most important feature of the orange histograms is that they show a low-mass tail consistent with dn/d ln m ∝ m^2/3 as discussed earlier. In fact, the histograms are well described by the WDM version of the ESP calculation (solid black) of Paranjape et al. (2013) who used a stochastic barrier adjusted to match CDM simulations, as well as a similar ESP calculation with the deterministic barrier (equation 33) which we discuss below.

The overall number of ESPeaks identified by our algorithm (1261 when using equation 28 and 1335 when using equation 33) is reasonably close to the total number of protohaloes, which is 1522. The lower numbers of ESPeaks could partially be because we stop our algorithm at the lower mass limit of 10¹⁰ h⁻¹ M_⊙. For comparison, integrating the ESP prediction using equation (33) (dashed line in Fig. 1) above the free-streaming scale M_fs gives a prediction of 1380 objects in the simulation volume.

Matching ESPeaks and haloes

If the excursion set picture is valid, the ESPeaks we identify should be correlated with the actual protohaloes. Could the mass function mismatch simply be because our low-mass ESPeaks are not associated with protohaloes? To assess this, we match the protohalo catalogue and the ESPeaks catalogue as follows.

For every protohalo, we find the ESPeaks contained inside of a sphere of its Lagrangian radius, and associate the ESPeak of highest mass with the halo. This procedure matches 978 (i.e. 64 per cent) of the protohaloes to ESPeaks. We then repeat the same procedure matching ESPeaks to haloes using the filter radius on which the ESPeak was identified, and in this case we can match 935 (i.e. 74 per cent) of the ESPeaks to protohaloes. (Including the subpeaks in the analysis makes these numbers 68 and 72 per cent, respectively.) We discuss possible reasons for the relatively large fraction of mismatched objects below (see also Ludlow & Porciani 2011a). In Appendix C, we also discuss the effect of repeating the exercise with the ellipsoidal collapse barrier (equation 29).

As a visual example, we note that the well-behaved protohalo in the top panel in Fig. 3 was assigned a matching ESPeak while the other two distorted objects were not. While the majority of the protohaloes we can match to ESPeaks look like the object in the top panel and the majority of unmatched protohaloes are distorted, there are also a number of examples of well-behaved protohaloes that are not matched, as well as distorted protohaloes that are.

The top panel of Fig. 6 shows the masses of the protohaloes that we could match to ESPeaks compared to the corresponding ESPeak masses. The points are coloured by the spherically averaged protohalo sphericity S (equation 25). There is a strong trend of S with halo mass: low-mass haloes are decidedly aspherical. This is consistent with the CDM results of Ludlow & Porciani (2011b). Additionally, the scatter in mass assignment also correlates strongly with S, with low-mass, aspherical haloes having the largest scatter. However, at a given halo mass, the mass mismatch on average does not seem to correlate strongly with halo shape.

Figure 6.

Top panel: the masses M_halo of protohaloes that could be matched to ESPeaks, compared to the corresponding ESPeak mass M_walk, coloured by the spherically averaged protohalo sphericity S (equation 25). Low-mass protohaloes are clearly more aspherical than high-mass ones, and the scatter in mass assignment also increases significantly for the more aspherical objects. At a given halo mass, however, the average mass mismatch does not correlate strongly with S. Bottom panel: the orange (blue) histogram shows the mass function of protohaloes that could (not) be matched to ESPeaks. The grey histogram shows the total halo mass function (same as the red histogram in Fig. 1). The matched fraction is quite large at high masses and falls significantly at low masses.

Open in new tabDownload slide

The histograms in the bottom panel of the figure show the mass function of matched (orange) and unmatched (blue) protohaloes, with the grey histogram showing the total halo mass function (same as the red histogram in Fig. 1). The matched fraction is quite large at the highest masses (reaching 100 per cent for m ≳ 3 × 10¹³ h⁻¹ M_⊙), remains approximately constant at intermediate masses m ∼ M_hm and falls significantly at low masses m ≲ 10¹² h⁻¹ M_⊙.

In the top panel of Fig. 7, we show the ratio of ESPeak mass to protohalo mass as a function of ESPeak mass, for ESPeaks that we could match to haloes, coloured in this case by the spherically averaged protohalo curvature x. We see a weak trend of curvature with mass mismatch: ESPeaks matching shallower protohaloes appear to have larger mass mismatches.

Figure 7.

Top panel: the masses M_walks of ESPeaks that could be matched to protohaloes, compared to the corresponding protohalo mass M_halo, coloured by the spherically averaged protohalo curvature x (equation 15). We see a weak trend of the mass mismatch with curvature: shallower protohaloes tend to have larger mismatches between M_halo and M_walk. Bottom panel: the orange (blue) histogram shows the mass function of ESPeaks that could (not) be matched to protohaloes. The grey histogram shows the total halo mass function (same as the red histogram in Fig. 1). Most of the unmatched ESPeaks were assigned dramatically lower masses than any protohalo in the sample. The sum of the orange and blue histograms is consistent with analytical ESP predictions (compare the smooth black curves and solid orange histogram in Fig. 1).

Open in new tabDownload slide

The histograms in the bottom panel of the figure show the mass function of matched (orange) and unmatched (blue) ESPeaks, while the grey histogram is the same as in Fig. 6 and shows the mass function of all haloes. We clearly see that most of the unmatched ESPeaks were assigned dramatically lower masses than any protohalo in the sample. This could indicate that the ESP picture is, in fact, not appropriate for these objects. However, the presence of a small fraction of low-mass ESPeaks that do have matching protohaloes, which in turn have larger true masses and low curvatures, suggests that the explanation of this trend could be more subtle. We therefore explore the properties of the mismatched objects in more detail in the next subsection. Recall that the sum of the mass functions of matched and unmatched ESPeaks is consistent with analytical ESP predictions (compare the smooth black curves and solid orange histogram in Fig. 1).

Failures of the protohalo ↔ ESPeak matching

We have already seen above that the masses of the ESPeaks that cannot be matched to protohaloes are too low compared to the mass function of haloes. As illustrated in Fig. 8, these unmatched protohaloes tend to have significantly lower peak curvatures and higher densities (solid blue curves in the top and bottom panels, respectively) than those that could be matched to ESPeaks (solid orange curves in the respective panels). Similar trends are also seen, albeit to a lesser extent, in the densities and curvatures of ESPeaks that (do not) have matching protohaloes, as seen in the dashed orange (blue) curves in the two panels. The inset in the top panel of the figure shows the joint distribution of curvature and mass for the matched (orange) and mismatched (blue) protohaloes. The absence of a significant mass dependence of x indicates that the difference in peak curvature is not a result of the two populations occupying different mass ranges.

Figure 8.

Cumulative distributions of spherically averaged peak curvature (top) and density (bottom) for the protohaloes that could be matched to ESPeaks (solid orange) and for protohaloes for which no corresponding ESPeak could be found (solid blue). The dashed orange (blue) curves show the corresponding quantities for ESPeaks that could (not) be matched to protohaloes. The inset in the top panel shows the joint distribution of protohalo mass and curvature in the matched (orange) and unmatched (blue) cases. The absence of a significant mass dependence of x shows that the trend seen in the cumulative distribution – unmatched objects have smaller curvatures – is not entirely caused by the preferentially low masses of the unmatched objects.

Open in new tabDownload slide

Fig. 9 shows the distribution of mass and sphericity S for the matched (orange) and unmatched (blue) protohaloes. There is a strong trend of sphericity with mass, as noted in Fig. 6. Apart from this, however, there is no significant difference between the sphericities of matched and unmatched protohaloes. The dashed line, which is somewhat shallower than the measured mass trend, shows the linear fit to corresponding measurements in CDM simulations presented by Ludlow & Porciani (2011b).

Figure 9.

Joint distribution of protohalo mass and sphericity S (equation 25) in the matched (orange) and unmatched (blue) cases. There is a strong trend of sphericity with mass, as noted in Fig. 6. Apart from this, however, there is no significant difference between the sphericities of matched and unmatched protohaloes. The dashed line, which is somewhat shallower than the measured mass trend, shows the linear fit to corresponding measurements in CDM simulations presented by Ludlow & Porciani (2011b).

Open in new tabDownload slide

Discussion

Our comparison of the outcomes of simulation and the empirical walks suggests that, at least at masses m ≳ M_hm, the empirical approach more or less correctly predicts both the locations and masses of collapsed haloes. At smaller masses, however, the algorithm is able to predict the location of a collapsed object only for a relatively small fraction (∼55 per cent of protohaloes below 10¹² h⁻¹ M_⊙ have an associated ESPeak, and the fraction falls below 50 per cent quickly below this mass scale) and almost always predicts too small a mass in these cases. The large fraction of ESPeaks that cannot be matched to protohaloes are also predicted to have dramatically lower masses than any protohalo in the simulation.

We have also seen that, on average, the mass and position mismatches seem to be uncorrelated with the shapes of the protohaloes. That is to say, although smaller protohaloes are decidedly aspherical (with a scatter in mass mismatch that correlates with sphericity), there seems to be no trend between the average sphericity and the ratio m_ESPeak/m_halo in the matched cases (Fig. 6) and no significant difference between sphericities of matched and unmatched protohaloes with m_halo < M_hm (Fig. 9).

The unmatched protohaloes do have smaller curvatures than the matched ones. In principle, this could simply be because of their lower masses. However, the absence of a significant mass dependence of x in the inset in the top panel of Fig. 8 indicates that this is not the case. Additionally, in the case of ESPeaks matched to protohaloes, the protohalo curvature correlates with the mass mismatch (Fig. 7).

Note that the objects identified in the simulation of AHA13 were split into different types (see Section 2.2), the most important being ‘type-1’ (virialized haloes) and ‘type-2’ (objects in late stages of formation). The results above are for the 1522 ‘type-1’ objects, while we have ignored the 438 ‘type-2’ objects. If the latter also form from peaks, our choice of ‘type-1’ could be a cause for concern since a peaks-based analysis such as ESP might simply not be able to distinguish between them. Indeed, we do not find a significant difference between the protohalo regions of ‘type-2’ and ‘type-1’ objects. Moreover, we find that including the ‘type-2’ objects in the analysis on the same footing as ‘type-1’ leads to a larger number (1097) of matched ESPeaks, meaning that 87 per cent of our ESPeaks can be matched to some object that is either about to or has completely virialized. (The fraction of unmatched protohalo patches are now ∼44 per cent, as compared with 33 per cent when using only ‘type-1’. This is largely simply because the combined number of ‘type-1’ and ‘type-2’ objects (1960) is significantly larger than that of the ESPeaks, which can only be accommodated in the ESP calculation by lowering the collapse threshold.)

Interestingly, AHA13 found that the transition between ‘type-2’ and ‘type-1’ occurs fast and is associated with a rapid mass growth, bringing a ‘type-2’ object to a mass around or above the half-mode mass by the time it has virialized and thus turned into a ‘type-1’ object. This is consistent with the picture that power spectra with steeper (effective) slopes show enhanced accretion rates (Lacey & Cole 1993, 1994). These observations suggest that the excursion set calculation could be failing because it is unable to capture the quick mass growth that ‘type-1’ objects experience around the half-mode mass scale, possibly due to an incorrect prediction of collapse-time for a given peak patch. The rapid transition between these two types of objects means that even small errors in predicting the collapse time could dramatically alter the predicted locations and masses of fully virialized haloes. As we discuss in the next section, such an error can also account for the correlations we find between protohalo curvature and mass/location mismatches.⁵

ANALYTICAL RESULTS

In this section we use the results of our numerical study to motivate an analytical approximation which captures the sharp cut-off in the mass function better than the standard ESP calculation.

A possible explanation for the mismatch between ESPeaks and protohaloes

The behaviour discussed above might be explained if, at small masses, the algorithm systematically overpredicts the time at which a given peak patch should collapse. This is because a patch that collapses earlier than predicted will have time to accrete mass by the time of interest, and will consequently have a larger mass than predicted. As noted by AHA13, low-mass WDM haloes tend to grow much more rapidly than their high-mass counterparts, so even a small error in collapse time could have a dramatic impact on the predicted mass function. This is further corroborated by our observation in Section 4.3 that the assignment of peaks to either virialized haloes or objects in the late stages of formation is somewhat uncertain.

Let us suppose that there is in fact such a systematic uncertainty in collapse time. This is not an unreasonable assumption; similar effects have been noticed and discussed by other authors (Monaco 1999; Giocoli et al. 2007) in the case of CDM. Such effects could arise due to simplifying choices made in models such as ellipsoidal collapse (see Appendix B for a justification), as well as due to other physical mechanisms such as assembly bias.⁶ Can this explain the orders-of-magnitude mass increases that are required to go from the ESPeaks mass function to the halo mass function? To see why this is indeed the case, consider the following.

A collapse-time uncertainty can be interpreted as introducing a second barrier in the problem, in the sense that if a patch has been identified at mass m using the barrier B (in this case the one in equation 28) it must then be allowed to accrete until mass M when its density δ_patch reaches the ‘mass re-assignment’ barrier B_mra = B − ΔB, where ΔB is positive if the original collapse-time prediction was an overestimate. Fig. 10 illustrates the point. Since top-hat/Gaussian filtering induces strong correlations between walk heights, the scale M at which δ_patch reaches B_mra is essentially determined by the walk height and slope at scale m. The height at m is just B(m), and in the excursion set peaks picture the walk slope is strongly correlated with the peak curvature, so that dδ/dσ₀ ≈ x/γ (this relation is exact for Gaussian filtering). The model therefore says that sharper peaks will tend to have more closely matched masses, while shallower peaks will have larger mismatches.

Figure 10.

An illustration of mass reassignment: an ESPeak is identified at scale m with height δ = B_m and curvature x. If the collapse-time prediction in the dynamical model is a systematic overestimate, the ESPeak is allowed to accrete mass for a time Δt which translates into a lowering of the barrier by ΔB. Since the random walk describing the ESPeak has strongly correlated steps, the accretion occurs along essentially a straight line of slope |$\dot{\delta }=x/\gamma$| in the δ–σ₀ plane. This sets the mass scale M > m which is ultimately assigned to this ESPeak.

Open in new tabDownload slide

Note also that the WDM transfer function leads to a ‘freeze-out’ of all power spectrum integrals as m → 0 (cf. Section 2.3). This will amplify the above effect at small masses where a small change in σ₀ will imply a huge change in m. Additionally, peaks of lower significance will tend to be shallower on average, and this will also systematically enhance the mismatch at low masses.

If this idea is correct, then we should see two effects. First, matched protohalo patches with lower curvatures should have preferentially larger mass mismatches and vice versa; and secondly, unmatched protohaloes must have preferentially low values of curvature (a shallow walk that ‘freezes’ before crossing an incorrect barrier will not register as a potential halo). Fig. 7 is consistent with the first effect, while the second effect is seen quite clearly in Fig. 8 (see also the discussion in Section 4.5).

In the following, we will therefore assume that the idea of a collapse-time overprediction is correct, and leave for future work a more detailed modelling of the scatter of the mass mismatch by including e.g. the effect of protohalo sphericity. We can implement the notion of a second barrier in the ESP calculation as follows. We start by recapitulating the calculation of Paranjape et al. (2013), which we refer to as standard ESP.

Standard excursion set peaks

Re-assigning mass

An explicit example

The distribution of ΔB will in general depend on the scale m at which the ESPeak is originally identified; we know that at large m the mass assignment is essentially correct, with small scatter, while there is a trend towards underpredicting masses at small m. Since there is little guidance from theory for the actual values of |$\overline{\Delta B}$| and σ_ΔB, we have left these as free parameters, except for requiring that they become numerically small for large m or small σ_0m. The following is intended as a proof of principle, and we leave a more detailed analysis and estimate of |$\overline{\Delta B}$| to future work.

Figure 11.

Halo mass functions: the solid red histogram and dashed black curve are the same as in Fig. 1 and show, respectively, the mass function of haloes in the simulation and the ESP prediction (equation 34) using the deterministic barrier (equation 33). The solid blue curve shows the result of re-assigning masses using equation (41), with parameter values from equation (43). For comparison, the dot–dashed blue curve shows the fit based on sharp-k filtering presented by Schneider et al. (2013, their spherical collapse fit with q = 1 and c = 2.7), which is a good description of the mass function measured in their simulation. The dotted black curve shows the Tinker et al. (2008) fitting form.

Open in new tabDownload slide

The top panel of Fig. 12 shows our analytical mass function at z = 0 with masses re-assigned using the same parameter values as in equation (43) for three different values of dark matter particle mass m_dm (solid curves). We also show the corresponding curves fit by Schneider et al. (2013), again using their spherical collapse fit (dot–dashed curves). In the bottom panel we show our analytical prediction for m_dm = 0.25 keV, now for three different redshifts. The solid curves marked 0.25 keV in both panels are identical, and the same as the solid blue curve in Fig. 11.

Figure 12.

Top panel: analytical halo mass functions for different values of m_dm: 0.25 keV (red); 0.5 keV (yellow) and 0.75 keV (blue). The solid curves show the result of re-assigning masses in the ESP calculation using equation (41), with parameter values from equation (43). For comparison, the dot–dashed curves show the fit based on sharp-k filtering presented by Schneider et al. (2013, their spherical collapse fit with q = 1 and c = 2.7), which is a good description of the mass function measured in their simulation. Bottom panel: ESP mass functions with re-assigned masses at three different redshifts.

Open in new tabDownload slide

Consequences for CDM

The predictions of the ellipsoidal collapse model, augmented by a systematic uncertainty in mass assignment, accurately describe the WDM mass function. By the logic discussed earlier, the same expressions with the CDM transfer function should describe the CDM mass function. In particular, the half-mode mass scale for CDM is small enough that, in practice, every halo has a virialized progenitor. This means that the effects of collapse-time uncertainty – which were very pronounced around the half-mode mass of WDM due to the rapid growth of those objects – are now essentially an uncertainty in the time of major mergers, and consequently the associated mass mismatch must be significantly smaller. We see in Fig. 13 that this is indeed the case for masses m ≳ 10¹¹ h⁻¹ M_⊙. It is also reassuring to note that changing the value of |$\overline{\Delta B}$| in the CDM case has much less effect on the mass function at m ≳ 10¹¹ h⁻¹ M_⊙ than in the WDM case. For example, we have checked that increasing the amplitude in equation (43) by a factor of 1.5 or changing the exponent from 3 to 2 both lead to ≲3 per cent changes in the CDM mass function for m ≳ 10¹¹ h⁻¹ M_⊙.

Figure 13.

ESP mass function predictions for CDM using the deterministic barrier (equation 33) (dashed black) and after re-assigning masses (solid blue). The dotted black curve shows the Tinker et al. (2008) fitting form, which was calibrated over a range between ∼10^10.5 and ∼10^15.5 h⁻¹ M_⊙. The solid red curve shows the prediction, after mass re-assignment, for a ‘realistic’ WDM particle candidate with m_dm = 3.3 keV. In all cases the mass re-assignment parameters are the same as in equation (43). The lower panel shows the ratio of the CDM predictions to the Tinker et al. fit.

Open in new tabDownload slide

At lower masses our specific implementation of mass re-assignment predicts a factor of ∼2 larger number of haloes than expected from the mass function fit by Tinker et al. (2008). This behaviour of the re-assigned mass function at low masses is not very robust; however, it is sensitive to the specific numerical choice in equation (42). It is possible to adjust this number (and those in equation 43) to simultaneously get a good match to the CDM and WDM simulations, although we have not pursued this exercise here. One must also keep in mind that the Tinker et al. (2008) fit was calibrated for masses between ∼10^10.5 and ∼10^15.5 h⁻¹ M_⊙.

In other words, our proposed modifications to the ESP calculation not only correctly describe the sharp turn in the WDM mass function, but also describe the CDM mass function with the same accuracy as the standard ESP calculation of Paranjape et al. (2013). For CDM we have also checked that the linear Lagrangian halo bias predicted by this model (not shown) matches measurements in CDM simulations (Tinker et al. 2010) with the same accuracy as the Paranjape et al. (2013) calculation. However, as noted earlier, low mass in WDM does not mean low significance, and, in principle, low significance CDM haloes could be different from low-mass WDM haloes. Testing this would need high-resolution CDM simulations, or WDM simulations with a slightly larger mass such as m_dm ∼ 0.5 keV.

DISCUSSION AND CONCLUSIONS

Can we predict where and when haloes form? In this paper, we have thoroughly evaluated our ability to predict the abundance of collapsed objects by performing an in-depth analysis of the properties of the initial density field at the locations where collapse occurs in numerical simulations. To accomplish this, we used a perturbation spectrum with a small-scale cut-off such as those arising in WDM cosmologies. As discussed in Section 2, the resulting suppression of low-mass haloes provides powerful additional leverage which is absent in the CDM case.

Numerical simulations have traditionally had great difficulty in making a prediction for the abundance of haloes in such a scenario due to the artificial fragmentation of filaments – a problem that has only recently been overcome by AHA13. As a consequence, we were in the unique situation of being able to perform a thorough comparison between this numerical experiment and the mass function predicted from excursion set theory, by analysing the properties of haloes on an object-by-object basis. We summarize our results below, and discuss some outstanding issues.

It is well known that the standard excursion set approach predicts a mass function that is completely inconsistent with the numerical results (see e.g. Schneider et al. 2012, and also our Fig. 1). We showed that the inclusion of the peaks constraint in excursion sets (Musso & Sheth 2012; Paranjape & Sheth 2012; Paranjape et al. 2013) leads to a turnover in the mass function as well as an overall number of collapsed objects that is consistent with the simulation results. However, it also predicts masses around and below the half-mode mass scale that are significantly smaller than those measured in the simulation, leading to a small-mass slope of the mass function (dn/d ln m ∼ m^2/3) that is inconsistent with that found from the simulation (cf. Fig. 1). This prediction is remarkably robust against changing details of the calculation such as the shape and stochasticity of the barrier.

We next investigated the origin of this discrepancy between simulation and theoretical predictions. In particular, we analysed the Lagrangian properties of ‘protohaloes’ (the initial locations of groups of particles which are eventually identified as haloes in the simulation), and also performed empirical excursion set peak walks in the initial density field used in the simulation. We can summarize our findings as follows. Based on these results, we argued that the likely cause for the observed mass mismatch between ESPeaks and protohaloes is a systematic overprediction of the collapse time for a given perturbation. We then showed how such an uncertainty can be accounted for and corrected in the excursion set language (Section 5.1 and Fig. 10), and presented an explicit example of such a correction which describes the numerical WDM results very well (Figs 11 and 12). As an important consistency check, we also showed that the same model gives an accurate description of the CDM mass function, with the simple replacement of the WDM initial power spectrum with that of CDM (Fig. 13).

1. All haloes in the simulation are consistent with forming near peaks in the initial density field (Fig. 4).

2. The overdensities of protohaloes are strongly correlated with their shear ellipticities, but show no correlation with the shear prolateness (Fig. 5). The former is expected from arguments based on ellipsoidal collapse dynamics, while the latter is not (compare equation 28 with equation 29). The fact that the protohalo overdensity has no correlation with its prolateness is intimately connected with the distribution of individual shear eigenvalues and hence with the dynamical ordering of the collapse times of each (Ludlow & Porciani 2011b; Despali et al. 2013). It will be interesting to find a dynamical model that is consistent with our results, perhaps along the lines presented by Ludlow & Porciani (2011b).

3. The number of ‘ESPeaks’ (1261) identified by our empirical algorithm (Section 4.1) is reasonably close to the actual number of protohaloes (1522).

4. A significant fraction (74 per cent) of ESPeaks can be matched to actual protohaloes, while 64 per cent of the protohaloes can be matched to ESPeaks (details in Section 4.3).

5. The curvatures of these matched objects are significantly higher, and their overdensities significantly lower, than those of protohaloes and ESPeaks that could not be matched to each other (Fig. 8).

6. Most strikingly, the masses of ESPeaks are systematically lower than the protohalo masses. This is true for both matched and unmatched objects (respectively, top and bottom panels of Fig. 7). Since the ESPeak mass function is very well described by the ESP calculation (Fig. 1), this fully accounts for the discrepancy between the ESP halo mass function and that measured in the simulation.

7. For matched objects, the mismatch in mass assignment correlates with protohalo curvature (Fig. 7), while the scatter in the mass mismatch correlates with protohalo shape (Fig. 6).

8. We have checked that, apart from having larger overdensities and lower curvatures than their matched counterparts, the unmatched protohaloes do not appear to be special in any other property related to the velocity shear, density Hessian or moment of inertia.

9. If we also include in the analysis objects in late stages of formation (‘type-2’ in AHA13), 87 per cent of the ESPeaks can be matched to protohaloes (although the fraction of unmatched protohaloes are now larger, largely because the total number of protohaloes increases). These ‘type-2’ objects are known to be undergoing rapid mass growth (AHA13).

We emphasize that our solution works because it explicitly alters the mass assignment step of the ESP calculation, in our case by introducing a second barrier. Simply introducing new statistical variables defined by smoothing the initial density field in a single-barrier calculation (say, by setting B → δ_c + Y + xσ₀) would not work, because the predicted mass function in this case would still behave as dn/d ln m ∼ m^2/3 at low masses, as discussed in Section 2.3. At the heart of this issue is the difference between the physics of individual halo formation and the statistics of the initial density field: mismatches in collapse time predictions are primarily a physical, not statistical, problem. In CDM, since the σ₀(m) relation is always steep, errors in the physical collapse model can be accommodated by altering the statistical modelling (e.g. by changing the barrier shape as a function of σ₀, or by introducing stochasticity in the barrier). WDM, on the other hand, presents us with a situation where such solutions no longer work since the σ₀(m) relation ‘freezes out’ at small masses (Fig. 2). A full solution of the problem would likely involve a single barrier with an explicit dependence on mass (rather than σ₀) which is fixed by an accurate model of collapse.

Our work can be extended in several directions which could yield clues towards building a more predictive model. Although we argued that uncertainties in collapse time can easily arise in toy models such as ellipsoidal collapse, or due to assembly-bias-like effects, we have not provided conclusive evidence pinpointing a specific physical mechanism. Furthermore, the protohalo curvatures we measure are significantly lower than the simplest ESP prediction (Fig. 4). While this might be partially accounted for by the mass re-assignment, it is not clear whether there is a deeper reason for this discrepancy. Also, the scatter in ESPeak mass at fixed halo mass for matched objects correlates strongly with the protohalo shapes (Fig. 6). Our mass re-assignment currently incorporates only curvature, and it will be interesting to additionally account for protohalo shapes. The overdensities of unmatched protohaloes have a significantly higher tail than that of matched objects (Fig. 8). This could be indicating that low-mass WDM haloes can form near density peaks without necessarily satisfying the excursion set peaks constraints; e.g. there could be some other criterion for a sufficiently overdense patch to virialize. One way forward would be to analyse the local environments of these unmatched objects to look for peculiarities. Finally, while our chosen particle mass of m_dm = 0.25 keV is untenably warm, any realistic massive dark matter candidate would lead to very similar phenomenology around its half-mode mass scale, and our analysis would be relevant whenever this scale can be numerically resolved. As mentioned earlier, however, for a particle much colder than our present choice, low-mass haloes will also have low significance (ν ≪ 1) which was not possible in our case, and this could lead to additional effects. Clearly, it will be of great interest to investigate these aspects further in future work.

We thank Raul Angulo for providing us with the halo catalogues published in AHA13. We wish to thank Ravi Sheth, Tom Abel and Cristiano Porciani for insightful discussions and comments on the draft. OH acknowledges support from the Swiss National Science Foundation (SNSF) through the Ambizione fellowship.

REFERENCES

APPENDIX A: FULLY STOCHASTIC EXCURSION SET PEAKS

Here we present some of the formal details of the excursion set peaks calculation. This will highlight the difficulty in working with the full stochasticity that must be dealt with when using a barrier such as equation (28), and motivate the simpler, deterministic approximation (equation 33) used in the main text.

Formal expression for mass function

In what follows, an overdot denotes a derivative with respect to σ₀. All variables are assumed to be Gaussian filtered on a scale R which is related to the mass m through |$m=(4\pi /3)\bar{\rho }R_{\rm TH}^3$|⁠, where R_TH satisfies 〈 δ_G(R)δ_TH(R_TH) 〉 = 〈 δ_TH(R_TH)² 〉 with the subscripts ‘G’ and ‘TH’ on δ denoting Gaussian and top-hat filtering, respectively (Paranjape et al. 2013, see also Footnote 7). In practice this gives R ≈ 0.46R_TH with a slow variation. We will drop the filter subscripts below.

The raw number density of peaks |$\mathcal {N}_{\rm pk}( {\boldsymbol X})$| defined in equation (A2) consists of the following quantities: the joint (Gaussian) distribution function p(∂_ijψ, ∂_ijδ, ∇δ) of the shear, density gradient and density Hessian smoothed on scale R; the peaks constraint Pk(∂_ijδ, ∇δ) which enforces ∇δ = 0 and ζ_i < 0, where ζ_i are the eigenvalues of ∂_ijδ; and the excursion set constraint ES(σ₀, {δ, B}) which enforces up-crossing⁸ of the barrier by the random walk at the scale σ₀(R). The compact notation in ES(⋅⋅⋅) hides the fact that the up-crossing condition will introduce a dependence on |$\dot{\delta }=x/\gamma$| and possibly other stochastic quantities through |$\dot{B}$|⁠. The Dirac delta in equation (30) then sets the mass to be |$m=\bar{m}(\sigma _0)$|⁠, where |$\bar{m}(\sigma _0)$| is the inverse function of σ₀(R(m)) as discussed above, which is straightforward to compute numerically.

Consequences of misalignment

The 12 degrees of freedom in the two tensors ∂_ijψ and ∂_ijδ are most conveniently organized as the three eigenvalues of each, the three relative Euler angles between their respective eigenvectors and three additional Euler angles that fix the orientation of one of them with respect to the chosen basis.

The correlation structure between the tensors ∂_ijψ and ∂_ijδ implies that, in a generic realization of the field, they will be misaligned. That is, although they are strongly correlated (with correlation coefficient γ), their respective eigenvectors will not be parallel (see e.g. Desjacques 2008; Lavaux & Wandelt 2010; Despali et al. 2013). The misalignment is captured by the three relative Euler angles, and these must be marginalized over (the other three angles never appear in the distribution due to statistical isotropy, and can be trivially marginalized over). The effect of this marginalization is to introduce a non-trivial coupling between the anisotropic combinations of eigenvalues of the two tensors. So if we define y and z as in equations (19) and (20) then, as shown in an elegant calculation by Desjacques (2008), the marginalization over Euler angles will introduce terms involving products between at least one of {(y + z), (3y − z)} and at least one of {(Y + Z)/σ₀, (3Y − Z)/σ₀}, where Y and Z were defined in equations (17) and (18), respectively. This is quite different from the coupling between the isotropic variables x and μ through the familiar Gaussian term |${\sim }{\rm e}^{-(x-\gamma \mu )^2/2(1-\gamma ^2)}$| that appears in the BBKS calculation.

This anisotropic coupling is irrelevant if we were to ignore the excursion set constraint ES(⋅⋅⋅) in equation (A2) or if the barrier B in ES(⋅⋅⋅) did not depend on any of the anisotropic variables {y, z, Y, Z} (an example would be the barrier used by Paranjape et al. 2013), since in these cases one could marginalize over, e.g. Y and Z to recover expressions similar to those analysed by BBKS (note that ignoring the excursion set constraint altogether would give back the BBKS calculation exactly). In the case of the barrier (equation 28), however, the calculation will involve integration of a term like |${\sim } p(\lambda |\zeta ){\rm e}^{-(\mu +Y/\sigma _0-\gamma x)^2/2(1-\gamma ^2)}$| over Y and Z, where the distribution p(λ|ζ) of the eigenvalues of ∂_ijψ given the eigenvalues of ∂_ijδ is given by equation (22) of Desjacques (2008).

Figure A1.

Accounting for the effects of stochasticity: ESP mass function predictions for WDM using the stochastic barrier (equation 28) and the approximation |$\dot{Y} = y/\gamma$| (see text for details). From bottom to top, the solid curves show the result of keeping the zeroth-, second- and fourth-order terms in equation (A9), while the dashed curves (almost indistinguishable from the upper solid curves) include the third- and fifth-order terms. The upper short-dashed black curve is the same as in Fig. 1 and shows the ESP calculation using the deterministic barrier (equation 33), which gives a very good description of the ESPeaks mass function (solid orange histogram in Fig. 1) obtained with the stochastic barrier (equation 28). While the odd order terms do not contribute significantly to the stochastic calculation, the even order terms have clearly not converged – the fourth-order term gives a significant contribution compared to the second-order calculation.

Open in new tabDownload slide

A simpler way out

As the results of the previous section show, it rapidly becomes very complicated to deal with the full stochasticity inherent in even a simple barrier prescription like equation (28). Luckily, approximating this barrier with the deterministic expression (33) leads to an excellent description of the mass function of ESPeaks found using the stochastic barrier (equation 28) (compare the dashed black curve in Fig. 1 with the solid orange histogram).

The main message here is that the barrier stochasticity induced by ellipsoidal effects is a technical detail that does not address the main problem – that of mass re-assignment – we are faced with. Rather, the simplicity of the deterministic barrier solution allows us to tackle this problem in a computationally straightforward way, as discussed in the main text. In principle, one could imagine having a physically better motivated model of mass re-assignment; provided it can be expressed in the formal language of equation (37), the hurdles in using it to make mass function predictions would be purely technical.

APPENDIX B: ELLIPSOIDAL DYNAMICS AND COLLAPSE-TIME UNCERTAINTY

To understand why one might expect a small systematic uncertainty in the predicted collapse time in ellipsoidal dynamics, it will help to first recapitulate some of the basic features of spherical collapse.

The key conceptual point here is that the model itself reaches the virial radius at t_vir = t(θ = 3π/2) which occurs slightly before the collapse time t_coll = t(θ = 2π); in particular, t_vir ≈ 0.91t_coll. However, one can argue that virialization actually occurs over a dynamical time-scale t_dyn (the free-fall time at virialization; see e.g. Padmanabhan 1993) which, coincidentally, is of order 10 per cent of the collapse time. So in practice it is quite reasonable to use t_coll(≈t_vir + t_dyn) in place of t_vir when computing the linearly extrapolated overdensity, which gives the well-known result δ_L = 1.686. (Had we instead evaluated δ_L at θ = θ_vir, we would get δ_L = 1.583, a ∼6 per cent decrement from the traditional value; see e.g. Bond & Myers 1996, hereafter BM96.)

In the ellipsoidal model discussed by BM96 (see also White & Silk 1979; Eisenstein & Loeb 1995; Monaco 1999), it is more difficult to identify conserved quantities since there is no spherical symmetry and non-linear tides can have a complicated influence. The definition of virialization is therefore somewhat ambiguous. BM96 settled on a simple prescription in which each principle axis is evolved non-linearly until it shrinks to a predetermined fraction f_c of the global scale factor, after which it is assumed to remain fixed at that value. Collapse is defined as the time at which the longest axis (i.e. the smallest shear eigenvalue) satisfies this condition. Choosing f_c = 0.178 ensures that the overdensity at virialization for a spherical configuration is |$f_{\rm c}^{-1/3}\simeq 178$|⁠, the traditional value. Since there is no compensation for dynamical time-scales, it seems reasonable that this prescription underpredicts collapse times. In the spherical limit, this would be a ∼10 per cent effect as discussed above.

The fit (equation 29) by Sheth et al. (2001, hereafter SMT01) to the resulting linear overdensity values at collapse as a function of ellipticity and prolateness reduces to the traditional spherical collapse result when Z = 0 = Y, meaning that it rescales the BM96 prescription for the barrier by a factor of 1.686/1.583 = 1.06. This would cause a corresponding ∼10 per cent increase in collapse time. While this increase is of the correct magnitude in the spherical limit (see above), it is less clear whether this is also true for significantly triaxial configurations. In the EdS case where the growth factor is proportional to the scale factor, one can check that a barrier rescaling as above does, in fact, rescale the collapse time by a constant at any ellipticity. However, the situation is different e.g. in flat ΛCDM with Ω_m < 1. In this case, the growth factor is different from the scale factor and one can show that a barrier rescaling leads to a slightly larger effect on the collapse time for e_v > 0 than it does for e_v = 0. In other words, if one were aiming for a ∼10 per cent increase in collapse time for all ellipticities, then rescaling the barrier by a constant will tend to overcompensate at large ellipticities, and would go in the direction of explaining the results of Section 5.3.

Of course, the specific prescription for virialization itself is somewhat ad hoc, and BM96 discuss several modifications, all of which lead to few per cent changes in collapse time (see also Angrick & Bartelmann 2010; Ludlow & Porciani 2011b). Additionally, the specific choice of halo finder in the simulation will also lead to (probably unquantifiable) systematics. It is therefore not unreasonable that uncertainties in dynamical modelling lead to a small (but, in the WDM case, important) systematic overprediction of collapse times.

APPENDIX C: PEAK SHAPES AND THE COLLAPSE THRESHOLD

Fig. C1 shows the distributions of ellipticity and prolateness defined using the tidal tensor (top panel) and using the Hessian of the density (bottom panel), coloured by the halo mass in each case. These quantities were defined in equations (17)–(20) above.

Figure C1.

Top panel: the distribution of patch-averaged ellipticity e_v (equation 17) and prolateness p_v (equation 18) of the tidal tensor ∂_ijψ, coloured by halo mass. Both the median value and scatter of e_v increase with decreasing halo mass, whereas we find no trend in p_v with halo mass. Bottom panel: corresponding patch-averaged quantities e_pk (equation 19) and p_pk (equation 20) defined for the density Hessian ∂_ijδ. (Note that the scale on the axes is different from the top panel.) There is a weak preference for low-mass objects to have p_pk > 0, which is consistent with the BBKS results for peak shapes.

Open in new tabDownload slide

In the top panel we see that low-mass haloes have a weak preference for larger values of ellipticity e_v, while there is no trend of prolateness p_v with mass. This is consistent with the results of Ludlow & Porciani (2011a) and Despali et al. (2013) for CDM haloes. In the bottom panel (note the difference in axes scales) we see that there is a large scatter in values of e_pk at any mass, while there is a weak trend for p_pk > 0 at low masses. The latter is consistent with the BBKS results for peak shapes.

We noted earlier that the barrier (equation 29) associated with ellipsoidal dynamics does not describe the measured densities of the protohalo patches; in particular, we found no correlation between the measured protohalo density and prolateness. We have repeated the exercise of finding ESPeaks and matching them with protohaloes using the full ellipsoidal collapse barrier (equation 29). While we find that a similar fraction (∼62 per cent) of protohaloes have matching ESPeaks, the scatter in the assigned masses, especially at low masses, is somewhat larger in this case as compared to using equation (28). Fig. C2 compares the mass distributions of matched objects in these two cases.

Figure C2.

Comparing ESPeak masses assigned by the algorithm described in Section 4.1 when using the simple stochastic barrier (equation 28; blue points) and when using the fully stochastic SMT01 barrier (equation 29; orange points). The points show the masses M_halo of protohaloes that could be matched to ESPeaks, against ESPeak mass M_walks. In both cases the matched fraction was comparable (64 per cent when using equation 28 and 62 per cent when using equation 29); however, the scatter when using the SMT01 barrier is somewhat larger than when using equation (28), especially at low masses, which can be traced to the fact that the SMT01 barrier introduces a dependence of the peak-centred overdensity on prolateness which is not present in the case of the actual protohaloes.

Open in new tabDownload slide