Excursion set peaks in energy as a model for haloes

TL;DR

This paper introduces a new model for dark matter halo formation based on energy peaks and gravitational flow convergence, offering a more physically motivated and divergence-free approach that aligns with simulation results.

Contribution

It replaces the traditional density maxima assumption with energy minima and gravitational flow convergence, providing a more physically grounded and divergence-free model for halo formation.

Findings

Predicts scatter in protohalo overdensities consistent with simulations

Removes divergences present in standard models for cosmological power spectra

Shows that energy-based peaks can serve as a robust alternative to density maxima

Abstract

The simplest models of dark matter halo formation are based on the heuristic assumption, motivated by spherical collapse, that virialized haloes originate from initial regions that are maxima of the smoothed density field. Here, we replace this notion with the dynamical requirement that protohalo patches be regions where the local gravitational flow converges to a point. For this purpose, we look for spheres whose gravitational acceleration at the boundary -- relative to their center of mass -- points towards their geometric center: that is, spheres with null dipole moment. We show that these configurations are minima of the total energy, i.e. the most energetically bound spheres. For this reason, we study peaks of the energy overdensity field, and argue that the approach shows considerable promise. This change simply requires that one modify the standard top-hat filter, with the added important benefit that, for power spectra of cosmological interest, the resulting model is no longer plagued by divergences. Although the formalism is no more complicated than the overdensity based approach, the model is richer in the sense that it naturally predicts scatter in the overdensities of protohalo patches that are destined to form haloes of the same mass, in qualitative agreement with simulations of halo formation.

Equations

$${\delta }_R\equiv {\int }_{V_R}\frac{\mathrm{d}\mathbf{r}}{V_R}\delta (\mathbf{r}),$$

$$U={\int }_{V_R}\mathrm{d}\mathbf{r}\rho (\mathbf{r})(\mathbf{r}-{\mathbf{r}}_{\mathrm{cm}})\cdot \mathbf{g},$$

$${\mathbf{r}}_{\mathrm{cm}}\equiv \frac{1}{M}{\int }_{V_R}\mathrm{d}\mathbf{r}\rho (\mathbf{r})\mathbf{r}=\frac{\bar{\rho }}{M}{\int }_{V_R}\mathrm{d}\mathbf{r}\delta (\mathbf{r})\mathbf{r},$$

$$\bm{g}=-4\pi G\bar{\rho}\bigg{[}\frac{\bm{r}-\bm{r}_{\mathrm{cm}}}{3}+\nabla\phi(\bm{r})-\int_{V_{R}}\frac{\mathrm{d}\bm{r}^{\prime}\rho(\bm{r}^{\prime})}{M}\nabla\phi(\bm{r}^{\prime})\bigg{]},$$

$$U=-4\pi G\bar{\rho }\frac{M{R}_I^2}{5}(1+{ϵ}_R),$$

$$R_I^2\equiv \frac{5}{3M}{\int }_{V_R}\mathrm{d}\mathbf{r}\rho (\mathbf{r})|\mathbf{r}-{\mathbf{r}}_{\mathrm{cm}}{|}^2\equiv \frac{5}{3M}T$$

$${ϵ}_R=\frac{5}{M{R}_I^2}{\int }_{V_R}\mathrm{d}\mathbf{r}\rho (\mathbf{r})(\mathbf{r}-{\mathbf{r}}_{\mathrm{cm}})\cdot \nabla \varphi (\mathbf{r})$$

$${ϵ}_R\equiv \frac{15}{4\pi R^5}{\int }_{V_R}\mathrm{d}\mathbf{r}\mathbf{r}\cdot \nabla \varphi (\mathbf{r}).$$

$${ϵ}_R=\frac{5}{R^5}{\int }_0^R\mathrm{d}r{r}^4{\delta }_r$$

$${\delta }_R(\mathbf{x})=\int \frac{\mathrm{d}\mathbf{k}}{(2\pi {)}^3}\delta (\mathbf{k})W_1(kR){\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{x}},$$

$${ϵ}_R(\mathbf{x})=\int \frac{\mathrm{d}\mathbf{k}}{(2\pi {)}^3}\delta (\mathbf{k})W_2(kR){\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{x}},$$

$${\mathbf{D}}_R=\frac{3}{R}\int \frac{\mathrm{d}\mathbf{k}}{(2\pi {)}^3}\frac{i\mathbf{k}}{k^2}\delta (\mathbf{k})j_2(kR).$$

$$g_{i}=-\frac{GM}{R^{2}}\bigg{[}\hat{r}_{i}+D_{R,j}(3\hat{r}_{i}\hat{r}_{j}-\delta_{ij})+\dots\bigg{]}$$

$${\mathbf{D}}_R=\frac{R}{5}\int \frac{\mathrm{d}\mathbf{k}}{(2\pi {)}^3}i\mathbf{k}\delta (\mathbf{k})W_2(kR)=\frac{R}{5}\nabla {ϵ}_R.$$

$${\nabla }_i{\nabla }_j{ϵ}_R=\frac{5}{R}\int \frac{\mathrm{d}\mathbf{k}}{(2\pi {)}^3}\frac{k_i{k}_j}{k^2}\delta (\mathbf{k})\frac{\mathrm{d}{W}_1(kR)}{\mathrm{d}R},$$

$$E\equiv\frac{5}{3M}\bigg{(}\!K+U-\frac{\Lambda}{6}T\bigg{)},$$

$$E=\frac{{\dot{R}}_I^2}{2}+k{R}_I^2-\frac{G{M}_I}{R_I}-\frac{\Lambda }{6}{R}_I^2,$$

$$E\simeq -\frac{G{M}_I}{R_I}{ϵ}_R.$$

$${\ddot{R}}_I=2k{R}_I-\frac{G{M}_I}{R_I^2}+\frac{\Lambda }{3}{R}_I$$

$$W_n(z)\equiv (2n+1)!!\frac{j_n(z)}{z^n},$$

$$\mathrm{d}{W}_n(z)/\mathrm{d}\ln z=-z^2{W}_{n+1}(z)/(2n+3),$$

$$W_n(z)-W_{n-1}(z)=z^2{W}_{n+1}(z)/[(2n+3)(2n+1)],$$

$${\sigma }_{jn}^2(R)\equiv \int \frac{\mathrm{d}k}{k}{k}^{2j}\frac{k^{3}P(k)}{2{\pi }^2}{W}_n^2(kR).$$

$${\mathcal{W}}_n(r)=\{\begin{array}{ll}\frac{\Gamma (n+3/2)}{\Gamma (n)\Gamma (1/2)}\frac{[1-(r/R{)}^2{]}^{n-1}}{\pi R^3}& \text{mbox}ifr\le R\\ 0& \text{mbox}otherwise\end{array}.$$

$${\nu }_R(\mathbf{x})\equiv \frac{{\delta }_R(\mathbf{x})}{{\sigma }_{01}(R)}\mathrm{and}{\omega }_R(\mathbf{x})\equiv \frac{{ϵ}_R(\mathbf{x})}{{\sigma }_{02}(R)}.$$

$${\eta }_{i,R}\equiv \frac{{\nabla }_i{ϵ}_R}{{\sigma }_{12}(R)}=\int \frac{\mathrm{d}\mathbf{k}}{(2\pi {)}^3}\frac{i{k}_i\delta (\mathbf{k})}{{\sigma }_{12}(R)}{W}_2(kR){\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{x}}$$

$${\zeta }_{ij,R}\equiv -\frac{{\nabla }_i{\nabla }_j{ϵ}_R}{{\sigma }_{22}(R)}=\int \frac{\mathrm{d}\mathbf{k}}{(2\pi {)}^3}\frac{k_i{k}_j\delta (\mathbf{k})}{{\sigma }_{22}(R)}{W}_2(kR){\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{x}},$$

$$x_R\equiv -\frac{{\nabla }^2{ϵ}_R}{{\sigma }_{22}(R)}=\int \frac{\mathrm{d}\mathbf{k}}{(2\pi {)}^3}\frac{k^2\delta (\mathbf{k})}{{\sigma }_{22}(R)}{W}_2(kR){\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{x}},$$

$$y_R\equiv -\frac{7\mathrm{d}{ϵ}_R/\mathrm{d}\ln R}{R^2{\sigma }_{23}(R)}=\int \frac{\mathrm{d}\mathbf{k}}{(2\pi {)}^3}\frac{k^2\delta (\mathbf{k})}{{\sigma }_{23}(R)}{W}_3(kR){\mathrm{e}}^{i\mathbf{k}\cdot \mathbf{x}}.$$

$$x_R=-\frac{5\mathrm{d}{\delta }_R/\mathrm{d}\ln R}{R^2{\sigma }_{22}(R)},$$

$$y_R=35\frac{{ϵ}_R-{\delta }_R}{R^2{\sigma }_{23}(R)}.$$

Full text

††pubyear: 2021††pagerange: Excursion set peaks in energy as a model for haloes–C

1 Introduction

Virialized dark matter haloes play an important role in models of nonlinear strucure formation (Cooray & Sheth, 2002). The most basic quantity of interest is the comoving number density of haloes as a function of halo mass. Following Press & Schechter (1974), this halo abundance retains memory, i.e. depends on the statistics, of the primordial fluctuation field. The simplest models, inspired by spherical collapse, assume that the initial overdensity within a spherical region must equal (exceed) a critical value if it is to become a virialized halo today (before today). Energy conservation arguments are used to estimate its final density (Gunn & Gott, 1972).

This assumption has led to extensive study of the abundance of such critically overdense patches in the suitably smoothed initial overdensity fluctuation field. To avoid double-counting patches which are overlapping with or contained within larger ones, there are two key approximations: one is that such patches must be local peaks of the smoothed field (Bardeen et al., 1986), and the other is that they must not be as overdense on larger scales (Bond et al., 1991). These requirements can be combined (Bond, 1989; Appel & Jones, 1990; Manrique & Salvador-Solé, 1995), and together they define the excursion set peaks approach of Paranjape & Sheth (2012) and Paranjape et al. (2013), which provides a good description of halo statistics (mass function and bias) once some input from simulations is included to tune the value of the critical overdensity.

However, some problems with the peak based approach remain open. First, when investigating the predictions of this model on a halo by halo basis, at small masses not all protohaloes are peaks of the initial density field. About 25% of the small haloes are missed by this prescription (Ludlow & Porciani, 2011). Secondly, for a $\Lambda$CDM power spectrum and a top-hat filter in real space (the most physically motivated smoothing), the actual calculation of the peak abundance contains divergences and cannot be carried out. Thirdly, as we discuss in this paper, a density peak does not coincide with the center of mass of the peak patch. Therefore, the motion relative to the center of mass of the particles around it may not be well approximated by spherical collapse.

From a dynamical point of view, the total energy density enclosed in some reference volume is often a more relevant quantity than the mass density. If the initial density profile is not flat, then the relation between the two overdensities is stochastic (e.g. Bond & Myers, 1996). As a result, a peak in energy overdensity is not the same as one in matter. In this paper, we study peaks in the (absolute value of the) enclosed energy density. That is, we look for the most energetically bound regions in the initial density field.

As we describe below, this approach shows considerable promise for at least three reasons. First, the mean energy overdensity governs the evolution of the moment of inertia of the collapsing patch (e.g. Chandrasekhar, 1969), which is commonly used to describe virialization. Second, for a sphere, the gradient of the mean energy is proportional to the dipole moment: when it vanishes, the center of mass coincides with the center of the sphere, and the acceleration of particles at the surface points towards it. This is therefore the spherical surface whose actual evolution (described e.g. by perturbation theory) most closely matches spherical collapse. I.e., the position around which the spherical collapse approximation is most accurate. Third, the statistics of the mean energy overdensity have better convergence properties than those of the matter overdensity: one can thus build a self-consistent peak theory without having to tweak the top hat filter as is usually done (e.g. Manrique & Salvador-Solé, 1995; Paranjape et al., 2013; Chan et al., 2017). Thus, characterizing protohaloes as peaks of the energy overdensity field addresses the problems of the matter density based excursion set peaks approach.

Section 2 discusses the physical motivation for working with energy rather than density as the primary variable. After first setting up notation and highlighting a number of interesting relations between derivatives of energy and overdensity, Section 3 derives our expression for the comoving number density of excursion set peaks in energy and describes a number of generic features of the approach. Section 4 compares a number of these predictions with the properties of protohaloes identified in simulations. A final section summarizes some consequences of our analysis and discusses a number of interesting directions for further study. Appendix A contains some technicalities of the multipole expansion, Appendix B provides expressions for some of the important correlations between matter and energy overdensities, and Appendix C discusses a modification of our energy peaks approach: a hybrid model which uses both energy and density to identify protohaloes.

2 Matter, energy and motion

2.1 Matter vs energy overdensity

Calling $\rho(\bm{r})$ the matter density at $\bm{r}$ and $\bar{\rho}$ its background value, the mean matter overdensity within a sphere of physical radius $R$ centered at the origin is usually defined as

[TABLE]

where $\delta(\mathbf{r})\equiv\rho(\mathbf{r})/\bar{\rho}-1$ and $V_{R}\equiv 4\pi R^{3}/3$. The potential energy in the sphere due to matter is (e.g. Chandrasekhar, 1969; Binney & Tremaine, 1987)

[TABLE]

where $\bm{r}_{\mathrm{cm}}$ is the center of mass position, defined as

[TABLE]

where $M$ is the total mass in $V_{R}$, and $\bm{g}\equiv-\nabla\Phi+[\nabla\Phi]_{\mathrm{cm}}$ is the acceleration at $\bm{r}$, induced by matter only, relative to the acceleration of the center of mass.

The acceleration splits into background and peculiar: $\nabla\Phi=4\pi G\bar{\rho}(\bm{r}/3+\nabla\phi)$, with the potential perturbation normalized so that $\nabla^{2}\phi=\delta$ in physical coordinates. The relative acceleration becomes

[TABLE]

the last term in the bracket being the peculiar acceleration of the center of mass $-[\nabla\phi]_{\mathrm{cm}}$. For comparison, the acceleration from the cosmological constant is simply $(\Lambda/3)(\bm{r}-\bm{r}_{\mathrm{cm}})$. The potential energy $U$ can then be written as

[TABLE]

where

[TABLE]

is the (square of the) inertial radius, and

[TABLE]

is the potential energy overdensity associated with $V_{R}$. It plays for $U$ the same role that the matter overdensity $\delta_{R}$ plays for the mass $M=\bar{\rho}V_{R}(1+\delta_{R})$. The contribution from the peculiar acceleration of the center of mass is constant over the sphere and so drops out of equation (7).

At leading order in perturbations, in equation (7) one can replace $\rho(\bm{r})/M\simeq 1/V_{R}$ and, for a spherical volume, $R_{I}\simeq R$ and $\bm{r}_{\mathrm{cm}}\simeq 0$. We can thus redefine the linearized mean energy overdensity as

[TABLE]

If $V_{R}$ is a sphere, only the radially directed monopole term of the peculiar gravitational acceleration, $\nabla\phi\sim(\delta_{r}/3)\,\bm{r}$, survives after integrating over the angles, and thus

[TABLE]

(e.g. Bond & Myers, 1996, equation 2.27b). For a homogeneous sphere, $\delta({\bm{r}})$ is constant, $R_{I}=R$ and $\delta_{R}=\epsilon_{R}$. This returns the familiar result $U=-(5/3)GM^{2}/R$. In general, however, $\delta_{R}\neq\epsilon_{R}$. Note also that, although we defined $\delta_{R}$ and $\epsilon_{R}$ using physical coordinates, they have exactly the same expression in comoving coordinates.

At a generic position ${\bm{x}}$, and for an arbitrary configuration of the density field, the matter overdensity in a spherical volume $V_{R}$ is a weighted sum over the Fourier modes:

[TABLE]

where $W_{1}(z)\equiv 3j_{1}(z)/z$. The expression for the mean linear energy overdensity at the same position follows inserting the Fourier modes $(-i\bm{k}/k^{2})\delta(\bm{k}){\rm e}^{i{\bm{k}\cdot\bm{r}}}$ of the peculiar acceleration $\nabla\phi$ in equation (8); integrating over $\bm{r}$ gives

[TABLE]

where $W_{2}(z)\equiv 15j_{2}(z)/z^{2}$. Replacing $\delta_{R}$ with $\epsilon_{R}$ simply corresponds to a change of filter. At $z\gg 1$ this filter decays faster than $W_{1}$ by one power of $z$.

As we noted in the introduction, an approach based on peaks in $\delta_{R}$ leads to divergences. Therefore, it is interesting to explore what happens if one studies peaks in $\epsilon_{R}$ instead. The more convergent asymptotic behaviour of $W_{2}$ compared to $W_{1}$ suggests that they will be mathematically better behaved. In the next subsections we argue that peaks in $\epsilon_{R}$ are also physically more reasonable than peaks in $\delta_{R}$.

2.2 Haloes as energy overdensity peaks

Any initial region is displaced and deformed by gravity. Its evolution separates as center of mass displacement plus collapse onto the center of mass. Only the latter creates a high density region, and it is this process that the spherical collapse model approximates. In this subsection, we discuss how to choose locations where this approximation is optimal.

Let us define the dimensionless dipole moment of the sphere $\bm{D}_{R}\equiv\bm{r}_{\mathrm{cm}}/R$. From equation (3), its Fourier expression at leading order (at $\bm{x}=0$ for convenience) is

[TABLE]

Having $\bm{D}_{R}\neq 0$ means that the sphere’s center of mass, towards which the acceleration converges, does not coincide with the geometric center. Hence, the collapse of the object is not isotropic, but one side collapses faster than the other. To avoid this, and find the locus of the particles that would reach the center at the same time (neglecting higher multipoles), one must shift the sphere until $\bm{D}_{R}=0$.

One can see this directly by looking at the acceleration $\bm{g}$ of a particle at the surface of the sphere relative to the center of mass, given by equation (4) with $\rho(\bm{r})/M$ replaced by $1/V_{R}$ at leading order. As we show in Appendix A, this relative acceleration at $\hat{\bm{r}}R$ can be written as

[TABLE]

using a multipole expansion. This expression contains:

(i)  a monopole term that we treat with spherical collapse;

(ii)  one anisotropic dipole term, which we ask to vanish;

(iii)  a double infinite series of higher order multipoles, which we suppose neglible in this work.

This expansion has only one dipolar anisotropy (since the center of mass acceleration was subtracted off), which can be set to zero by a translation. Dealing with higher order anisotropies would require a different transformation, like a deformation of the spherical boundary.

Combining the $k^{2}$ and $j_{2}$ factors in equation (12) gives

[TABLE]

The final equality indicates that a sphere with null dipole moment must also have $\bm{\nabla}\epsilon_{R}=0$: a convergence point of the local gravitational flow must be a stationary point in the energy overdensity field $\epsilon_{R}$. Moreover,

[TABLE]

and in particular $\nabla^{2}\epsilon_{R}=(5/R)\mathrm{d}\delta_{R}/\mathrm{d}R$, which describes the change in infall velocity as $R$ grows. For it to be slower from any direction, so that larger (outer) shells collapse later than inner ones, the Hessian $\nabla_{i}\nabla_{j}\epsilon_{R}$ must be negative definite. Thus, we seek stationary points of $\epsilon_{R}$ that are peaks.

2.3 The collapse of the inertial radius

The total energy contained inside the sphere, normalized like in equation (6) i.e. in physical coordinates, is

[TABLE]

where $K\equiv\int_{V_{R}}\mathrm{d}\bm{r}\rho(\bm{r})|\dot{\bm{r}}-\dot{\bm{r}}_{\mathrm{cm}}|^{2}/2$ is the kinetic energy and $T$ is the moment of inertia defined in equation (6). It can be conveniently written as

[TABLE]

where $k\equiv(5/6M)\!\int_{V_{R}}\!\!\mathrm{d}\bm{r}\rho(\bm{r})v_{I}^{2}$, with $\bm{v}_{I}\equiv[(\bm{r}-\bm{r}_{\mathrm{cm}})/R_{I}]\,\dot{}$, is the peculiar kinetic energy with respect to $R_{I}$ as a local scale factor, and $M_{I}\equiv 5UR_{I}/3GM=(4\pi/3)\bar{\rho}(1+\epsilon_{R})R_{I}^{3}$ is an effective mass. The true mass is $M=(4\pi/3)\bar{\rho}(1+\delta_{R})R^{3}$. Since $\delta_{R}\ll 1$ initially, $M\simeq(4\pi/3)\bar{\rho}R^{3}$.

Zel’dovich initial conditions, relating peculiar velocities to the potential perturbation, give $\dot{\bm{r}}\simeq H(\bm{r}+\nabla\phi)$ in an EdS universe. When inserting this in $K$, assuming zero curvature, the background contributions to $E$ cancel, leaving

[TABLE]

In a $\Lambda$CDM cosmology, a factor of $(3/5+2f/5\Omega_{m})$ multiplies this expression, where $f=\mathrm{d}\ln D/\mathrm{d}\ln a$ and $D$ is the growth function of linear matter perturbations ($f=1$ in EdS). For an isolated object, or for a portion of a spherically symmetric density profile, this energy is conserved.

In terms of $R_{I}$, the virial equation for the evolution of the moment of inertia, $\ddot{T}/2=2K+U+\Lambda T/3$, is

[TABLE]

Since $k$ is of second order in perturbations, and at first order $M_{I}$ and $E$ are conserved111Assuming $\dot{M}_{I}\simeq 0$ is self-consistent, since a spherical collapse solution has $R_{I}\simeq a\,(1-\epsilon_{R}/3)$ at linear order. It also follows from $\dot{U}=-\int_{V_{R}}\!\mathrm{d}\bm{r}\,\rho(\bm{r})\,\dot{\bm{r}}\cdot\bm{g}$ (note the sign!) for an isolated body (Chandrasekhar, 1969; Binney & Tremaine, 1987), so that $\dot{M}_{I}\propto-\int_{V_{R}}\!\mathrm{d}\bm{r}\,\rho(\bm{r})\,\bm{v}_{I}\cdot\bm{g}\simeq 0$ ar first order. The coupling of a sphere with its environment appears only at second order., equations (17), (18) and (19) coincide (up to first order) with the equations of spherical collapse for $R$, except that one replaces $\delta_{R}\to\epsilon_{R}$. Thus, $R_{I}$ approximately follows a spherical collapse solution for overdensity $\epsilon_{R}$, with corrections starting at second order in perturbation theory. Hence, $\epsilon_{R}$ plays for the inertial radius $R_{I}$ the role that $\delta_{R}$ plays for $R$: it indicates if and how fast the moment of inertia of an extended region is destined to shrink. Therefore, $\epsilon_{R}$ may be at least as good an indicator as $\delta_{R}$ for inferring if an initial patch is destined to become a nonlinear overdense halo. Although we remarked earlier that $\delta_{R}\neq\epsilon_{R}$ in general, we expect they will be correlated with one another. We quantify this and a number of other relevant correlations, in the next section.

3 Excursion set peaks of the mean energy overdensity

Having laid out the physical motivation for studying peaks in $\epsilon_{R}$, we now discuss their statistics.

3.1 Normalized variables

In this subsection we lay out the notation that will be needed in what follows. We begin by defining the generic filter

[TABLE]

such that $W_{n}(z)\to 1$ as $z\to 0$. The recursion relations of the spherical Bessel functions imply that

[TABLE]

which we use below. We also define the variances

[TABLE]

Because $W_{2}$ decays faster than $W_{1}$ at large $k$, $\sigma_{j2}^{2}$ may converge even when $\sigma_{j1}^{2}$ does not. We will see that, as a result, peaks theory for $\epsilon_{R}$ is better behaved than for $\delta_{R}$. The exact behavior of the different filters is shown in Fig. 1. Setting $z=kR$ and Fourier transforming $W_{n}(kR)$ yields

[TABLE]

Evidently, in real space, filters with larger $n$ weight the central regions more than the outer parts.

Secondly, we introduce the normalized fields

[TABLE]

Next, we need the (normalized) gradient

[TABLE]

and the (negative of the) normalized Hessian

[TABLE]

where $\sigma_{12}^{2}(R)=\langle{\bm{\nabla}}\epsilon_{R}\cdot{\bm{\nabla}}\epsilon_{R}\rangle$ and $\sigma_{22}^{2}(R)=\langle[\nabla^{2}\epsilon_{R}]^{2}\rangle$ were defined in equation (23). These definitions differ from what is commonly used in the literature (e.g. Bardeen et al., 1986) only because we use the filter $W_{2}$.

We then define the curvature $x_{R}\equiv{\rm Tr}(\zeta_{R})$ as the trace of the (negative of the) normalized Hessian:

[TABLE]

normalized so that $\langle x_{R}^{2}\rangle=1$. Finally, we introduce the (normalized) ‘slope’ of the mean energy overdensity

[TABLE]

The last equality follows from case $n=2$ of equation (21).

For what is to follow, it is useful to note that the curvature $x_{R}$ can also be written as

[TABLE]

(equation 21 with $n=1$). The derivative $-{\rm d}\delta_{R}/{\rm d}\ln R$ on the right hand side is sometimes called the ‘slope’ of the mean matter overdensity. In the standard analysis, the ‘slope’ and ‘curvature’ $-\nabla^{2}\delta_{R}/\sigma_{21}(R)$ of the overdensity are the same for Gaussian smoothing filters. The expression above shows that the slope of the tophat smoothed overdensity is the same as the curvature of the mean energy overdensity.

This is significant because the slope of the mean matter overdensity plays an important role in the ‘upcrossing approximation’ of the traditional excursion set approach (Musso & Sheth, 2012, although their ‘slope’ variable is really ${\rm d}\delta_{R}/{\rm d}\sigma_{01}^{2}$). The simplest versions of this approach require ${\rm d}\delta_{R}/{\rm d}\ln R\leq 0$ as this implements the request that the outer shells collapse after the inner ones. Equation (30) shows that by constraining this slope, upcrossing constrains the curvature of the mean energy overdensity $x_{R}$.

We had argued that we are interested in stationary points of $\epsilon_{R}$: positions where ${\bm{\nabla}}\epsilon_{R}=0$ (c.f. equation 14). The ‘upcrossing’ constraint $\mathrm{d}\delta_{R}/\mathrm{d}\ln R<0$ of equation (30) adds the additional requirement that $\mathrm{Tr(\zeta_{R})}>0$. Since the same argument (about outer shells collapsing later than the inner ones) actually applies to the full $\zeta_{ij}$, not just to its trace, $\nabla_{i}\nabla_{j}\epsilon_{R}\propto-\zeta_{ij}$ must be negative definite. Thus, we are interested in those stationary points of $\epsilon_{R}$ that are peaks.

Furthermore, for the ‘slope’ of the energy one has that

[TABLE]

Evidently, requiring ${\rm d}\epsilon_{R}/{\rm d}\ln R\leq 0$ imposes that $\delta_{R}\leq\epsilon_{R}$. Physically, this constraint guarantees that the inertial radius $R_{I}$, which is initially equal to $R$, gradually becomes smaller than $R$ as the density profile develops a denser core. We exploit this in the next subsection.

In what follows, we will be interested in whether or not $\epsilon_{R}$ exceeds a critical value. With some abuse of notation, we use $\epsilon_{c}(R)$ to denote the fact that this critical value may be different for each $R$. The associated normalized height and slope of this critical value are given by

[TABLE]

3.2 Peaks in the $\epsilon_{R}(x)$ field

Since $\epsilon_{R}(\bm{x})$ is a Gaussian field, the comoving number density of peaks of height $\epsilon_{c}$ at fixed comoving scale $R$ can be computed following the steps laid out by Bardeen et al. (1986) for $\delta_{R}$ (see also Lazeyras et al., 2016, for a shorter derivation based on rotational invariants). Neglecting for simplicity the subscript $R$ on variables, the result is formally written as

[TABLE]

where $V_{*}\equiv[6\pi(\sigma_{12}/\sigma_{22})^{2}]^{3/2}$ defines a characteristic volume, $p(x|\omega_{c})$ is the conditional distribution of $x$ given $\omega=\omega_{c}$, and

[TABLE]

as given by equation (A15) of Bardeen et al. (1986), follows from integrating over the traceless part of $\zeta_{ij}$. Therefore, $[f(x)/V_{*}]\,p(x,\omega_{c})$ is the number density of peaks of height $\omega_{c}$ and curvature $x$. We note in passing that, had we not integrated over the traceless part of $\zeta_{ij}$, we would have expressions for peak ‘ellipticity’ and ‘prolateness’, exactly as in Appendix C of Bardeen et al. (1986), but with our $W_{2}$ smoothing window.

Bayes’ rule can be used to show that $p(x|\omega_{c})$ is the Gaussian distribution

[TABLE]

where the cross-correlation coefficient $\gamma_{x\omega}$ is

[TABLE]

Equation (33) describes the differential distribution of the number density of energy peaks of height $\epsilon_{c}$ when the smoothing scale is $R$. However, it does not say how $R$ must be chosen. Since a given $R$ contains mass $\propto 4\pi R^{3}/3$, ${\rm d}n/{\rm d}\epsilon_{c}$ for fixed $R$ is not obviously related to a distribution of halo masses. To remedy this shortcoming of peaks theory, we must combine it with the excursion set approach.

3.3 Excursion sets in $\epsilon_{R}$

To set the smoothing scale, we add the ‘excursion set’ requirement that at a fixed position $\bm{x}$ the smoothed field $\epsilon_{R}$ equals the critical value $\epsilon_{c}$ on scale $R$, but is smaller on the next larger smoothing scale (the ‘upcrossing’ approximation of Musso & Sheth, 2012). As noted in Section 2.1, this constraint on the slope $y_{R}$ ensures that the inertial radius $R_{I}$ shrinks relative to $R$ as the object collapses and its density profile steepens and peaks around its center of mass. As $\epsilon_{c}$ may depend on $R$ in general, the upcrossing condition is $\epsilon_{R}\geq\epsilon_{c}(R)$ but $\epsilon_{R+\Delta R}\leq\epsilon_{c}(R+\Delta R)$. To first order in $\Delta R$, and using the normalized slope from equation (29), this means $\omega_{c}\leq\omega_{R}\leq\omega_{c}+(R\sigma_{23}/7\sigma_{02})(y-y_{c})\Delta R$ and $y\geq y_{c}$, where $y_{c}$ is given by equation (32).

Integrating over the allowed ranges and dividing by $\Delta R$ gives the distribution of $R$ values at which $\epsilon_{R}$ upcrosses:

[TABLE]

where

[TABLE]

is the Gaussian conditional distribution for $y$, with

[TABLE]

This result closely resembles that of equation (33), even more so when $y_{c}=0$ (for constant $\epsilon_{c}$).

Using $Rf_{\mathrm{up}}(R)=\sigma_{02}\,f_{\mathrm{up}}(\sigma_{02})\,|{\rm d}\ln\sigma_{02}/{\rm d}\ln R|$ and integrating over $y$ gives

[TABLE]

where

[TABLE]

reduces to $\omega_{c}$ for a constant barrier, and

[TABLE]

Equation (40) describes the distribution of smoothing radii $R$ satisfying the upcrossing constraint when the smoothing filters are centered on randomly chosen positions in the field. In particular, these positions need not be peaks.

3.4 Excursion set peaks of $\epsilon_{R}$

The excursion set peaks approach (Paranjape & Sheth, 2012; Paranjape et al., 2013) posits that the comoving number density of protohalo patches of initial size $R$ (that will become haloes of mass $M=\bar{\rho}\,4\pi R^{3}/3$) can be estimated by combining the excursion sets and peaks analyses of the previous two subsections. I.e., we are interested in the distribution in $R$ of points that on scale $R$ are both upcrossing $\epsilon_{c}$ and are peaks. This comoving density is given by

[TABLE]

with $\tilde{y}\equiv y-y_{c}$. This expression combines equations (33) and (37).

The conditional probability $p(\tilde{y},x|\omega_{c})$ can be factorized using Bayes’ rule as $p(\tilde{y}|x,\omega_{c})\,p(x|\omega_{c})$, where $p(x|\omega_{c})$ equals the probability of $x-\gamma_{x\omega}\omega$ with $\omega=\omega_{c}$ (c.f. equation 35), and $p(\tilde{y},x|\omega_{c})$ is the conditional distribution of $y-\gamma_{y\omega}\omega$ given $x-\gamma_{x\omega}\omega$, up to a shift by $y_{c}$. That is,

[TABLE]

where

[TABLE]

and

[TABLE]

Integrating over $\tilde{y}$, and using the fact that ${\rm d}n/{\rm d}\ln M=(R/3)\,{\rm d}n/{\rm d}R$ because $M\propto R^{3}$, the excursion set peaks expression for the mass function is

[TABLE]

where

[TABLE]

and the function $F$ is defined by equation (42).

In the high peaks limit, $p(x|\omega_{c})$ becomes sharply peaked around its mean $\gamma_{x\omega}\omega_{c}$ and, since $f(x)\sim x^{3}$ for $x\gg 1$, the integral over $x$ tends to $(\gamma_{x\omega}\omega_{c})^{3}\gamma_{y\omega}\,\omega_{*}$. Moreover,

[TABLE]

as one can see by combining equations (41) and (90) for the quantities on the left hand side, and using $n_{*}V_{*}$ from equation (50). As a result,

[TABLE]

The last expression shows that although the full calculation (equation 49) depends on $\sigma_{22}$, the high peaks limit does not.

Since ${\rm d}\ln M/{\rm d}\ln R^{3}=1$, this expression has the same form as the high peaks limit of equation (33) for ${\rm d}n/{\rm d}\omega_{c}$. However, there $\omega_{c}$ denotes peak height on a fixed scale whereas here, $\omega_{c}=\epsilon_{c}/\sigma_{02}(R)$ is a proxy for the variable scale $R$ (and hence, for $M$).

Formally, peaks and excursion set peaks in the enclosed matter overdensity are still described by equation (33) and (43) respectively, replacing $\epsilon_{R}$ with $\delta_{R}$. That is, $\omega_{c}\to\nu_{c}$ and all $\sigma_{j2}\to\sigma_{j1}$. We say ‘formally’ because in practice, for many $P(k)$ of current interest, the integral defining $\sigma_{21}$ does not converge (whereas that for $\sigma_{22}$ does), so the calculation is inconsistent. (E.g., the resulting expressions are ill defined since $V_{\star}=0$.) For practical calculations, one has to regularize this divergence somehow, for instance by replacing the top-hat filter $W_{1}$ with a more convergent one (e.g. Gaussian). This tweak is not needed for our energy based approach. (In this respect, the high peaks limit is less problematic than the full calculation, since the most divergent quantity is not required. Physically, this happens because high stationary points are likely all peaks, so the constraint on the sign of the Hessian is irrelevant.)

3.5 A more formal derivation

We give here a more elegant derivation of the results of the previous subsections, using only the transformation rule of Dirac’s delta functions and treating excursion sets and peaks on the same footing. The number density of the stationary points of a single realisation of the field $\epsilon_{R}(\bm{x})$ is

[TABLE]

where $\bm{x}_{\mathrm{sp}}$ denotes each solution of the constraint equation $\eta_{i}(\bm{x})=0$, and $\det(\nabla\bm{\eta})=(-\sigma_{22}/\sigma_{12})^{3}\det(\zeta)$ is the Jacobian of the transformation from $\bm{\eta}$ to $\bm{x}$. Similarly, the number density in $[0,+\infty)$ of the solutions $R_{c}$ (‘crossings’) of the constraint equation $\epsilon_{R}(\bm{x})=\epsilon_{c}$ is

[TABLE]

where $|\mathrm{d}(\epsilon_{R}-\epsilon_{c})/\mathrm{d}R|$ transforms the number density in $\epsilon_{R}$ into one in $R$.

To make equation (53) a density of peaks, one must count only the solutions where the eigenvalues $\zeta_{i}$ of $\zeta_{ij}$ are all positive. To make it a differential density of peaks of fixed height $\omega_{c}$, one must further multiply by $\delta_{\rm D}(\omega_{R}-\omega_{c})$. Hence, the differential mean density of peaks of height $\epsilon_{c}$ is

[TABLE]

where the angle brackets denote an average over the joint distribution of $\omega,{\bm{\eta}},{\bm{\zeta}}$ on scale $R$. This average gives equation (33) for a fixed smoothing scale $R$. Integrating this over $\epsilon_{c}$ gets rid of $\delta_{\rm D}(\epsilon_{R}-\epsilon_{c})$, since integration and mean commute.

Similarly, to make equation (54) a density of upcrossing points, the derivative $\mathrm{d}(\omega_{R}-\omega_{c})/\mathrm{d}R$ must be negative. The mean distribution in $R$ of such points is thus

[TABLE]

where ${}^{\prime}=\mathrm{d}/\mathrm{d}R$, and the angle brackets denote an average over the joint distribution of $\omega_{R}$ and $\omega_{R}^{\prime}$. Note that $\langle\omega\omega^{\prime}\rangle=0$ (since $\omega$ is normalized) and direct calculation shows that $\omega^{\prime}\propto y-\gamma_{y\omega}\omega$ is the relevant quantity for $p(y|\omega_{c})$. Thus, the result is equation (37) up to a total derivative $\mathrm{d}\sigma_{02}/\mathrm{d}R$ to change variables from $R$ to $\sigma_{02}(R)$.

When combining peaks and excursion sets, the problem becomes that of finding the density of points in the 4-dimensional space $\{R,\bm{x}\}$ that solve $\{\omega,\bm{\eta}\}=\{\omega_{c},\bm{0}\}$ (plus the peak and upcrossing constraints). That is,

[TABLE]

where $J=\det[\partial\{\omega_{R}-\omega_{c},\bm{\eta}\}/\partial\{R,\mathbf{x}\}]$ is the 4-dimensional Jacobian determinant. Since $\nabla(\omega_{R}-\omega_{c})=0$ at peaks, $J$ factorizes into $|J|=(\sigma_{22}/\sigma_{12})^{3}|\det(\zeta)||\mathrm{d}(\omega_{R}-\omega_{c})/\mathrm{d}R|$. Equation (57) is the same as equation (43).

3.6 Peak velocities and velocity bias

We have emphasized the fact that our energy-peaks excursion set approach is a simple but well-motivated modification of the density-based approach. As Bardeen et al. (1986) noted, the peaks constraint modifies the distribution of velocities. This is because, although velocities do not correlate with (even space derivatives of) the density, they do correlate with odd derivatives.

For Gaussian statistics, the constrained mean is linear in the constrained variables. Since the gradient vanishes at a peak, the mean peak velocity is zero (in each direction). However, the dispersion around this mean is modified by the constraint on the gradient. Since the Fourier transform of the velocity is given by ${\bm{v}}=i{\bm{k}}\,\delta({\bm{k}})/k^{2}\,W_{1}(kR)$ and the gradient of $\epsilon_{R}$ is ${\bm{\nabla}}\epsilon=i{\bm{k}}\,\delta({\bm{k}})\,W_{2}(kR)$, their cross-correlation is

[TABLE]

and the peak velocity dispersion is

[TABLE]

Thus,

[TABLE]

Since $\sigma^{2}_{\rm vpk}\neq\sigma^{2}_{-11}$ this manifests as ‘velocity bias’ in Fourier space (Desjacques & Sheth, 2010). ‘Velocity bias’ is unfortunate nomenclature because it suggests that the velocity of a peak identified on scale $R$ differs from the mean motion of the dark matter within $R$, when in fact this is the definition of the speed of the peak patch. The ’bias’ is really with respect to the statistics of dark matter speeds and arises from two separate effects: the dependence on smoothing scale $R$, and the peak constraint, which gives rise to the terms which correct $\sigma_{-11}$.

3.7 Scale dependence and scatter in $\delta_{R}$

Because of equation (31), $y$ and $\delta$ are deterministically (and linearly) related at fixed $\epsilon$. Thus, a scatter in $y$ automatically produces a scatter in $\delta$. In particular, the integral over $\tilde{y}$ in equation (43) shows explicitly that peaks which are upcrossing on the same $R$ can have a range of $\delta$, weighted by $p(\tilde{y}|x,\epsilon_{c})$. That is to say, a model in which haloes of the same mass all have the same $\epsilon_{c}$ generically predicts scatter in the enclosed overdensities of the protohalo patches, with no additional work (i.e., equation 43 is no more complicated than the excursion set peaks model for $\delta$, in which $\delta_{c}$ is deterministic rather than stochastic).

The moments of the conditional distribution of $\delta$ given the ESP constraint (peaks plus upcrossing) on scale $R$ can be obtained inserting powers of

[TABLE]

in equation (43), and then dividing by $\mathrm{d}n/\mathrm{d}R$. For instance, the mean value and variance of $\delta$ at fixed $R$ are

[TABLE]

where

[TABLE]

and $p(\tilde{y}|x,\epsilon_{c})$ is given by equation (44).

For any $n$, the integral over $\tilde{y}$ in equation (64) equals

[TABLE]

where ${}_{1}F_{1}(a,b,z)$ is the confluent hypergeometric function, and $\mu$ and $\Sigma$ are the mean and variance of $p(\tilde{y}|x,\epsilon_{c})$, given by equations (46) and (47) respectively. For $n=1$ this returns $\mu\,F(\mu/\Sigma)$, with $F$ defined as in equation (42); for $n=2$ it gives

[TABLE]

and for $n=3$

[TABLE]

The high peak limit of these expressions is recovered by setting $\mu\simeq\gamma_{y\omega}\omega_{*}\gg\Sigma$, in which case one gets that $\langle\delta|\mathrm{ESP}\rangle\to\epsilon_{c}(1+\mathrm{d}\ln\sigma_{02}/\mathrm{d}\ln R^{5})<\epsilon_{c}$ and $\mathrm{Var}(\delta|\mathrm{ESP})\to(R^{2}\sigma_{23}\Sigma/35)^{2}$.

We show in the Appendix that $R^{2}\sigma_{23}\propto\sigma_{02}$, with a correction factor that is independent of $R$ for power law power spectra. In addition, $\sigma_{02}\propto\sigma_{01}$ for such spectra, so we generally expect rms$(\delta|\mathrm{ESP})$ to be approximately $\propto\sigma_{01}$, which is in agreement with the scatter around the mean overdensity of protohaloes identified in simulations (Dalal et al., 2008; Robertson et al., 2009; Elia et al., 2012; Despali et al., 2013). Figure 2 shows a plot of $\langle\delta|\mathrm{ESP}\rangle$ with the rms spread around this mean, for $\epsilon_{c}=1.686$. Notice that the mean is smaller than $\epsilon_{c}$, a consequence of the upcrossing constraint that $y_{R}\propto(\epsilon_{R}-\delta_{R})\geq 0$. Moreover, the mean is a decreasing function of $\sigma_{02}$, and hence of $\sigma_{01}$. In contrast, the mean overdensity within protohalo patches identified in simulations increases approximately as $1.686+0.4\sigma_{01}$ (Sheth et al., 2001, 2013). Therefore, we might naively expect $\epsilon_{c}$ to increase even more strongly with $\sigma_{02}$. We test this prediction below.

4 Comparison with simulations

This section describes a few simple tests of our formalism using 5378 haloes identified in the $z=0$ outputs of the Flora simulation, the largest box in the SBARBINE suite (Despali et al., 2016). Our halo set is composed as follows: all the 1378 halos more massive than $10^{15}h^{-1}M_{\odot}$, 2000 randomly chosen haloes with mass between $10^{14}$ and $10^{15}h^{-1}M_{\odot}$, and 2000 randomly chosen haloes with mass between $4\times 10^{13}$ and $10^{14}h^{-1}M_{\odot}$.

The simulation followed the gravitational evolution of $1024^{3}$ dark matter particles each of mass $6.35\times 10^{11}h^{-1}M_{\odot}$ in a periodic cube of side $L_{\mathrm{box}}=2h^{-1}$Gpc with a Planck13 background cosmology: $\Omega_{m}=0.307$, $\Omega_{\Lambda}=0.693$, $\sigma_{8}=0.829$ and $h=0.677$.

Haloes were identified using a Spherical Overdensity halo finder with threshold of $319\times$ the background density. We use ‘protohalo’ to refer to the patch defined by a halo’s particles in the initial conditions. This patch is not spherical, whereas our formalism is based on spherically symmetric smoothing filters. To measure spherically averaged quantities in the initial conditions, for each protohalo having e.g. $N_{\rm halo}$ particles we identify a ‘protosphere’: this is the set of all particles closer than $R=(3N_{\rm halo}/4\pi)^{1/3}(L_{\rm box}/1024)$ to the protohalo’s center of mass.

4.1 Estimators

In practice, we compute averages by explicitly summing over real space positions, rather than using Fourier methods. This means that we estimate

[TABLE]

where the sum is over the particles in the protohalo or in the protosphere (compare equation 7, and recall that $\nabla\phi=-v/fH$ initially). Note that the denominator is close to $(3/5)\,N_{\mathrm{halo}}R^{2}$. Similarly, we estimate

[TABLE]

where now $j$ runs only over the particles in a thin spherical shell of radius $R$ (e.g. $0.95R<|\bm{r}-\bm{r}_{\mathrm{cm}}|<1.05R$). (We only retain haloes with more than 10 particles in this shell.) We have checked that this measurement agrees with $\hat{\epsilon}_{R}+({\rm d}\hat{\epsilon}_{R}/{\rm d}\ln R)/5$ (see equations 29 and 31) where the derivative is estimated by first evaluating equation (68) for many narrowly spaced $R$ values (while keeping $\bm{r}_{\mathrm{cm}}$ and $\bm{v}_{\mathrm{cm}}$ fixed to their values for the full sphere).

4.2 Peaks in the initial $\epsilon_{R}$ and $\delta_{R}$ distributions

Our first test is to check if $\epsilon_{R}$ is (close to) a local maximum. For this, we generate a fine cartesian grid with spacing $0.02R$, out to $\pm 0.2R$ from the center of mass position. We measure $\epsilon_{R}$ at each grid point using equation (68), where $j$ labels all particles within $R$ from the grid point. Figure 3 shows the $\epsilon_{R}$ field in the central regions (i.e. within $0.2R$ in the $x-y$ plane at $z=0$) of 25 randomly selected protohaloes, ordered by mass. The integers in the top left panels show the protohalo’s mass in units of $10^{13}h^{-1}M_{\odot}$. The color scale is normalized to the maximum value of $\epsilon$ in each panel. It is clear that the most massive objects (top) are noticeably better centered and smoother than the lower mass ones (bottom).

Maps of $\delta_{R}$ for these same regions are shown in Figure 4. While it is again true that the most massive objects are better centered and smoother, these maps are clearly much less smooth, and not nearly as well centered as the corresponding maps of $\epsilon_{R}$. The difference is particularly striking at the lower masses shown in the bottom two or three rows.

To emphasize this point, Figure 5 shows the ratio $\epsilon_{R}(\Delta x)/\epsilon_{R}(0)$ (left) and $\delta_{R}(\Delta x)/\delta_{R}(0)$ (right) along the $x$-axis, where $\Delta x$ is the offset from the protohalo’s center of mass, for 100 randomly selected protohaloes. Note that the $\epsilon_{R}$ profiles are more likely to curve downwards. They also tend to be smoother; the $\delta_{R}$ profiles can have several peaks within a small distance of the center of mass position, consistent with the fact that the variance of $\nabla^{2}\delta_{R}$ diverges. This illustrates why $\epsilon_{R}$ is better suited for identifying protohaloes.

Finally, Figure 6 shows an estimate of the correspondence between the local maximum in $\epsilon_{R}$ or $\delta_{R}$ and the protohalo center. For each protohalo, we first smooth $\epsilon$ with a filter of scale $R$, at a number of positions within $0.3R$ of the protohalo center. If our model is correct, then the position of the largest of these smoothed $\epsilon_{R}$ values should coincide with the protohalo center. The orange curves in Figure 6 show the cumulative distribution of the distance (in units of $R$) between the position of the largest value within $0.3R$ and the protohalo center, for the halos in our three broad mass bins. (We only show scales $r\leq 0.25R$ because all the curves must reach unity at $0.3R$, by definition.) There is a clear trend with mass – massive halos are better centered – so this quantifies the obvious mass dependence in Figure 3. Note that $r/R=0.215$ corresponds to the innermost 1% of the volume, so the fact that more than 80% of the massive objects have their maximum $\epsilon_{R}$ in the innermost 1% of the volume is remarkable. For lower masses, this fraction is approximately 50%. As a converence check, we have repeated the exercise, but now restricting our search for a local maximum to $0.2R$. In this case, all the curves must reach unity at $0.2R$, so the symbols only show the distributions out to $0.1R$. They lie on the smooth curves, indicating convergence.

The blue curves in Figure 6 show a similar analysis of $\delta_{R}$. Comparison with the orange curves shows that, for all three mass bins, the maximum in $\delta_{R}$ tends to lie substantially further from the center than it does for $\epsilon_{R}$. The blue crosses show the result of restricting our search for a local maximum to $0.2R$ (rather than $0.3R$), for the halos in the middle mass bin. It is clearly offset from the dashed blue curve, indicating that the procedure has not converged. Had we repeated the analysis using maxima within $0.4R$ instead, the curves would all have smaller values of $p(<r/R)$, at least for $r/R<0.1$ or 0.2. These lower values, and the lack of convergence, are yet another way of illustrating why $\epsilon_{R}$ is better suited for identifying protohaloes than is $\delta_{R}$.

4.3 Collapse thresholds, slope and stochasticity

Our next test is to check if the excursion set quantity $y_{R}\propto\epsilon_{R}-\delta_{R}$ (the slope of the trajectory, c.f. equation 31) is strictly positive. Figure 7 shows that it is: $y_{c}\approx 2\pm 0.7$ for all protospheres.

Following Section 3.7, yet another test is to check if $\epsilon$ is a stronger function of protohalo mass than is $\delta$. The top panel of Figure 8 shows $\delta_{R}$ in the protosphere of radius $R$ centered on the protohalo center of mass for a number of protohaloes, each labeled by its $\sigma_{01}(R)$ (recall that $M\propto R^{3}$ increases as $\sigma_{01}(R)$ decreases). The solid line shows the linear regression

[TABLE]

which provides the best fit to the mean trend. This is qualitatively consistent with, although steeper than, the 0.48 scaling reported in previous literature (e.g. Sheth et al., 2001; Robertson et al., 2009).

The rms scatter around the mean scales approximately as $0.32\,\sigma_{01}$, and is also consistent with previous work (Robertson et al., 2009). The dark and light shaded regions show the mean plus and minus one or two standard deviations (i.e., the 68% and 95% prediction bands). The bottom part of this panel shows the data in the format used by Despali et al. (2013); this removes the scaling of the mean and rms with $\sigma_{01}$ so the mean height, 0.63, equals the slope of equation (70) and the amplitude of the rms scatter around it gives the scaling of the rms with $\sigma_{01}$, 0.32. In this panel, darker colors show objects that have larger values of $\epsilon_{R}$: as expected, $\epsilon_{R}$ and $\delta_{R}$ are strongly correlated.

In the lower panel of Figure 8, the orange symbols show a similar analysis of $\epsilon_{R}$ as a function of $\sigma_{02}$. The solid line shows the linear regression

[TABLE]

This is steeper than equation (70), in agreement with the discussion in Section 3.7. The rms scatter around this relation increases as mass decreases (qualitatively like for $\delta_{R}$: in this case it scales as $0.32\,\sigma_{02}$. The blue symbols in the panel show the result of estimating $\epsilon_{R}$ using the protohalo (rather than protosphere) particles. They tend to lie above the orange ones, consistent with the fact that each protosphere likely contains some particles that did not make it into the halo and are therefore likely to be less bound than the full set of protohalo particles. The bottom part of the panel shows that the slope of the $\epsilon_{R}-\sigma_{02}$ relation is 0.8, and the differently colored symbols show that the scatter correlates with $\delta_{R}$. As a result, at fixed $\delta_{R}$ the scatter in $\epsilon_{R}$ is smaller.

4.4 Velocities

Finally, we have compared the velocities of protosphere and protohalo patches with our predictions. The rms velocity is smaller for more massive protospheres, and is always smaller than $\sigma_{-11}$. This is in qualitative agreement with previous measurements which showed that the prediction based on $\delta_{R}$-peaks (set all $W_{2}\to W_{1}$ in our equation 59) is quite accurate (Sheth & Diaferio, 2001). The $\delta_{R}$- and $\epsilon_{R}$-peak predictions only differ by a few percent so, to increase signal to noise, Figure 9 shows the distribution of $v_{\rm 1d}/\sigma_{\rm 1d}$ in protospheres, where $v_{\rm 1d}$ represents each of the three cartesian components of the velocity, and $\sigma_{\rm 1d}^{2}=\sigma_{\rm vpk}^{2}/3$. If the mass dependence of $\sigma_{\rm 1d}$ is correct, the resulting distribution should be Gaussian with zero mean and unit rms.

The smooth black curve shows such a unit variance Gaussian. Orange triangles and blue squares show our measurements when $\sigma_{\rm 1d}$ is that for $\epsilon_{R}$- and $\delta_{R}$-peaks; they are very similar to the Gaussian. (Statistical errors on the measurements are similar to the symbol sizes shown.) The black dotted curve uses $\sigma_{\rm 1d}=\sigma_{-11}$ which includes the effect of smoothing the velocities but not the additional terms coming from the peaks constraint: it is clearly inconsistent with the Gaussian shape. As results for protohaloes are almost identical, we conclude that the peaks-based approach provides a good description of protohalo speeds, but the differences between $\epsilon_{R}$- and $\delta_{R}$-peaks are too small to matter.

4.5 Halo abundances: An illustrative example

This subsection presents a comparison of the halo mass function measured in the $z=0$ output of the simulation with that predicted by our approach. The simplest prediction, equation (49), requires knowledge of how $\epsilon_{c}$ depends on smoothing scale $R$ (and hence on mass). We could use equation (71) for this, but, as Figure 8 shows, there is substantial scatter (rms$=0.3\,\sigma_{02}$) around the mean trends with mass. Moreover, this scatter correlates with $\delta_{R}$, and hence with $y_{R}\propto\epsilon_{R}-\delta_{R}$. Figure 10 shows this explicitly: objects with steeper slopes tend to have larger $\epsilon_{R}$. At fixed $y_{R}$, the scatter is reduced slightly to rms$=0.25\,\sigma_{02}$.

To include this correlation we assume that

[TABLE]

where $q_{R}$ is a random variate that we assume is uncorrelated with either of $\epsilon_{R},y_{R},x_{R}$ and accounts for the scatter in $\epsilon_{c}$ at fixed $y_{R}$ shown in Figure 10. This is motivated by previous work with $\delta_{R}$, in which $q_{R}$ is related to the amplitude of the shear field (Sheth & Tormen, 2002; Sheth et al., 2013). The shear statistics are $\chi^{2}$ with 5 degrees of freedom. To highlight the fact that we have not yet built shear statistics into our model, rather than using $\chi^{2}_{5}$ statistics, we assume that $q_{R}=\sigma_{02}(R)\,q$, where $q$ is drawn from a Lognormal with unit mean and rms chosen to match Figure 10; this makes our $q_{R}$ more like the parameter $\beta$ in previous ESP work with $\delta_{R}$ (e.g. Paranjape et al., 2013). Matching Figure 10 requires $(a_{0},a_{1},a_{2})=(1.71,0.21,0.42)$ and the rms of $q$ to be 0.25.

Then, equation (43) – or (49) – should be replaced by

[TABLE]

where the mass function at fixed $q$ is

[TABLE]

$\omega_{c}$ is $\epsilon_{c}$ of equation (72) divided by $\sigma_{02}$, and $Y$ is the domain where $y-y_{c}>0$. Note that now equation (32) gives

[TABLE]

with222These expressions assume that the derivative acts only on the explicit $R$-dependence in $y_{R}(\bm{x})$ and $q_{R}(\bm{x})$, and not on the scale dependence of the peak position $\bm{x}$. We show in Appendix C that this is equivalent to assuming that $\bm{x}$ is also a peak for $\delta_{R}$; otherwise, an additional term $(-7/R^{2}\sigma_{23})(5/\sigma_{22})\,(\nabla_{i}\delta\,\zeta_{ij}^{-1}\nabla_{j}\delta$) would appear (see equation 114, and recall that $\nabla_{j}y\propto\nabla_{j}(\epsilon-\delta)\propto-\nabla_{j}\delta$). We further set $\nabla_{i}\delta\,\zeta_{ij}^{-1}\nabla_{j}\delta\approx 1.5\sigma_{11}^{2}(1-\gamma^{2}_{\nabla_{\delta}\nabla_{\epsilon}})/(x/3)$, where the numerator is the conditional variance of $\nabla\delta$ given that $\nabla\epsilon=0$, we approximated $\zeta_{ij}\approx\delta_{ij}x/3$ (retaining only the trace part), and 1.5 is a fudge factor to account for a similar term in equation (77).

[TABLE]

from equations (31), (29) and (30), and

[TABLE]

Equation (75) shows that assuming $\epsilon_{c}(y,q)$ brings additional dependence on $x$ and $y$ (as well as $q$) into $y_{c}$, which affects the predicted abundances as well as mean and variance of the typical overdensity of protohalo patches (equations 62 and 63). Some integrals in this example actually have analytical expressions. However, we do not give them here as the main point is simply to illustrate that our energy-peaks approach can be extended in a relatively transparent manner to include more involved collapse criteria. Exploring what these criteria are is work in progress. To motivate this exploration, Figure 11 shows that equations (73) and (72), with $(a_{0},a_{1},a_{2})=(1.75,0.21,0.42)$, provides a good description of halo abundances in our simulation.

5 Discussion and conclusions

The assumption that haloes form from the spherical collapse of a homogeneous sphere has motivated the development of models in which haloes are identified with sufficiently dense peaks in the initial overdensity fluctuation field. For a homogeneous sphere, the enclosed overdensity (equation 1) and energy (equation 17) are the same. However, if the sphere is not homogeneous – either because of substructure within it, or simply because its radial density profile, while smooth, is not flat – then they are not.

We argued that a sphere in a realistic (non-spherical) density field follows spherical collapse as closely as possible if centered at a minimum of the enclosed total energy, rather than maximum of the mean matter density. At such minima, the center of convergence of the gravitational flow and its geometrical center coincide (equation 12). We showed how to modify the excursion set peaks approach if haloes are identified with peaks in energy (equation 11) rather than enclosed density (equation 10). In this respect, it is natural to call ours the $\epsilon$ESP approach, and the traditional approach $\delta$ESP. However, to ease notation, we simply used ESP throughout.

Our analysis and result (equation 43) are no more complicated than the usual overdensity-based analysis. The model naturally predicts scatter in the matter overdensity $\delta_{R}$ of protohalo patches (equation 63) around a mean value. This mean decreases with $R$ (i.e., at lower masses) (equation 62 and Figure 2) if protohalo patches all have the same critical value of the energy overdensity $\epsilon_{c}$ whatever their size $R$. However, both N-body simulations and theoretical arguments suggest that $\delta_{R}$ must be larger at smaller masses. Therefore, our model predicts that $\epsilon_{c}$ must increase even more strongly as protohalo patch size decreases.

Measurements in simulations confirm many of the generic predictions of our $\epsilon$ESP approach. Maps of $\delta_{R}$ and $\epsilon_{R}$ show that protohalo patches are much more likely to be centered on $\epsilon_{R}$- than $\delta_{R}$-peaks (Figures 3–6). The excursion set prediction that $\epsilon_{R}\geq\delta_{R}$ (equation 31) is also confirmed (Figure 7), as is the prediction that $\epsilon_{R}$ should be a steeper function of mass than is $\delta_{R}$ (Figure 8 and equations 70 and 71). Although we only tested our formalism on halos more massive than $10^{13}h^{-1}M_{\odot}$, and find that it works better at higher masses, our results strongly suggest that $\epsilon_{R}$ is much better than $\delta_{R}$ for modeling smaller masses as well. Our $\epsilon$ESP approach provides an excellent description of the speeds of protohalo patches (equations 59 and Figure 9). It also provides a good description of halo abundances (Figure 11) provided we include the fact that there is some stochasticity in $\epsilon_{R}$ values, even at fixed mass (Figure 10 and equations 72 and 73). This motivates further study of our energy-based excursion set peaks model.

• •

An obvious direction for future work is to check if energy- rather than density-based halo finders are less stochastic in the late time field.

• •

A more ambitious goal is to determine $\epsilon_{c}$ from first principles, including the mean trends and stochasticity shown in Figures 8 and 10, and how these depend on the algorithm used to identify haloes in the simulations. This necessarily requires modeling the virialization process (e.g., using energy conservation arguments).

• •

A closely related effort is to extend the analysis to allow triaxiality. This ingredient is required because some (but not all!) of the stochasticity in protohalo matter overdensity correlates with the initial tidal shear field (Sheth et al., 2001; Hahn et al., 2009; Borzyszkowski et al., 2017). In models of triaxial collapse (e.g. Bond & Myers, 1996; Monaco, 1997; Angrick & Bartelmann, 2010; Ludlow et al., 2014), the leading order departure from sphericity is quadrupolar. In this respect, our energy-based approach is attractive because choosing locations that are stationary points of the energy sets the dipole to zero but leaves all other multipoles unconstrained. I.e., $\epsilon$ESP provides a natural setup for including deviations from sphericity. These effects would add to the stochasticity in $\delta_{R}$ discussed in this paper, which exists even in a purely spherical configuration.

• •

Another direction is to estimate the implications for the spatial distribution of the protohalo patches: halo ‘bias’. An excursion set-based model for halo abundances (like any analytical model) carries with it a description of halo bias (Sheth & Tormen, 1999; Musso et al., 2012; Castorina et al., 2017; Modi et al., 2017). However, protohalo patches appear to be better described by a smoothing window that oscillates less strongly than $W_{\rm TH}=W_{1}$ (Chan et al., 2017); our $W_{2}$ (equation 11) is indeed more damped than $W_{\rm TH}$ (equation 10), so it has qualitatively the right behaviour. Moreover, $W_{2}$ has zeroes on slightly larger scales $x=kR$ than $W_{\rm TH}$ (Figure 1), which is also qualitatively consistent with the measurements. This, with the modified velocity statistics and associated ‘velocity bias’ of equations (59), provides an additional test of our approach.

• •

Excursion sets also provide a natural framework for describing the phenomenon known as assembly bias: haloes of fixed mass may cluster differently depending on their assembly history (Sheth & Tormen, 2004; Gao et al., 2005; Faltenbacher & White, 2010; Lazeyras et al., 2017; Borzyszkowski et al., 2017). Models based on a smoothing filter that is a tophat in $k$-space (rather than $W_{1}$), have $\delta_{R}$ that is Markov: in effect, the memory of larger scales is absent in such models. For Markov models to exhibit assembly bias, other variables than $\delta_{R}$, such as the shear, must matter for halo formation (Keselman & Nusser, 2007; Castorina & Sheth, 2013). For density-based models that use a Top-Hat filter $W_{1}$, only $\delta_{R}$ and its slope can be used to explain assembly bias phenomena, because statistical moments with more powers of $k$ (like the peak curvature or the second $R$-derivative) diverge (Musso & Sheth, 2014), effectively making the slope of $\delta_{R}$ Markov. Gaussian smoothing allows density and second derivatives to be used (Zentner, 2007; Dalal et al., 2008), but because slope and curvature are the same for this filter, there is, in effect, only one additional parameter that can be used to model assembly bias phenomena. Moreover, the connection to the physics of spherical collapse is lost. In our energy-based approach – i.e. smoothing with $W_{2}$ – the second derivatives are statistically well behaved, and slope and curvature are distinct. This potentially allows physically realistic models of additional halo properties.

• •

The excursion set troughs model uses an overdensity-based approach to study abundance and clustering of voids in the large scale matter distribution (Sheth & van de Weygaert, 2004; Paranjape et al., 2012; Jennings et al., 2013; Achitouv et al., 2015; Massara & Sheth, 2018). It is natural to study if our energy-based approach is useful for voids and other constituents of the cosmic web. If so, incorporating it into the ‘skeleton’ framework (Hanami, 2001; Sousbie et al., 2008; Codis et al., 2015; Musso et al., 2018) or other approaches (Bond et al., 1996; Shen et al., 2006; Neyrinck, 2016), should be as straightforward as it was for peaks.

• •

Finally, our approach is relevant in the context of modeling primordial black hole abundances (Nakama et al., 2014; Germani & Musco, 2019) because there too, the condition for formation involves multiple variables for which the total Jacobian is no longer the simple product of the peak determinant times the upcrossing term, so the extra complications treated in Appendix C are necessary (Germani & Sheth, 2020; Young & Musso, 2020).

Acknowledgements

We are grateful to the ICTP for its hospitality over the years, but especially during the summer of 2014 when we first discussed the main idea presented here, and to the IFPU and the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe” for their hospitality during July 2019 when this manuscript was completed. We thank Giulia Despali for sharing the data of the Flora simulation with us, and Corentin Cadiou for many helpful discussions on numerical aspects.

Data availability statement

The data of the Flora simulation, and those underlying this article, can be shared on reasonable request to the corresponding author.

Appendix A Multipole expansion

The Fourier modes of the rescaled peculiar acceleration $-\nabla\phi$ are $(i\bm{k}/k^{2})\delta({\bm{k}})\mathrm{e}^{i\bm{k}\cdot\bm{r}}$. Those of its average over $V$ have $\mathrm{e}^{i\bm{k}\cdot\bm{r}}$ replaced by $3j_{1}(kR)/kR$. Those of $\bm{r}_{\mathrm{cm}}/3=\bm{D}_{R}R/3$ have it replaced by $j_{2}(kR)$. Since $j_{2}(x)-3j_{1}(x)/x=-j_{0}(x)$, equation (4) gives

[TABLE]

Expanding the plane wave in spherical harmonics returns

[TABLE]

where $Y_{i_{1}\dots i_{l}}(\bm{\hat{r}})=[r^{l+1}/(2l-1)!!]\partial^{l}(1/r)/(\partial r_{i_{1}}\dots\partial r_{i_{l}})$ is a modified version of the real spherical harmonics. Their indices are by construction totally symmetric and traceless. The first few of them are $Y=1$, $Y_{i}=\hat{r}_{i}$, $Y_{ij}=\hat{r}_{i}\hat{r}_{j}-\delta_{ij}/3$ and $Y_{ijk}=\hat{r}_{i}\hat{r}_{j}\hat{r}_{k}-(\delta_{ij}\hat{r}_{k}+\delta_{ik}\hat{r}_{j}+\delta_{jk}\hat{r}_{i})/5$.

For $r=R$, the $l=0$ term is just $j_{0}(kR)$, and drops out of equation (78). The $l=1$ term gives $-\hat{k}_{i}\hat{k}_{j}\delta({\bm{k}})W_{1}(kR)r_{j}$, and thus $-[(\delta_{ij}/3)\delta_{R}+q_{ij}]r_{j}$, where $q_{ij}$ is the traceless shear. In general, the $Y$’s obey $\hat{k}_{j}Y_{i_{1}\dots i_{l}}(\bm{\hat{k}})Y_{i_{1}\dots i_{l}}(\bm{\hat{r}})=Y_{ji_{1}\dots i_{l}}(\bm{\hat{k}})Y_{i_{1}\dots i_{l}}(\bm{\hat{r}})+[l/(2l+1)]Y_{i_{1}\dots i_{l-1}}(\bm{\hat{k}})Y_{ji_{1}\dots i_{l-1}}(\bm{\hat{r}})$. We can then write equation (78) as

[TABLE]

where

[TABLE]

Multiplying by $4\pi G\bar{\rho}$ gives back equation (13), with $M=(4\pi R^{3}/3)\bar{\rho}(1+\delta_{R})$, where $\delta_{R}$ is actually needed only in the first term and is of second order elsewhere.

Appendix B Correlations between variables

In what follows, we will drop the explicit dependence on $R$. We compute the cross-correlation coefficients of the normalized variables $\omega=\epsilon/\sigma_{02}$, $\nu=\delta/\sigma_{01}$, $x$ and $y$ defined in Section 3.1. Since

[TABLE]

then in the same order

[TABLE]

Notice that equations (90) and (91) also imply that

[TABLE]

Equation (31) implies $\mathrm{d}\epsilon/\mathrm{d}\ln R=-(R^{2}\sigma_{23}/7)y=5(\delta-\epsilon)$. Therefore

[TABLE]

so

[TABLE]

The final term in square brackets is a measure of the stochasticity between $\delta_{R}$ and $\epsilon_{R}$. The normalized correlation coefficient between the two is

[TABLE]

this is closer to unity when the large scale power dominates over that from smaller scales. Similarly,

[TABLE]

which is therefore equivalent to equation (48).

In addition, since $\langle(\mathrm{d}\epsilon/\mathrm{d}\ln R)^{2}\rangle=(R^{2}\sigma_{23}/7)^{2}$, then

[TABLE]

For power-law spectra, the term in square brackets is a constant which decreases as the slope of $P(k)$ decreases (becomes more negative). Equations (90) and (48) then give

[TABLE]

Figure 12 shows how a number of these cross-correlation coefficients for a $\Lambda$CDM power spectrum depend on smoothing scale $R$, parametrized by $\sigma_{02}$. To understand this scale dependence, it is useful to study the case in which $P(k)\propto k^{n}$. Then ${\rm d}\ln\sigma_{jl}^{2}/{\rm d}\ln R=-(2j+n+3)$ (independently of the filter), and, for $n=(0,-1,-2)$,

[TABLE]

where ND indicates that the coefficient contains divergent integrals.

Appendix C A threshold for $\delta_{R}$

We discuss here how our formalism should be modified if one wanted to introduce a critical value for the density $\delta_{R}$, while still fixing the position $\bm{x}$ through the constraint $\nabla\epsilon_{R}=0$. In this case, using the formalism introduced in Section 3.5, the number density of stationary points would be

[TABLE]

where $\tilde{J}=\det[\partial\{\delta_{R}-\delta_{c},\bm{\eta}\}/\partial\{R,\mathbf{x}\}]$ is the 4-dimensional Jacobian determinant. However, now $\nabla_{i}(\delta_{R}-\delta_{c})\neq 0$, since the peak constraint is enforced on $\epsilon_{R}$ and not on $\delta_{R}$, and $\tilde{J}$ no longer factorizes in a simple way. Assuming for simplicity that $\nabla_{i}\delta_{c}=0$, one has

[TABLE]

The total determinant is no longer the simple product of the peak determinant times the upcrossing term $\partial\delta/\partial R$.

We recall that the position $\bm{x}_{\mathrm{pk}}$ of the peak (or in general of the stationary point) depends on the smoothing scale $R$. Therefore, we introduce the “convective” derivative

[TABLE]

which accounts for this displacement. In order to preserve the constraint $\eta_{i}=0$ across scales, one must impose that $\mathrm{d}\eta_{i}/\mathrm{d}R=0$, which implies that

[TABLE]

Hence, the Jacobian determinant can be written as

[TABLE]

which describes the crossing along the trajectory that follows the stationary point of $\epsilon_{R}$. Moreover, since $\eta_{i}=\nabla_{i}\epsilon/\sigma_{12}$ and $\partial W_{2}/\partial R=(5/R)(W_{1}-W_{2})$, we also have that

[TABLE]

where the first term vanishes because $\eta_{i}=0$. Therefore, equation (110) becomes

[TABLE]

where we have used equations (26) and (27) to set $\sigma_{12}\nabla_{i}\eta_{j}=\sigma_{22}\,\zeta_{ij}$. Equation (114) shows that if $\bm{x}_{\mathrm{pk}}$ (a stationary point of $\epsilon_{R}$) is also a stationary point of $\delta_{R}$, i.e. if both $\bm{\nabla}\epsilon_{R}=0$ and $\bm{\nabla}\delta_{R}=0$, then convective and partial derivatives of any field are equal.

Since $\zeta_{ij}$ is positive definite due to the peak constraint, so is $\zeta_{ij}^{-1}$, and the second term from the above equation applied to $\delta$ is always positive. It can then happen that $\mathrm{d}\delta/\mathrm{d}R>0$ even if $\partial\delta/\partial R<0$. Hence, requesting upcrossing along the peak trajectory imposes the stricter constraint

[TABLE]

with $x_{R}$ as in Equation (30). If $\lambda_{1}>\lambda_{2}>\lambda_{3}>0$ are the eigenvalues of $\zeta_{ij}$, the inequality above is equal to

[TABLE]

where we denoted $x$, $y$ and $z$ the Cartesian coordinates in the frame in which $\zeta_{ij}$ is diagonal. This constraint reduces to the ordinary upcrossing condition $x_{R}>0$ if $\bm{\nabla}\delta=0$.