# Estimating major merger rates and spin parameters ab initio via the clustering of critical events

Corentin Cadiou, Eric Pichon-Pharabod, Christophe Pichon, Dmitri Pogosyan

► To cite this version:

Corentin Cadiou, Eric Pichon-Pharabod, Christophe Pichon, Dmitri Pogosyan. Estimating major merger rates and spin parameters ab initio via the clustering of critical events. Monthly Notices of the Royal Astronomical Society, 2024, 531 (1), pp.1385-1397. 10.1093/mnras/stae1220. hal-04246812

HAL Id: hal-04246812

https://hal.science/hal-04246812v1

Submitted on 26 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Estimating major merger rates and spin parameters ab initio via the clustering of critical events

Corentin Cadiou<sup>1,2,*</sup> Eric Pichon-Pharabod,<sup>3</sup> Christophe Pichon<sup>2,3,4</sup> and Dmitri Pogosyan<sup>5</sup>

Accepted 2024 May 3. Received 2024 April 15; in original form 2023 September 21

## ABSTRACT

We build a model to predict from first principles the properties of major mergers. We predict these from the coalescence of peaks and saddle points in the vicinity of a given larger peak, as one increases the smoothing scale in the initial linear density field as a proxy for cosmic time. To refine our results, we also ensure, using a suite of ~400 power-law Gaussian random fields smoothed at ~30 different scales, that the relevant peaks and saddles are topologically connected: they should belong to a persistent pair before coalescence. Our model allows us to (a) compute the probability distribution function of the satellite-merger separation in Lagrangian space: they peak at three times the smoothing scale; (b) predict the distribution of the number of mergers as a function of peak rarity: haloes typically undergo two major mergers (>1:10) per decade of mass growth; (c) recover that the typical spin brought by mergers: it is of the order of a few tens of per cent.

**Key words:** large-scale structure of Universe – cosmology: theory.

## 1 INTRODUCTION

On large scales, the galaxy distribution adopts a network-like structure, composed of walls, filaments, and superclusters (Geller & Huchra 1989; Sohn et al. 2023). This network is inherently tied to the cosmic microwave background, the relic of the density distribution in the primordial Universe. The non-uniformity of this initially quasi-Gaussian field evolved under the influence of gravity into the so-called cosmic web (Bond, Kofman & Pogosyan 1996) we now observe. One can therefore hope to predict the evolution of the cosmic web by studying the topological properties of the initial density field. From its evolution, one should be able to predict the rate of mergers of dark haloes and their geometry hence their contribution to halo spin.

The classical method to study mergers is to run cosmological simulations (e.g. Bertschinger 1998; Vogelsberger et al. 2020), compute where haloes are located at each time increment and construct that way their merger tree (e.g. Lacey & Cole 1993; Moster, Naab & White 2013).

The theory of merger trees for dark haloes has a long-standing history starting from the original Press-Schechter theory (Press & Schechter 1974), excursion set (Bardeen et al. 1986; Peacock & Heavens 1990; Bond et al. 1991), and peak patch theory (Bond & Myers 1996) or related formalisms (Manrique & Salvador-Sole 1995, 1996; Hanami 2001; Monaco, Theuns & Taffoni 2002; Salvador-Solé, Manrique & Botella 2022). One notable recent variation is the suggestion to use peaks of ‘energy’ field as such progenitors

(Musso & Sheth 2021). One of the prevailing theoretical ideas is to consider peaks of the initial density field in initial Lagrangian space, smoothed at scales related to halo masses, to be halo progenitors (Bardeen et al. 1986). The statistical properties of mergers can then be predicted analytically through extensions of the excursion set theories (Lacey & Cole 1993; Neistein & Dekel 2008), or, alternatively, it can be measured in peak-patch simulations (Stein, Alvarez & Bond 2019). In the first approach, all the information is local, preventing us from computing the geometry of mergers. In the second approach, the geometry of mergers is accessible as a Monte Carlo average.

In this paper, we provide an alternative framework that specifically takes into account the geometry of mergers while remaining as analytical as possible. The goal is not for accuracy, but rather to provide a simple framework to access information such as the statistics of the number of mergers or the spin brought by mergers.

The paper is organised as follows: in Section 2, we present our model together with analytical estimates of the number of direct mergers with a given halo In Section 3, we extend our model to take into account the topology of the density field using persistence pairing via DisPerse. This allows us to draw mock merger trees and to predict the probability density function (PDF) of the spin parameter that each merger contributes as a function of rarity. Section 4 wraps up.

Section A introduces the relevant spectral parameters. Section B revisits event statistics while relying on the clustering properties of peaks and events computed from first principles. Section D recalls the critical event PDF and provides a fit to it. Section E discusses the evolution of peak rarity with Gaussian smoothing.

<footnote n="1"><sup>1</sup>Lund Observatory, Division of Astrophysics, Department of Physics, Lund University, Box 43, SE-221 00 Lund, Sweden</footnote>

<footnote n="2"><sup>2</sup>CNRS and Sorbonne Université, UMR 7095, Institut d'Astrophysique de Paris, 98 bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>IPHT, DRF - INP, UMR 3680, CEA, Orme des Merisiers Bat 774, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="4"><sup>4</sup>Department for Physics, Korea Institute of Advanced Studies (KIAS), 85 Hoegiro, Dongdaemun-gu, Seoul 02455, Republic of Korea</footnote>

<footnote n="5"><sup>5</sup>Department of Physics, University of Alberta, 11322-89 Avenue, Edmonton, Alberta T6G 2G7, Canada</footnote>

<footnote>*E-mail: corentin.cadiou@astro.lu.se</footnote>
# 2 ANALYTICAL MODEL OF MERGERS PER HALO

## 2.1 Model of halo mergers: cone of influence of a halo

The starting point of our approach is the peak picture (Bardeen et al. 1986). In this framework, peaks in the linear density field are the seeds for the formation of haloes. We rely on the proposition that the overdensity $\delta(R)$ of a given peak smoothed at scale $R$ can be mapped to a formation redshift $z$, using the spherical collapse model, $\delta(R_{\text{TH}})D(z) = \delta_c$, where $\delta_c = 1.686$ is the critical density, $D(z)$ is the growing mode of perturbations, and $R_{\text{TH}}$ is the top-hat scale. The mass of the halo formed at this time is inferred from the top-hat window scale, $M = \bar{\rho}_m^{4/3}\pi R_{\text{TH}}^3$.

While the peak-patch theory provides useful insight into the origin of haloes of a given mass at a given time, it becomes singular during mergers. Indeed, by construction, the merging halo disappears into the larger one together with its corresponding peak. Here, we instead rely on the analysis of the geometry of the initial linear density field in (N + 1)D, where N is the dimension of the initial Lagrangian space and the extra dimension is the smoothing scale. At each smoothing scale, we can formally define critical points (peaks, saddles, and minima) following Bardeen et al. (1986). As smoothing increases, some of these critical points eventually disappear when their corresponding halo (for peaks) merge together. This approach allows to capture mergers in space-time from an analysis of the initial conditions in position-smoothing scale employing the critical event theory (Hanami 2001; Cadiou et al. 2020), where critical events of coalescence between peaks and saddle points serve as proxies for merger events. The process is sketched in Fig. 1 and illustrated in (1 + 1)D in Fig. 2. In the following, we will rely on the use of Gaussian filters, which allow us to advance further the analytical description. First, this choice allows us to obtain analytical results with the critical event theory. Second, it yields peaks whose density decreases monotonically with smoothing scale – and hence the collapse time of the associated haloes grows under the spherical collapse model. While exploring the effect of different filters is beyond the scope of this paper, we note that using any filter that is positive with sufficiently smoothly tapered boundaries does not pose any fundamental challenge to the numerical analysis we perform later on. Tophat filter by itself is not of this class, having sharp boundaries leading to ill-defined second- and third-derivatives of the smoothed field with cold dark matter-like power spectrum.

Let us track the Lagrangian history with decreasing smoothing of one peak first identified at a smoothing scale $R_0$ and position $x_{pk}(R_0)$. Since we employ Gaussian filters, we need to match the Gaussian and Top-Hat smoothing scales, $R$ and $R_{\text{TH}}$ to assign mass to haloes. The criteria of equal mass encompassed by the filter<fnref n="1" /> gives $R_{\text{TH}} \approx 1.56R$ in 3D. At $R_0$ the peak describes a halo that has collected its mass from a spherical cross-section of (N + 1)D space of volume $\propto R_0^N$; we can call this sphere a Lagrangian patch of the halo. At smaller $R$, the Lagrangian position of a peak changed to $x_{pk}(R)$ and the volume of its Lagrangian patch decreased to $\propto R^N$. The history of the peak in the (N + 1)D space, including mass accumulation, now consists of its trajectory $x_{pk}(R)$ and a cone of cross-section volume $\sim R^N$ around it as shown in Fig. 3 for (2 + 1)D example.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Typical density field around a critical event in 1D. From left to right: (1) We see two distinct objects (two peaks in red), separated by a minimum (in blue) in the density field. (2) The critical event occurs (inflection point in green), but the density profile still shows two different objects (highlighted with dashed lines). (3) The merger has completed and there remains only one object (peak in red).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Representation of the smoothing of a 1D Gaussian random field representing density fluctuations in the primordial universe. Red lines are maxima of the field, blue lines are the minima. Green points are critical events, where a minimum/maximum pair meets and annihilates itself. By counting how many critical events there are in the vicinity of a peak line, one can count the number of subhaloes merging into the corresponding halo.</figcaption></figure>

A critical event marks the end of a trajectory of a peak that disappeared when its scale reached $R_{ce}$ and is absorbed into a surviving peak. Counting all critical events within a (N + 1)D straight cylinder with a spatial cross-section of radius $R_0$ around the final surviving peak position $x_{pk}(R_0)$ will give the number of all mergers that ever happened within this peak Lagrangian patch. For instance, if two small haloes have merged together before merging with a larger halo, that would count as two events. However, we are interested in counting only the last direct merger event that brought the combined mass of the two small haloes into the main one. Physical intuition tells us that to count only those, we need to count critical events within the history cone of the peak mass accumulation. Indeed, the Lagrangian patch of the peak grows in size by absorbing the layers along its boundary, and if that layer contains a critical event, it is a direct merger.

In (1 + 1)D, as illustrated in Fig. 4, direct mergers correspond to critical events (green points) that are not separated by any intervening minimum (blue line) from the surviving peak. In (N + 1)D, this is generalized by only counting a critical event as a direct merger if, at fixed smoothing $R = R_{ce}$, it is connected to the main peak in ND space by a filamentary ridge with no other saddles in-between. Fig. 4 confirms that, indeed, one can find a 'cone of influence' around a peak line which spatial cross-section radius is $\alpha R$, where $\alpha$ is between 2 and 4, that contains most and, vice versa, almost exclusively the critical events that are directly connected to this peak. The very fact that such a cone can be defined, at a statistical precision, even in

<footnote n="1">¹Different criteria modify the relation between Gaussian and Top-Hat filters, for instance matching the variance of the perturbations leads to $R_{\text{TH}} \approx 2.1R$. Here, we instead rely on matching masses within the filtering window $W$, which implies $\int W_{\text{TH}}(r, R_{\text{TH}}) d^N r = \int W_G(r, R_G) d^N r$.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Visualization of the action of Gaussian smoothing on the critical events of a 2D field. The vertical axis is the smoothing scale, increasing upwards. The horizontal cross-section is a 2D space. The various vertical lines represent the tracks of extrema positions (maxima in pink, minima in blue, saddles in green) as one changes the field smoothing. The red and blue squares represent the corresponding critical events (in red are points of peak-saddle and in blue that of minima-saddle coalescence). The grey cones show the volume within some fraction of the smoothing scale (here chosen arbitrarily to be 1.2 times the smoothing scale) around each maxima track, which contain all the past physical history of a given peak. This paper aims to characterize these cones and the properties of the critical events within them so as to compute major merger rates as a function of final halo mass.</figcaption></figure>

principle, is non-trivial. The reason for that is, as Fig. 4 demonstrates again, the presence of the exclusion zone around the main peak devoid of any other critical events. This causes direct merger events to be most likely located in the first layer of critical events from the main peak. This exclusion effect is expected from peak-peak correlation studies (e.g. Baldauf et al. 2021) and can be used as a quantitative way to determine the boundary of the ‘cone of influence’ that we develop in Sections 2.3.2 and B.

Counting critical events within the ‘influence cone’ of the peak that contribute to the mass and spin growth of this surviving halo is the main analytical tool of this paper.

The evolution of a halo in smoothing direction tells us directly its history in terms of mass accumulation, as $M \propto R^N$. That is, we can describe what happened with the halo as its mass increased, say, 10-fold. In the next section, we apply this picture to count merger events. We shall limit ourselves to the (3 + 1)D case as it bears the most relevance to dark matter halo formation. We however note that there is no additional theoretical difficulty in deriving the general (N + 1)D case.

## 2.2 Number of major merger events within a mass range

Here and in the following, we are interested in counting the number of objects that directly merged into an object of final mass $M_0$ as its mass grew from $f^3 M_0$ to $M_0$, where $f < 1$. Since smoothing scale maps directly on to mass, this amounts to counting the number of critical events between two scales $R_0 \propto \sqrt[3]{M_0}$ and $R_1 \propto f \sqrt[3]{M_0}$.

Let us now count direct major mergers that we define as mergers with satellites that bring at least $f^3 M_0$ mass in the merger event. First, note that we define here the mass ratio with respect to the final mass of the peak at a fixed time, rather than at the time of the merger. The number of such mergers, as halo grew from $f^3 M_0$ to $M_0$, is given by the number of critical events in the section of the cone of influence of the halo contained between $fR_0$ and $R_0$. It can be obtained with the following integral

$$N_{\text{merger}} = \int_{f R_0}^{R_0} dR \int_{0}^{\alpha R} d^3 r n_{\text{ce}}(R, r), \quad (1)$$

where $n_{ce}(R, r)$ is the number density of critical events at the point $(R, r)$ in the extended (3 + 1)D space of positions-smoothing scale and the adjustable parameter $\alpha$, introduced in the previous section, sets the radial extent of the $R = \text{const}$ cross-section of the influence cone. Note that in the above formula, we are able to avoid the dependence on the past trajectory $x_{pk}(R)$ of the halo, by treating each slab of constant $R$ independently, and by evaluating the radial distance $r$ to the critical event found of at $R$ from the main peak position, $x_{pk}(R)$, defined at the same smoothing.

As a first estimate, let us approximate $N_{merger}$ by taking the density of the critical events inside the cone of influence to be equal to its

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Our model associates peaks (in red) to haloes, and critical events (in green) to mergers as a function of smoothing scale $R$ which encodes mass evolution. Critical events correspond to the coalescence of a peak with a minimum (in blue). We count any critical event directly connected to the central peak (i.e. not separated by another minimum) as a merger. Here, this definition would give six mergers (shown with grey dashed lines) into the central peak. Alternatively, we can estimate this number by counting the number of critical events within the shaded grey ‘influence cone’ of spatial width $\alpha R$. Here, $\alpha = 3.3$ is chosen so that the cone includes only directly connected events capturing five events. Note that the left-most event indicated by an arrow is missed, but it lies at large smoothing $R$ where the sample variance is large. The second arrow gives an example of an event that is not directly connected to the central peak, instead merging with the peak on its left. The cone boundary separates this event from the directly connected one above it.</figcaption></figure>
mean value. Thus, we first neglect any spatial correlation between the existence of a peak (corresponding to the surviving object in a merger) and the critical event (corresponding to the absorbed object in a merger). In Cadiou et al. (2020), we determined that the average density of critical events in (3 + 1)D that correspond to peak mergers is given by

$$ \bar{n}_{ce} = \frac{R}{R_*^5} \frac{1 - \tilde{\gamma}^2}{\tilde{\gamma}^2} \frac{29\sqrt{15} - 18\sqrt{10}}{600\pi^2}, \quad (2) $$

where the spectral scale $R_*$ and parameter $\tilde{\gamma}$ are given in Section A, together with other parameters that characterize the statistics of a density field. Note that $R_* \propto R$, so the number density of critical events scales as $n_{ce} \propto R^{-4}$. The cubic part of this dependence, $R^{-3}$, reflects the decrease of spatial density of critical points with increasing smoothing scale; the additional scaling, $R^{-1}$, reflects that the frequency of critical events is uniform in $\log R$.

Assuming a power-law density spectrum with index $n_s$ and Gaussian filter, equation (2) gives

$$ \bar{n}_{ce} \approx 0.0094 \left( \frac{5+n_s}{2} \right)^{3/2} R^{-4}. \quad (3) $$

Using this mean value $\bar{n}_{ce}$ for critical density in equation (1) gives us a first rough, but telling estimate of the merger number that a typical halo experiences

$$ N_{merger} \approx 0.039(-\ln f) \left(\frac{5+n_s}{2}\right)^{3/2} \alpha^3. \quad (4) $$

The final results depend on the specific values of $f$ and $\alpha$. The value of the fraction $f$ controls down to which mass ratio mergers should be included and has thus a straightforward physical interpretation. The parameter $\alpha$ controls the opening of the cone of influence around the peak from which critical events should be considered as direct mergers and is to be determined.

Using equation (4) with $n_s = -2$, which corresponds to the typical slope of a $\Lambda$ cold dark matter ($\Lambda$CDM) power-spectrum at scales ranging from Milky-Way-like systems to clusters, we can count the number of direct mergers with mass ratio larger than 1:10 ($f = 1/\sqrt{10}$). For the choice of $\alpha$, a natural idea would be to count only critical events within the halo Top-Hat Lagrangian radius, then $\alpha = R_{TH}/R \approx 1.56$ and we find $N_{merger} \approx 0.2$. However, this is clearly an underestimation for $\alpha$, as it corresponds to the situation where the centre of the merging halo has already been incorporated into the main one when the latter was of equal mass. A more sensible choice, as will bear out in the analysis presented in the next section is to extend the opening of the cone to twice that ratio $\alpha = 2R_{TH}/R \approx 3.1$ which corresponds to the two mass spheres of the surviving and merging peaks touching at the scale of the critical event, signifying the merger onset. We then find $N_{merger} \approx 1.7$. Overall, we see that the number of direct major mergers that a halo experiences while increasing its mass tenfold is small, of the order of one or two. For scale invariant history, each decade of mass accumulation contributes a similar number of mergers, so a cluster that grew from a galactic-scale protocluster thousand-fold in mass ($f = 1/\sqrt{1000}$), did it having experienced $N_{merger} \approx 5$ major mergers. These conclusions follow directly from first principles while studying the structure of the initial density field.

## 2.3 Accounting for rarity and clustering

Let us now refine our model so as to define $\alpha$ more rigorously by (i) requiring the merging object to have gravitationally collapsed before

<figure><img src="image_6.png" /><figcaption>Figure 5. Panel (a) shows the tracks of two peaks (solid line, red) and a saddle point (dashed line, blue) of the density field as a function of smoothing scale $R$ in 1 + 1D space. Panel (b) shows the field profile at three smoothing scales. At smoothing scale $R_0$, the peaks and saddle points are distinct. At smoothing $R_{ce}$, one peak and a saddle point create a critical event C, after each only one peak survives to larger smoothing scales. The merging of peaks is completed at the scale $R_{pk}$ when the overdensity of the surviving peak, now at point A, is equal to the overdensity of the critical event $\delta_c$ and thus can be viewed as reaching the threshold $\delta_c$ for halo formation at the same time. We can then interpret the critical event as a merger with mass ratio $M_A/M_C = (R_{ce}/R_{pk})^3$ at a redshift corresponding to $D(z) = \delta_c/\delta_A$.</figcaption></figure>

it merges and by (ii) taking into account the correlations between the central halo peak and the merger critical events.

Let us therefore consider the number density of critical events of given height $v$ as a function of their distance $r$ to a central peak, both defined at the same smoothing scale $R$

$$ n_{ce|pk}(v, R, r) = \bar{n}_{ce}(R)C(v, \gamma) (1 + \xi_{ce,pk}(v, R, r)), \quad (5) $$

where $C$ is the distribution function, $\int_{-\infty}^{\infty} C(v, \gamma) dv = 1$, of critical events in overdensity and is given by the analytical formula in Section D along with a useful Gaussian approximation. The clustering of critical events in the peak neighbourhood is described by the peak–critical event correlation function on a slice of fixed $R$, $\xi_{ce,pk}(v, R, r)$. The composite index pk refers to any peak parameters that may be specified as a condition. Our goal is to determine what range of $v$ and what extent of $r$ one needs to consider to count haloes merging into a surviving halo with a particular $v_{pk}$.

### 2.3.1 Heights of critical events

To establish the range of critical event heights that describe mergers of real haloes, we rely on the spherical collapse approximation to map the peak overdensity to its gravitational collapse time. We consider a physical halo at redshift $z$ to be described by a peak at the Top-Hat scale $R_{pk}$ such that $\delta_{pk}(R_{pk}) = \delta_c/D(z)$, where $\delta_c$ is the critical overdensity for collapse model and $D(z)$ is the linear growing mode value at redshift $z$. Thus, a critical event found at smoothing $R_{ce}$ describes the merger of a satellite halo of scale $R_{ce}$ into the main peak at scale $R_{pk}$ that corresponds to the same redshift, i.e. such that $\delta_{ce}(R_{ce}) = \delta_{pk}(R_{pk}) = \delta_c/D(z)$. The situation is demonstrated in Fig. 5. The ratio of scales and, correspondingly, masses of the two merging haloes is therefore determined by the condition of equal overdensities at the time of the merger.

Let us now consider a peak that has reached a scale $R_0$ and that experienced a merger in its past when its scale was $R_{pk}$, $R_{pk} \le R_0$. Requiring that, at the time of the merger, the surviving peak’s scale (mass) is larger than that of the satellite sets the relation $R_0 \ge R_{pk} \ge R_{ce}$ or, conversely, $\sigma_0(R_0) \le \sigma_0(R_{pk}) \le \sigma_0(R_{ce})$, where $\sigma_0$ is the root mean square (RMS) of the density field and is defined in Section A.
From which it follows that the rarity of the relevant critical events at scale $R_{ce}$,

$$v = \frac{\delta_{ce}(R_{ce})}{\sigma_{0}(R_{ce})} = \frac{\delta_{pk}(R_{pk})}{\sigma_{0}(R_{ce})} = v_{pk} \frac{\sigma_{0}(R_{pk})}{\sigma_{0}(R_{ce})} \quad (6)$$

is within the range

$$v_{pk} \frac{\sigma_{0}(R_{0})}{\sigma_{0}(R_{ce})} \le v \le v_{pk}. \quad (7)$$

The lower bound corresponds to mergers that have been completed at the very last moment, $R_{pk} = R_0$. The upper bound is achieved for mergers of two equal mass objects, $R_{pk} = R_{ce}$.

To obtain the total number of merger events in a halo history we now integrate the conditional event density in equation (5) over the physically relevant range of heights from equation (7), i.e.

$$N_{\text{merger}}(v_{\text{pk}}) = \int_{f_{R_0}}^{R_0} dR \int_{0}^{\infty} 4\pi r^2 dr \times \int_{v_{\text{pk}}}^{v_{\text{pk}}} \frac{\sigma_0(R_0)}{\sigma_0(R)} dv n_{\text{ce|pk}}(v, R, r). \quad (8)$$

### 2.3.2 Clustering of critical events around peaks

Ultimately, the conditional event density $n_{ce|pk}$ should include only critical events that will directly merge with the peak; this would make the density of critical events go to zero $n_{ce|pk}(v, R, r) \to 0$ far from the peak, i.e. when $r \to \infty$. While we cannot implement this condition analytically (as it is non-local), we can measure $n_{ce|pk}(v, R, r)$ numerically, as will be done in the upcoming Section 3. We can however approximate the conditional density *ab initio* by relaxing the conditions that the peak is linked to the critical event to obtain $\tilde{n}_{ce|pk}$ given by

$$\tilde{n}_{ce|pk}(v, R, r) \equiv \frac{\langle \text{Peak}(x) \text{ Event}(y) \rangle}{\langle \text{Peak}(x) \rangle}, \quad (9)$$

using formally straightforward analytical calculations of critical event–peak correlations, as described in Section B. Here, brackets denote an ensemble average, where $x$ and $y$ are the random vectors containing the density and its successive derivatives at the location of the peak and the critical event, respectively and ‘Peak’ and ‘Event’ enforce the peak and critical event conditions, respectively. We can expect $\tilde{n}_{ce|pk}$ to track the exact $n_{ce|pk}$ up to several smoothing length distances from the peak, but further away from the peak it just describes the mean unconstrained density of critical events: $\tilde{n}_{ce|pk}(v, R, r) \to \tilde{n}_{ce}(R)C(v)$ as $r \to \infty$. Therefore, in this approximation, the question of where to truncate the integration over the peak neighbourhood remains and we have

$$N_{\text{merger}}(v_{\text{pk}}) \approx \int_{f_{R_0}}^{R_0} dR \int_{0}^{\alpha R} 4\pi r^2 dr \int_{v_{\text{pk}}}^{v_{\text{pk}}} \frac{\sigma_0(R_0)}{\sigma_0(R)} dv \tilde{n}_{\text{ce|pk}}(v, R, r). \quad (10)$$

For scale-free spectra, introducing the dimensionless ratios $u = r/R$ and $w = R/R_0$ and changing the order of integration, equation (8) can be written as

$$N_{\text{merger}}(v_{\text{pk}}) = \int_{f}^{1} \frac{dw}{w} \int_{0}^{\infty} du \frac{d^2 N_{\text{ce|pk}}}{du \, d \ln w}, \\ \frac{d^2 N_{\text{ce|pk}}}{du \, d \ln w} = 4\pi u^2 \int_{v_{\text{pk}}}^{v_{\text{pk}}} dv w^4 n_{\text{ce|pk}}(v, w, u), \quad (11)$$

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Integration of the analytical expression of the clustering of critical events around peaks, equation (11), for different values of $v_{pk}$. Here $v_{ce}$ is integrated from $0.681 v_{pk}$ to $v_{pk}$, with power spectrum $n = -2$ and smoothing ratio $R_{pk}/R_{ce} = 0.95$. We detail the method used to compute this graph in Section C1. Rarer peaks have a stronger excess of events in their vicinity, which mirrors the known Kaiser bias (Kaiser 1984).</figcaption></figure>

where we note that the differential event count per logarithm of smoothing scales has only a dependence on the smoothing scale $w$ via the bounds of the height integration, since for power-law spectra $w^4 n_{ce|pk}(w)$ is $w$-independent.

In Fig. 6, we plot the differential version of equation (11) per radial spatial shell, integrated over $v$ in the range $0.681 v_{pk} \le v \le v_{pk}$ which corresponds to $(R/R_0)^3 = 1/10$ for $n = -2$ power-law spectrum. Fig. 6 shows that critical events are indeed preferentially clustered at distances $\sim 3R$ from a peak. It is natural to interpret this distance as the boundary of the cone influence of the peak, with more distant critical events constituting ‘the field’ of events that are not directly merging into the peak. Thus, the computation of the correlation function supports choosing $\alpha \approx 3$ in the estimate of equation (4), i.e. qualitatively double the value of $R_{TH}$. Note also that the probability of a critical event to be found close to a peak is suppressed, which corresponds to an exclusion from the interior of the peak’s cone of influence, as was qualitatively observed in Fig. 4.

## 3 NUMERICAL INTEGRATION OF MERGERS

We will now rely on a numerical approach to study the properties of the mergers a given halo undergoes through the properties of critical events that are absorbed by the peak representing the halo. The main numerical step is the identification of critical event–peak pairs that mark specific mergers. This in turn involves the identification in the initial field at a sequence of smoothing scales of connected peak–saddle–peak triplets that describe filamentary links between the merging peaks. At a particular smoothing scale, the saddle and one of the peaks of a given triplet merge into a critical event, while the remaining peak is a surviving halo. Such triplets are identified with DisPerse (Sousbie, Colombi & Pichon 2009) that performs Morse complex analysis of the random density field. We perform the analysis on an ensemble of scale-free Gaussian realizations of initial density fields.

### 3.1 Number of mergers per peak line

We identify proto-events as close-to-zero persistence-pairs of peaks and filament-saddles, and we study here their clustering relative
<figure><img src="image_8.png" /><figcaption>**Figure 7.** The radial number density of topologically linked critical event per peak as a function of distance to the peak on a slice of constant *R* and per logarithmic bin in smoothing scale (see the text for details of how critical events are identified and paired to a peak). We represent the mean separation with vertical lines. Critical events are preferentially clustered at *r* ~ 2.5 – 3.5*R*, depending on rarity. When compared to Fig. 6, the rarer peaks have relatively fewer connected events, suggesting that their satellites would have merged together before merging into the central object.</figcaption></figure>

to a given peak. We proceed as follows: we generate 430 scale-
invariant gaussian random fields of size 256³ with a power-law
power spectrum with spectral index *n* = -2. We smooth each
field with a Gaussian filter over multiples of a given smoothing
scale. Here we choose *R* = 1.1*^k* pixel with *k* running from 15 to
49. Each smoothing step corresponds to 33 per cent increase in the
mass associated with the smoothing scale, and eight steps give the
mass increasing approximately 10-fold. The skeleton is identified at
each scale using DisPerse (Sousbie et al. 2009) with a persistence
threshold of 1 per cent of the RMS of the smoothed fields. We
also measure at each critical point the value of the gravitational
acceleration. Since we rely on Gaussian filters, the overwhelming
majority of critical events encode the destruction of a peak with a
saddle with increasing smoothing, rather than the creation of such a
pair (the latter is 30 × less likely, see Cadiou et al. 2020, fig. F1). In the
following, we ignore the latter rare case as it has no straightforward
physical interpretation and is subdominant. We thus detect events by
following each pair of peak and saddle from one smoothing scale to
the next one. Each pair either survives at the next smoothing scale or
disappears at a critical event. In the latter case, we define the position
of the critical event as the middle point of the pair at the largest
smoothing scale where it survived. Given that DisPerse (Sousbie
2011) stores for each saddle the two maxima it is paired with, we
can therefore associate the critical event to the one peak that is not
involved in the event.

Using this data, we can now study the clustering of critical events
around peaks. In Fig. 7, we show the ratio of the number of critical
events per logarithmic smoothing scale bin per unit distance to the
number of peaks at the same scale and at the same rarity.<fnref n="2" /> Critical
events are counted in the interval of rarity 0.681 $v_{pk} \le v_{ce} \le v_{pk}$. This
range corresponds to the smallest scale cross-section of the influence
cone if we count major mergers of a peak that grows ten times in the
process, i.e. $R_{ce} = 1/\sqrt{10}R_0$ in equation (8). Note that our quantity
differs from two-point correlation functions in which one would
compute the distance between any critical event and any peak, as was

<table><caption>**Table 1.** Number of critical events per peak, both being measured at the same scale, per logarithmic bin of smoothing scale, as a function of the rarity of the peak.</caption>
  <thead>
    <tr>
      <th>v<sub>pk</sub> range</th>
      <th>[1.5, 2.0)</th>
      <th>[2.0, 2.5)</th>
      <th>[2.5, 3.0)</th>
      <th>[3.0, 3.5)</th>
      <th>[3.5, &infin;)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>dN<sub>ce|pk</sub></th>
      <td>0.90</td>
      <td>1.20</td>
      <td>1.40</td>
      <td>1.42</td>
      <td>1.27</td>
    </tr>
    <tr>
      <th>dln R</th>
      <td></td>
      <td></td>
      <td></td>
      <td></td>
      <td></td>
    </tr>
  </tbody>
</table>

done in the previous section (Fig. 6). Here, each critical event only
contributes to the density at a single distance to its associated peak,
and we thus expect the signal to drop to zero at infinite separation
since the probability for a critical event to be associated with a peak
far away vanishes.

Should the critical events be randomly distributed, their number
per linear *dr* would grow as *r*<fnref n="2" />. Instead, in Fig. 7 we find an exclusion
zone at small separations, an excess probability at *r*/*R* ~ 2–3, and a
cutoff at large separations as a consequence of our requirement for
the critical event to be paired to a peak.

We integrate the curves in Fig. 7 and report the mean number of
critical events per peak per logarithmic smoothing bin in Table 1.
All the effects accounted for in our numerical analysis give just
≈40 per cent lower values in Table 1 relative to equation (4) (per
lnf) with α = 3.1. This shows that the two main corrections to the
naive uniform density estimate – the restriction to only collapsed
satellite haloes (Section 2.3.1) on the one hand, and an attraction
of critical events of similar heights towards the peak influence zone
(Section 2.3.2) on the other hand – compensate each other to a large
extent.

As a function of peak rarity $v_{pk}$, the mean number of critical
events per peak in a constant $R$ slab first increases up to $v_{pk} \approx 3$
before decreasing. It is fairly flat in the range $v_{pk} = 2.5-3.5$, where
most of the physical haloes are, which argues for taking $\alpha$ parameter
independent on $v_{pk}$ if one uses the rough estimate as the global
mean density of events times effective volume as in equation (4),
Section 2.2.

So far, we have only studied properties of peaks and critical events
at the same scale, using the defining equation (7) to implicitly perform
a multiscale analysis. We can however track peaks from one slab of
smoothing scale to another to build peak lines and associate them to
critical events. This generalizes the procedure sketched on Fig. 4 in
3 + 1D. For each peak with mass $M_{pk}$, we follow its peak line
to find all associated critical events with a mass $M_{pk} \ge M_{ce} \ge
M_{pk}/10$ that have $\delta_{ce}(R_{ce}) \ge \delta_{pk}(R_{pk})$ where $R_{ce}$ and $R_{pk}$ are now
different. We only retain peaks that exist at scales $R_{pk} \ge R_{min}\sqrt{10}$
to ensure that we do not miss critical events below our smallest
scale.

Given our numerical sample, we are now in a position to study the
distribution of the number of mergers per peak, which we show on
Fig. 8. This measurement differs from the value one would obtain
by taking Table 1 multiplied by the logarithm of the range of scales.
Indeed, we account here for the decrease in the number of relevant
satellites as their mass approaches the final peak mass. We also
obtain our measurement by computing the number of mergers per
surviving peak at the scale $M_{pk}$. Compared to that, in Table 1, we
give the value per any peak at $1/10 M_{pk}$ scale, irrespective of whether
and how long it would survive at the further smoothing scales.
We find that the dependence of the number of mergers with peak
height is weak, with a mean number of mergers varying between
1.8 and 2. This shows that in the language of effective influence
cone radius of equation (4), $\alpha = 2R_{TH}/R = 3.1$ is a very good
proposition.

<footnote n="2">²Note that some peaks will not be associated to any critical event in the fixed interval of smoothing scales.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Distribution of the number of major mergers for peaks in different rarity bins, as labelled. The mean number of major mergers increases from 1.8 for low-v peaks to 2 for high-v peaks and is graphically represented with vertical lines with tilted caps for readability.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Merging objects bring in orbital spin; we quantify here its distribution. The distribution resembles a lognormal distribution with parameters $\mu_{\lambda} \approx 0.048$ and $\sigma_{\lambda} \approx 0.51$; the distribution resembles its N-body counterpart.</figcaption></figure>

## 3.2 Orbital spin parameter of mergers

Let us now define the spin parameter of an event at scale $R_{ce}$ relative to a peak of rarity $v$ at scale $R_{pk}$ as

$$ \lambda = \frac{R_{ce}^3 |\mathbf{r} \times (\nabla\psi - \nabla\psi_{pk})|}{R_{pk}^3 \sqrt{2} R_{pk} \sigma_{-1} (R_{pk})} \quad (12) $$

Here, we recall that $\sigma_{-1}$ is the variance of the gradient of the potential, whose definition is given in Section A. $\nabla\psi$ and $\nabla\psi_{pk}$ are the gravitational accelerations at the locus of the critical event and of the peak respectively, and $\mathbf{r}$ is the position of the critical event with respect to the peak. This definition reflects the spin definition of Bullock et al. (2001): under the Zel'dovich approximation, $R_{ce}^3 |\mathbf{r} \times (\nabla\psi - \nabla\psi_{pk})|$ is proportional to the mass times the cross product of the position with the velocity of the merging object relative to the peak at the time of its collapse, i.e. the numerator is proportional to the orbital angular momentum of the merger. Conversely, the denominator is proportional to the orbital angular momentum a merger of mass $\propto R_{pk}^3$ coming from a distance of $R_{pk}$ with velocity equal to the RMS of the velocity $\sqrt{2}\sigma_{-1}$ would have at the time of the collapse of the peak. We show on Fig. 9 the distribution of the spin parameter $\lambda$. The sample consists of mergers with mass ratios greater or equal to 1:10 around rare peaks. In practice, we select $v_{pk} > 2.5$; $R_{pk}/\sqrt{10} \le R_{ce} \le R_{pk}$ and $\delta_{ce} \ge \delta_{pk}$. The distribution is

roughly lognormal with a mean value of $\mu_{\lambda} \approx 0.048$ and a standard deviation of $\sigma_{\lambda} \approx 0.51$.

The novelty of our approach is to model mergers as punctual events. This has to be compared to previous theoretical works which treated angular momentum accretion around peaks as a continuous accretion process (e.g. White 1984; Ryden 1988). This allows direct comparisons to results about mergers obtained from N-body simulations. Remarkably, our estimate for the orbital spin of mergers is found to be very close to the values measured in N-body simulations, which are on the order of 0.04 (see e.g. Bullock et al. 2001; Aubert, Pichon & Colombi 2004; Danovich et al. 2015). This translates two effects. First, the mass ratio ($R_{ce}^3/R_{pk}^3$) is dominated by the more numerous smaller mergers, which drives $\lambda$ down. In addition, we find that critical events tend to have a mostly radial acceleration so that their tangential velocity is smaller than the RMS of the velocity field ($(\nabla\psi - \nabla\psi_{pk})/\sigma_{-1}$). This is qualitatively in line with findings from N-body simulations (Wetzel 2011; Jiang et al. 2015), where the tangential velocity of major mergers was found to be smaller than their radial velocity. Albeit simplistic, our model thus allows us to provide a natural explanation for the fact that mergers bring in a comparable amount of angular momentum to that of the full halo: gravitational tides alone (which is the only ingredient of our model) can funnel in a significant amount of angular momentum through mergers. This provides a theoretical motivation for the amplitude of the spin jumps during mergers employed in semi-analytical models (see e.g. Vitvitska et al. 2002; Benson, Behrens & Lu 2020).

## 3.3 Mass distribution of mergers

Finally, we compute the mass distribution of mergers as well as their density. The goal here is to obtain the distribution of the time and mass ratios of the mergers. In order to build the distribution of the mass ratio, we track peaks over multiple decades of smoothing scales (i.e. mass). Since the number density of peaks evolves as $R_{pk}^{-3}$, the sample size quickly decreases with smoothing scale; we select here peaks that exist at a scale larger than $R_{min}/\sqrt{100} \approx 4.6R_{min}$. The practical consequence is that our sample is only complete over two decades in mass ratio. We also only retain rare peaks that have $v_{pk} > 2.5$. Let us then estimate the time of the merger as follows: we associate both the peak and the critical event a time $t_{pk}$ and $t_{ce}$, respectively, using $\delta = \delta_c/D(z)$. Note that some of the selected peaks will have collapsed by $z=0$ while others will in the distant future. To aggregate peaks collapsing at different times, we compute the lookback time of the merger relative to the collapse time of the peak $\Delta = (t_{pk} - t_{ce})/t_{pk}$. We estimate the distribution of mergers in lookback time–mass ratio space using a 2D kernel density estimation which we show in Fig. 10, top panel. As expected, the more massive the merger, the more recently it happened. In the figure, we included mergers with mass ratios smaller than 0.01 in the shaded region rather than truncating the plot, but we remind the reader that our data set is not complete in this region.

We then show on the bottom panel of Fig. 10 the cumulative distribution function of the merger time, for mass ratios larger than 1:10 as a function of peak height. We recover the trend found in N-body simulations that rarer haloes have had more recent major mergers than lower mass ones. Cluster-like structures ($v_{pk} \gtrsim 3.5$) typically had 80 per cent of their mergers in the second half of their life (past 7 Gyr for a cluster at $z=0$), and had half their mergers in the last third of their life (last 5 Gyr for a cluster at $z=0$).

Our model reproduces qualitative trends observed in N-body simulations (figure 9 Genel et al. 2009; Fakhouri, Ma & Boylan-
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Top: distribution of the merger as a function of mass ratio and lookback time. We show the median with 68 per cent interval in black. The hashed area corresponds to regions of the parameter space that may not be complete, see the text for details. Bottom: the corresponding fraction of major merger as a function of lookback time relative to the peak collapse time, for different peak heights. Rarer objects have had more recent mergers.</figcaption></figure>

Kolchin 2010), namely that rarer peaks typically had their last major merger more recently than less rare ones.

Our model differs from the Extended Press–Schechter (EPS) approach (Bond et al. 1991; Bower 1991; Lacey & Cole 1993). In the EPS approach, random walks in $\sigma(R) - \delta(R)$ are mapped on to mass and time evolution, respectively. Jumps from one $R$ to another at fixed overdensity thus correspond to sudden changes in mass at a fixed time, i.e. mergers. While this approach yields statistically useful merger rates, it is unable to provide spatial information about mergers since neither the satellite nor the central halo are localized. Our model trades some of the analytical tractability of the EPS approach for additional information about the spatial location of mergers, all the while yielding a good agreement with $N$-body simulations.

# 4 CONCLUSION AND PERSPECTIVES

We built a model of mergers from an analysis of the initial conditions of the Universe. Following the work of Hanami (2001) and Cadiou et al. (2020), we relied on the clustering of critical events around peaks in the initial density field to study the properties of halo mergers. We started with a simple model that yields analytically tractable results, while further refinements presented in Section 3 allowed for more precise results at the cost of numerical integration.

We focused here on the analysis of merger events that bring at least 10 per cent of the final mass of the main halo, which we refer to as ‘major mergers’. We however note that our approach can be extended to the analysis of minor mergers, which we leave for future work.

We first obtained a zeroth-order analytical estimate of the mean number of mergers per decade in mass using the mean abundance

of critical event per peak in Sections 2.1 and 2.2. Our results are consistent with haloes having had one to two major mergers per relative decade of mass growth. We then refined our model in Section 2.3 by accounting for the timing of the collapse of the haloes involved in a merger candidate; a critical event should only be counted as a merger if its two associated haloes have collapsed before the merger happens. We showed this can be achieved semi-analytically by numerically computing the value of a cross-correlation function, equation (10), that reveals that critical events cluster at 2–3 times the smoothing scale of the peak (Fig. 6).

Finally, in Section 3, we addressed the double-counting issue, whereupon a given critical event may be associated with several peaks, by uniquely associating each to the one peak it is topologically connected to. To that end, we relied on multiscale analysis of Gaussian random fields using computational topology to restrict ourselves to the study of peaks that, up to the critical event, form persistent pairs. In this model, we found again that haloes of different rarities undergo about two major mergers. By tracking peaks in position-smoothing scale and by associating critical events with them, we were able to provide numerically cheap and easy-to-interpret data on the statistical properties of halo mergers. We found that mergers come from further away for rarer peaks (Fig. 7), but that the total number of major mergers only weakly depends on peak rarity (Fig. 8).

We then computed the gravitational tides at the location of the critical event to estimate the relative velocity of mergers and predict the orbital spin they bring in. We find that it has a lognormal distribution with a mean of $\mu_{\lambda} = 0.048$ and $\sigma_{\lambda} = 0.5$. These properties are remarkably close to the distribution of dark matter halo spins measured in hydrodynamical simulations ($\mu_{\lambda} = 0.038$, $\sigma_{\lambda} = 0.5$, Danovich et al. 2015). This suggests that our model captures the (statistical properties) of the orbital parameters of mergers, as is expected should they be driven by gravitational tides (Cadiou, Pontzen & Peiris 2021a, 2022). We also computed the distribution of the mass brought by mergers and their timing, which we found to be in qualitative agreement with results obtained in $N$-body simulations (Fakhouri et al. 2010).

While the aim of this model was not to compete with numerical simulations, it provides theoretical grounds to explain the properties of mergers observed in $N$-body simulations and efficient tools to predict their statistics and geometry *ab initio*. Our model could be improved with precision in mind, for example by taking into account deviations from spherical collapse under the effect of shears to improve our time assignments. It however reveals that the statistical properties of the merger tree of dark matter halo can be explained through a multiscale analysis of these initial conditions.

## 4.1 Perspectives

Statistics involving successive mergers could potentially be built on top of our model by using critical events associated with the same peak line, for example, to study the relative orientation of the orbital angular momentum of successive mergers. However, we found that such analysis was complicated by the fact that peaks move with smoothing scale. Different definitions of angular momentum (distance to the peak at the same scale, at the same density, or for a fixed peak density) yielded qualitatively different results. This should be explored in future work.

The model built in this paper relied on a linear multiscale analysis of the density field. This could be employed to provide control over the merger tree (timing of the merger, orientation) in numerical simulations through ‘genetic modifications’ of the initial field (Roth,
Pontzen & Peiris 2016; Rey & Pontzen 2018; Stopyra et al. 2021; Cadiou et al. 2021a, b). We also note that, as we did for tidal torque theory in Cadiou et al. (2022), such an approach would allow direct testing of the range of validity of the model.

In our analysis throughout the paper, we used Gaussian filters which provide closed analytical formulas for the critical event theory and also ensure that, along the peak trajectory, larger smoothing scales correspond to later redshifts of the collapse. In addition, Gaussian filters limit the extent to which peak-saddle point pairs can be created (rather than destroyed, see Cadiou et al. 2020), for which we have no physical interpretation yet. However, other choices of filter could be considered, as long as the filter is positive, so that the region associated with the peak at a larger smoothing incorporates the regions assigned with smaller smoothing, and has sufficiently smoothly tempered boundaries to define second and third derivatives of the field. For instance, the widely used Top-Hat filter has a compact support and allows to directly map the density to a collapse time and the smoothing scale to a mass through the spherical collapse model. But it has sharp boundaries and ill-defined derivatives beyond the first one for a density field with typical cosmological power spectra in ΛCDM-like hierarchical models. For critical event theory, one could consider a modification of Top-Hat that retains a compact support with nearly equal weight but has at least second derivatives continuous at the boundary. Closed analytical formulas for the critical event theory with such non-Gaussian filters are not known, however, we do not see any obstacle to numerical analysis using them. The focus of this work was to provide a theoretical framework to understand the properties of mergers, and we expect the qualitative results to hold for other reasonable filters.

This paper focused on mergers of peaks corresponding to the relative clustering of peak-saddle events. One could extend the analysis to the relative clustering of saddle-saddle events to provide a theoretical explanation for which filaments merge with which, thus impacting their connectivity or their length (Galárraga-Espinosa et al. 2023). Conversely, extending the model to the relative clustering of saddle-void events (which wall disappears when?) is also of interest, as the latter may impact spin flip, and is dual to void mergers, and as such could act as a cosmic probe for dark energy. One could compute the conditional merger rate subject to a larger-scale saddle-point as a proxy to the influence of the larger-scale cosmic web, following both Musso et al. (2018) and Cadiou et al. (2020) to shed light on how the cosmic web drives galaxy assembly (Kraljic et al. 2018; Laigle et al. 2018; Hasan et al. 2023). Eventually, such a theory could contribute to predicting the expected rate of starburst or AGN activity as a function of redshift and location in the cosmic web.

## ACKNOWLEDGEMENTS

We thank S. Codis and S. Prunet for their early contributions to this work via the co-supervision of EPP's master, and J. Devriendt and M. Musso for insightful conversations. We also thank the KITP, which was supported in part by grant NSF PHY-1748958, for hosting the workshop 'CosmicWeb23: connecting Galaxies to Cosmology at High and Low Redshift' during which this project was advanced. This work was partially supported by the grant Segal ANR-19-CE31-0017 of the French Agence Nationale de la Recherche and by the National Science Foundation under Grant No. NSF PHY-1748958. CC acknowledges support from the Knut and Alice Wallenberg Foundation and the Swedish Research Council (grant 2019-04659). We thank Stéphane Rouberol for running smoothly the infinity Cluster, where the simulations were performed.

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

Aubert D., Pichon C., Colombi S., 2004, MNRAS, 352, 376

Baldauf T., Codis S., Desjacques V., Pichon C., 2021, Phys. Rev. D, 103, 083530

Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, ApJ, 304, 15

Benson A., Behrens C., Lu Y., 2020, MNRAS, 496, 3371

Bertschinger E., 1998, ARA&A, 36, 599

Bond J. R., Myers S. T., 1996, ApJS, 103, 1

Bond J. R., Cole S., Efstathiou G., Kaiser N., 1991, ApJ, 379, 440

Bond J. R., Kofman L., Pogosyan D., 1996, Nature, 380, 603

Bower R. G., 1991, MNRAS, 248, 332

Bullock J. S., Dekel A., Kolatt T. S., Kravtsov A. V., Klypin A. A., Porciani C., Primack J. R., 2001, ApJ, 555, 240

Cadiou C., Pichon C., Codis S., Musso M., Pogosyan D., Dubois Y., Cardoso J. F., Prunet S., 2020, MNRAS, 496, 4787

Cadiou C., Pontzen A., Peiris H. V., 2021a, MNRAS, 502, 5480

Cadiou C., Pontzen A., Peiris H. V., Lucie-Smith L., 2021b, MNRAS, 508, 1189

Cadiou C., Pontzen A., Peiris H. V., 2022, MNRAS, 517, 3459

Danovich M., Dekel A., Hahn O., Ceverino D., Primack J., 2015, MNRAS, 449, 2087

Fakhouri O., Ma C.-P., Boylan-Kolchin M., 2010, MNRAS, 406, 2267

Galárraga-Espinosa D. et al., 2023, A&A, 684, A63

Geller M. J., Huchra J. P., 1989, Science, 246, 897

Genel S., Genzel R., Bouché N., Naab T., Sternberg A., 2009, ApJ, 701, 2002

Hanami H., 2001, MNRAS, 327, 721

Hasan F. et al., 2023, ApJ, 950, 114

Jiang L., Cole S., Sawala T., Frenk C. S., 2015, MNRAS, 448, 1674

Kaiser N., 1984, ApJ, 284, L9

Kraljic K. et al., 2018, MNRAS, 474, 547

Lacey C., Cole S., 1993, MNRAS, 262, 627

Laigle C. et al., 2018, MNRAS, 474, 5437

Manrique A., Salvador-Sole E., 1995, ApJ, 453, 6

Manrique A., Salvador-Sole E., 1996, ApJ, 467, 504

Monaco P., Theuns T., Taffoni G., 2002, MNRAS, 331, 587

Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121

Musso M., Sheth R. K., 2021, MNRAS, 508, 3634

Musso M., Cadiou C., Pichon C., Codis S., Kraljic K., Dubois Y., 2018, MNRAS, 476, 4877

Neistein E., Dekel A., 2008, MNRAS, 388, 1792

Peacock J. A., Heavens A. F., 1990, MNRAS, 243, 133

Pogosyan D., Gay C., Pichon C., 2009, Phys. Rev. D, 80, 081301

Press W. H., Schechter P., 1974, ApJ, 187, 425

Rey M. P., Pontzen A., 2018, MNRAS, 474, 45

Roth N., Pontzen A., Peiris H. V., 2016, MNRAS, 455, 974

Ryden B. S., 1988, ApJ, 329, 589

Salvador-Solé E., Manrique A., Botella I., 2022, MNRAS, 509, 5305

Sohn J., Geller M. J., Hwang H. S., Fabricant D. G., Utsumi Y., Damjanov I., 2023, ApJ, 945, 94

Sousbie T., 2011, MNRAS, 414, 350

Sousbie T., Colombi S., Pichon C., 2009, MNRAS, 393, 457

Stein G., Alvarez M. A., Bond J. R., 2019, MNRAS, 483, 2236

Stopyra S., Pontzen A., Peiris H., Roth N., Rey M. P., 2021, ApJS, 252, 28

Vitvitska M., Klypin A. A., Kravtsov A. V., Wechsler R. H., Primack J. R., Bullock J. S., 2002, ApJ, 581, 799

Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nat. Rev. Phys., 2, 42

Wetzel A. R., 2011, MNRAS, 412, 49

White S. D. M., 1984, ApJ, 286, 38
## APPENDIX A: NOTATIONS

Let us first introduce the dimensionless quantities for the density field, smoother over a scale R by a filter function W

$$ \delta(\mathbf{r}, R) = \int \frac{d^3 k}{(2\pi)^3} \delta(\mathbf{k}) W(kR) e^{ik \cdot \mathbf{r}}. \qquad (A1) $$

We will consider the statistics of this field and its derivatives in this paper. For practical purposes, let us introduce

$$ \sigma_i^2(R) = \frac{1}{2\pi^2} \int_0^\infty dk k^2 P(k) k^{2i} W^2(kR). \qquad (A2) $$

These are the variance of the density field $\langle \delta^2 \rangle = \sigma_0^2$, of its first derivative $\langle \nabla \delta \cdot \nabla \delta \rangle = \sigma_i^2$, etc.

Following Pogosyan, Gay & Pichon (2009), let us introduce the characteristic scales of the field

$$ R_0 = \frac{\sigma_0}{\sigma_1}, \quad R_* = \frac{\sigma_1}{\sigma_2}, \quad \tilde{R} = \frac{\sigma_2}{\sigma_3}. \qquad (A3) $$

These scales are ordered as $R_0 \ge R_* \ge \tilde{R}$. These are the typical separation between zero-crossing of the field, mean distance between extrema, and mean distance between inflection points (Bardeen et al. 1986; Cadiou et al. 2020). Let us further define a set of spectral parameters that depend on the shape of the underlying power spectrum. Out of the three scales introduced above, two dimensionless ratios may be constructed that are intrinsic parameters of the theory

$$ \gamma = \frac{R_*}{R_0} = \frac{\sigma_1^2}{\sigma_0 \sigma_2}, \quad \tilde{\gamma} = \frac{\tilde{R}}{R_*} = \frac{\sigma_2^2}{\sigma_1 \sigma_3}. \qquad (A4) $$

From a geometrical point of view, $\gamma$ specifies how frequently one encounters a maximum between two zero-crossings of the field, while $\tilde{\gamma}$ describes, on average, how many inflection points there are between two extrema. These scales and scale ratios fully specify the correlations between the field and its derivative (up to third order) at the same point. For power-law power spectra of index n, $P(k) \propto k^n$, with Gaussian smoothing at the scale R in 3D, we have $R_0 = R\sqrt{2/(n+3)}$, $R_* = R\sqrt{2/(n+5)}$ and $\tilde{R} = R\sqrt{2/(n+7)}$ while $\gamma = \sqrt{(n+3)/(n+5)}$ and $\tilde{\gamma} = \sqrt{(n+5)/(n+7)}$.

## APPENDIX B: PEAK-EVENT CORRELATION

In order to compute the number of mergers in the vicinity of a peak, we need to compute the two-point correlation function between critical events and peaks. We achieve this by evaluating the joint density of peaks and critical events

$$ n_{\text{peak,ce}} = \langle \text{Peak}(\mathbf{x}) \text{Event}(\mathbf{y}) \rangle, \qquad (B1) $$

where the average is taken over all 30 random variables defining the peak and critical event field up to the third derivative.

$$ \text{Peak}(\mathbf{x}) = |J_x| \Theta_H(-\text{Tr}(\mathbf{H}_x)) \Theta_H(\sum_i \mathbf{H}_{x,i}) \Theta_H(\det \mathbf{H}_x) \times \delta_D(x_1) \delta_D(x_2) \delta_D(x_3) \delta_D(x - \nu_{pk}), \qquad (B2) $$

$$ \text{Event}(\mathbf{y}) = |J_y| \Theta_H(-\text{Tr}(\mathbf{H}_y)) \Theta_H(\sum_i \mathbf{H}_{y,i}) \delta_D(\det \mathbf{H}_y) \times \delta_D(y_1) \delta_D(y_2) \delta_D(y_3) \delta_D(y - \nu). \qquad (B3) $$

Here, $\mathbf{x} = \{x, x_1, \dots, x_3, \mathbf{H}_x, \mathbf{H}_{x,1}, \dots, \mathbf{H}_{x,3}\}$ with $x, x_i, \mathbf{H}_x, \mathbf{H}_{x,i}$ the density, its gradient, its Hessian, and the minors of its Hessian at the peak location and $\mathbf{y} = \{y, y_1, \dots, y_3, \mathbf{H}_y, \mathbf{H}_{y,1}, \dots, \mathbf{H}_{y,3}\}$ with $y, y_i, \mathbf{H}_y, \mathbf{H}_{y,i}$ the density, its gradient, its Hessian and the minors of the Hessian at the critical event location. $|J_x| = |\det \mathbf{H}_x|$ and we provide an explicit formula for $|J_y|$ in equation (C2). This expression is the rotationally invariant equivalent of equation (18) of Cadiou et al. (2020).

The two-point correlation function can then be found as

$$ 1 + \xi_{\text{ce,pk}}(\nu_{\text{pk}}, \nu, R, r) = \frac{\langle \text{Peak}(\mathbf{x}) \text{Event}(\mathbf{y}) \rangle}{\langle \text{Peak}(\mathbf{x}) \rangle \langle \text{Event}(\mathbf{y}) \rangle}. \qquad (B4) $$

Note that, compared to Cadiou et al. (2020), we cannot perform the integration in the frame of the Hessian here. Indeed, the numerator involves cross-correlation between the peak and the critical event which breaks the rotational invariance assumption. While the exact integration cannot be carried out analytically, we can none the less compute it numerically.
APPENDIX C: COVARIANCE MATRICES

The covariance matrix of $x, x_i, x_{ij}$ is given by

$$ C = \begin{pmatrix} 1 & 0 & 0 & 0 & \frac{\gamma_2}{3} & \frac{\gamma_2}{3} & \frac{\gamma_2}{3} & 0 & 0 & 0 \\ 0 & \frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \frac{\gamma_2}{3} & 0 & 0 & 0 & \frac{1}{5} & \frac{1}{15} & \frac{1}{15} & 0 & 0 & 0 & 0 \\ \frac{\gamma_2}{3} & 0 & 0 & 0 & \frac{1}{15} & \frac{1}{15} & \frac{1}{15} & 0 & 0 & 0 & 0 \\ \frac{\gamma_2}{3} & 0 & 0 & 0 & \frac{1}{15} & \frac{1}{15} & \frac{1}{5} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{15} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{15} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{15} & 0 \end{pmatrix}. \quad (C1) $$

We provide in equation (C2) the expression for the jacobian required to compute the number density of critical events in a covariant form. Here $x_{ijk}$ stands for the derivative of the field of order $i$ w.r.t. the first direction, $j$ w.r.t. the second direction, and $k$ w.r.t. the third direction divided by the corresponding RMS defined by equation (A2).

$$ \frac{J_y}{R} = \det \mathbf{H}_y \times (\nabla^2 \nabla y) \cdot \mathbf{H}_y^{-1} \cdot \nabla \det \mathbf{H}_y = \begin{aligned} & \left[ y_{002} y_{020} y_{300} - 2y_{002} y_{110} y_{210} + y_{002} y_{120} y_{200} + 2y_{011} y_{101} y_{210} + 2y_{011} y_{110} y_{201} - 2y_{011} y_{111} y_{200} - y_{011}^2 y_{300} \right. \\ & \left. \qquad - 2y_{020} y_{101} y_{201} + y_{020} y_{102} y_{200} + 2y_{101} y_{110} y_{111} - y_{101}^2 y_{120} - y_{102} y_{110}^2 \right] \\ & \times \left( (y_{012} + y_{030} + y_{210})(y_{011} y_{101} - y_{002} y_{110}) + (y_{002} y_{020} - y_{011}^2)(y_{102} + y_{120} + y_{300}) + (y_{003} + y_{021} + y_{201})(y_{011} y_{110} - y_{020} y_{101}) \right) + \\ & \left[ y_{002} y_{020} y_{210} + y_{002} y_{030} y_{200} - 2y_{002} y_{110} y_{120} - 2y_{011} y_{021} y_{200} + 2y_{011} y_{101} y_{120} + 2y_{011} y_{110} y_{111} - y_{011}^2 y_{210} \right. \\ & \left. \qquad + y_{012} y_{020} y_{200} - y_{012} y_{110}^2 - 2y_{020} y_{101} y_{111} + 2y_{021} y_{101} y_{110} - y_{030} y_{101}^2 \right] \\ & \times \left( (y_{003} + y_{021} + y_{201})(y_{020} y_{200} - y_{110}^2) + (y_{012} + y_{030} + y_{210})(y_{101} y_{110} - y_{011} y_{200}) + (y_{102} + y_{120} + y_{300})(y_{011} y_{110} - y_{020} y_{101}) \right) + \\ & \left[ y_{002} y_{020} y_{201} + y_{002} y_{021} y_{200} - 2y_{002} y_{110} y_{111} + y_{003} y_{020} y_{200} - y_{003} y_{110}^2 - 2y_{011} y_{012} y_{200} + 2y_{011} y_{101} y_{111} \right. \\ & \left. \qquad + 2y_{011} y_{102} y_{110} - y_{011}^2 y_{201} + 2y_{012} y_{101} y_{110} - 2y_{020} y_{101} y_{102} - y_{021} y_{101}^2 \right] \\ & \times \left( (y_{102} + y_{120} + y_{300})(y_{011} y_{101} - y_{002} y_{110}) + (y_{002} y_{200} - y_{101}^2)(y_{012} + y_{030} + y_{210}) + (y_{003} + y_{021} + y_{201})(y_{101} y_{110} - y_{011} y_{200}) \right). \end{aligned} \quad (C2) $$

C1 Numerical implementation

In this section, we describe how Fig. 6 was obtained. We used Monte-Carlo integration to numerically evaluate equation (11). The statistical distribution of the aforementioned field variables is regulated by its 30 × 30 covariance matrix $\Sigma$, which we may compute symbolically for given separation distance $r$, smoothing scales $R_{pk}$ $R_{ce}$, and spectral index $n_s$. We aim to sample points following this distribution and evaluate the integrands of equation (B2) to obtain the expectancies $\langle \text{Peak}(\mathbf{x}) \text{Event}(\mathbf{y}) \rangle$. Ideally, we would like to evaluate the correlation function at identical smoothing scale. However, at equal smoothing scale and small separation, the covariance matrix becomes almost singular, resulting in unfavourable numerical artifacts plaguing the statistic of our result. To avoid this, we instead consider slightly different smoothing scales.

Of course, we can not naively do this directly, as $x, x_i, y, y_i$, $\det \mathbf{H}_y$ will rarely take the specific values that the Dirac deltas impose, and thus the integrands will almost always evaluate to zero. To circumvent this, we compute the distribution of the other field variables conditional to the Dirac deltas being satisfied, with two exceptions: as we will need to integrate $\nu_{pk}$ over a range to recover equation (11) anyway, we replace $\delta_D(y - \nu_{pk})$ by $\Theta_H(\nu_{peak}^\text{max} - y)\Theta_H(y - \nu_{peak}^\text{min})$; furthermore, as $\det \mathbf{H}_y$ depends on several field variables, computing the conditional distribution is not easily feasible. We therefore replace $\delta_D(\det \mathbf{H}_y)$ by a 'thick Dirac delta', $\delta_D^\text{th}(\det \mathbf{H}_y) = \frac{1}{2\epsilon}\Theta_H(\epsilon - \det \mathbf{H}_y)\Theta_H(\det \mathbf{H}_y + \epsilon)$, for some small $\epsilon$. The smaller $\epsilon$ is, the thinner the Dirac delta and the more faithful it is to the original integral. However, reducing $\epsilon$ comes at the cost of reducing the convergence rate of the Monte-Carlo approach.

Once the conditional distribution is computed, we simply draw random points following this distribution and average the evaluation of the integrand on these points, which converges to the value of the integral. This allows to compute $\int_{\nu_{pk}}^{\nu_{pt}} d\nu \langle \text{Peak}(\mathbf{x}) \text{Event}(\mathbf{y}) \rangle$ and $\langle \text{Peak}(\mathbf{x}) \rangle$, and results in Fig. 6. Practically, to obtain Fig. 6, we drew $2^{26}$ sample points. We used a smoothing ratio $R_{pk}/R_{ce} = 0.95$ and a thick Dirac delta of size $\epsilon = 10^{-3}$.
## APPENDIX D: EVENT DIFFERENTIAL DISTRIBUTION

The differential density of critical events with different field values $v$ are derived in Cadiou et al. (2020) and read

$$ n_{ce}(v, R) = \frac{3R}{\tilde{R}^2 R_*^3} \frac{1}{5} \left(\frac{3}{2\pi}\right)^{3/2} (1 - \tilde{\nu}^2) \sum_{i=5,6,9} c_{3,i}(v, \gamma) \exp\left(-\frac{\nu^2}{2(1 - 5\gamma^2/i)}\right), \quad (D1) $$

where the full expressions for $c_{3,i}(v, \gamma)$ were given in equation (39) of Cadiou et al. (2020) Factorizing the mean density of the critical events, $\tilde{n}_{ce}$, in this expression such that

$$ n_{ce}(v, R) = \tilde{n}_{ce}C(v), \quad \text{with} \quad \int_{-\infty}^{\infty} dv C(v) = 1, \quad (D2) $$

we introduce the normalized distribution of the events in height

$$ C(v, \gamma) = \frac{18\sqrt{10\pi}}{29 - 6\sqrt{6}} \sum_{i=5,6,9} c_{3,i}(v, \gamma) \exp\left(-\frac{\nu^2}{2(1 - 5\gamma^2/i)}\right), $$

which depends on the power spectrum (and on the smoothing scale if the spectrum is not scale invariant) only through the spectral parameter $\gamma$. In Fig. D1, we show the behaviour of $C(v, \gamma)$ for several choices of $\gamma$ that correspond to the spectral slopes of interest $n = -2.5, -2, -1.5$, see equation (A4). For comparison, the peak rarity distribution is represented as dashed lines. We see that critical events typically occur at lower $v$ than peaks, and are very rare for $v > 3$.

We note that $C(v, \gamma)$ for $v < 3$ can be approximated with an accuracy of better than $10^{-3}$ for $-2.5 \le n_s \le -1.5$ by a Gaussian with parameters

$$ \tilde{C}(v, \gamma) = \frac{1}{\sqrt{2\pi}\tilde{\sigma}} \exp\left(-\frac{(v - \tilde{\mu})^2}{2\tilde{\sigma}^2}\right), \\ \tilde{\mu} = 1.546\gamma, \quad \tilde{\sigma} = 1 - 0.5084\gamma^2 + 0.04140\gamma. $$

We show on the top panel of Fig. D1 the residuals of this fit.

<figure><img src="image_12.png" /><figcaption>**Figure D1.** Normalized differential distribution of critical events in density rarity value for several values of spectral slope, as marked. The top panel shows the absolute residuals when fitting $C(v, \gamma)$ with a Gaussian function. Parameters of the Gaussian are provided in the text. We show the peak rarity distribution as thin dashed lines.</figcaption></figure>

## APPENDIX E: MEAN EVOLUTION OF PEAK RARITY

Let us predict the evolution of peak rarity with Gaussian smoothing. The value of the density field smoothed with a Gaussian window at position $x$ changes with window size $R$ according to the diffusion equation

$$ \frac{\partial \delta(\mathbf{r}, R)}{\partial R} = R \Delta \delta(\mathbf{r}_{pk}, R). \quad (E1) $$
If we are interested in the change of the peak value along the peak track $r_{pk}(R)$, we notice that the partial derivative $\partial/\partial R$ can be replaced by the full derivative, since at the peak $\nabla\delta = 0$, thus

$$ \frac{d\delta(r_{pk}, R)}{dR} = R\Delta\delta(r_{pk}, R), \qquad (E2) $$

and for brevity, we can now drop the reference to the position of the peak in the argument. In terms of the rarity of the peak, $\nu(R) = \delta(R)/\sigma_0(R)$, the evolution equation becomes

$$ \frac{d\nu(R)}{d\log R} = -\nu(R) \frac{d\log \sigma_0}{d\log R} + R^2 \Delta\nu(R). \qquad (E3) $$

The process of changing the peak height with smoothing is stochastic since the Laplacian of the field is a random quantity. However, to estimate the mean change, let us approximate the stochastic Laplacian by its conditional mean value given the height of the peak

$$ \Delta\nu(R) \rightarrow \langle\Delta\nu(R)|\nu(R), \lambda_i < 0\rangle, \qquad (E4) $$

where the last conditional inequality, written compactly in terms of the eigenvalues $\{\lambda_i\}_{i=1,\dots,D}$ of the Hessian of the density, ensures that we are dealing with a maximum and not an arbitrary field point. In principle, one should also enforce the vanishing gradient of the field at peak position, but since the gradient of the field is uncorrelated with both the value and the second derivatives of the field at the same point, this condition is inconsequential to our problem.

Let us first evaluate the conditional mean in equation (E4) in 1D, where the peak condition is just $\Delta\nu < 0$. This gives us

$$ \langle\Delta\nu|\nu, \Delta\nu < 0\rangle = -\frac{\nu}{R_0^2} - \frac{1}{\zeta R_0^2} \sqrt{\frac{2}{\pi}} \frac{\exp\left(-\frac{\nu^2\zeta^2}{2}\right)}{1 + \operatorname{erf}\left(\frac{\nu\zeta}{\sqrt{2}}\right)}, \qquad (E5) $$

where $R_0$ is defined by equation (A3), while $\zeta = \gamma/\sqrt{1-\gamma^2}$ varies from zero at $\gamma = 0$, to infinity at $\gamma = 1$. The first term in equation (E5) is the general conditional response of the Laplacian to the field value, while the second correction comes from restricting the field to be at the local maximum. For a Gaussian window function, we also have that $d\log \sigma_0/dR + R/R_0^2 = 0$ for any power spectrum. We can then obtain the (mean) evolution equation for the peak rarity

$$ \frac{d\nu(R)}{d\log R} = -\frac{R^2}{\zeta R_0^2} \sqrt{\frac{2}{\pi}} \frac{\exp\left(-\frac{\nu^2\zeta^2}{2}\right)}{1 + \operatorname{erf}\left(\frac{\nu\zeta}{\sqrt{2}}\right)}, \qquad (E6) $$

where $R/R_0$ is a constant that depends on the power spectrum (cf Section A); it is equal to $\zeta$ for power-law spectra. Thus, the rarity of maxima is, on average, decreasing slowly with smoothing scale $R$. The decrease is less pronounced for rarer peaks and can be in first approximation neglected for $\nu\zeta > \sqrt{2}$.

For power-law spectra in $D$ dimensions, $\zeta = \sqrt{(n+D)/2}$. We obtain a criterion for insignificant rarity drift for peaks rarer than $\nu > 2/\sqrt{n+D}$. Note that this is based on generalizing the 1D result to higher dimensions. In 2D the drift term is more cumbersome and in 3D it does not have a closed analytical form, but the structure of the result remains the same.

This paper has been typeset from a TeX/LATEX file prepared by the author.