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

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 $\sim 400$ power-law Gaussian random fields smoothed at $\sim 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 percent.


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 et al. 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 et al. 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 et al. 2002;Salvador-Solé et al. 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 ⋆ E-mail: corentin.cadiou@fysik.lu.se 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 et al. 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 PDF of the spin parameter that each merger contributes as a function of rarity.Section 4 wraps up.
Appendix A introduces the relevant spectral parameters.Appendix B revisits event statistics while relying on the clustering properties of peaks and events computed from first principles.Appendix D recalls the critical event PDF and provides a fit to it.Appendix E discusses the evolution of peak rarity with Gaussian smoothing.
1 2 3 Increasing smoothing scale Figure 1.Typical density field around a critical event in 1D.From left to right: 1.We see two distinct objects (in red), separated by a minimum (in blue) in the density field.2. The critical event occurs (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 (in red).
Position Smoothing scale

Density
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.

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 δ(R) of a given peak smoothed at scale R can be mapped to a formation redshift z, using the spherical collapse model, δ(RTH)D(z) = δc, where δc = 1.686 is the critical density, D(z) is the growing mode of perturbations and RTH is the top-hat scale.
The mass of the halo formed at this time is inferred from the top-hat window scale, M = ρm 4 3 πR 3 TH .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 Figure 3. Visualisation 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 peaksaddle 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 contains all the past physical history of a given peak.This paper aims to characterise 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.
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 halos 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 R0 and position x pk (R0).Since we employ Gaussian filters, we need to match the Gaussian and Top-Hat smoothing scales, R and RTH to assign mass to haloes.The criteria of equal mass encompassed by the filter1 gives RTH ≈ 1.56R in 3D.At R0 the peak describes a halo that has collected its mass from a spherical cross-section of (N + 1)D space of volume ∝ R N 0 ; we can call this sphere a Lagrangian patch of

Cone of accretion
Critical event

Position
Directly connected critical events 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 6 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 αR.Here, α = 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.
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 ∝ 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 ∼ R N around it as shown in Fig. 3 for (2+1)D example.
A critical event marks the end of a trajectory of a peak that disappeared when its scale reached Rce 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 R0 around the final surviving peak position x pk (R0) 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 = Rce, it is connected to the main peak in N D 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 αR, where α 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 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 peakpeak 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 Section 2.3.2 and Appendix 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 ∝ 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.

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 M0 as its mass grew from f 3 M0 to M0, where f < 1.Since smoothing scale maps directly onto mass, this amounts to counting the number of critical events between two scales R0 ∝ 3 √ M0 and R1 ∝ f 3 √ M0.Let us now count direct major mergers that we define as mergers with satellites that bring at least f 3 M0 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 M0 to M0, is given by the number of critical events in the section of the cone of influence of the halo contained between f R0 and R0.It can be obtained with the following integral where nce(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 α, introduced in the previous section, sets the radial extent of the R = 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 Nmerger by taking the den-sity of the critical events inside the cone of influence to be equal to its 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 where the spectral scale R * and parameter γ are given in Appendix A, together with other parameters that characterize the statistics of a density field.Note that R * ∝ R, so the number density of critical events scales as nce ∝ 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 ns and Gaussian filter, equation (2) gives nce ≈ 0.0094 5 + ns 2 Using this mean value nce for critical density in equation ( 1) gives us a first rough, but telling estimate of the merger number that a typical halo experiences The final results depend on the specific values of f and α.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 α 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 ns = −2, which corresponds to the typical slope of a Λ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/ 3 √ 10).For the choice of α, a natural idea would be to count only critical events within the halo Top-Hat Lagrangian radius, then α = RTH/R ≈ 1.56 and we find Nmerger ≈ 0.2.However, this is clearly an underestimation for α, as it corresponds to the situation where the center 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 α = 2RTH/R ≈ 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 Nmerger ≈ 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/ 3 √ 1000), did it having experienced Nmerger ≈ 5 major mergers.These conclusions follow directly from first principles while studying the structure of the initial density field.
shows the tracks of two peaks (red) and a saddle point (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 Rce, 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 δ C and thus can be viewed as reaching the threshold δ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 = (Rce/R pk ) 3 at a redshift corresponding to D(z) = δc/δ A .

Accounting for rarity and clustering
Let us now refine our model so as to define α more rigorously by i) requiring the merging object to have gravitationally collapsed before 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 ν as a function of their distance r to a central peak, both defined at the same smoothing scale R where C is the distribution function, ∞ −∞ C(ν, γ) dν = 1, of critical events in overdensity and is given by the analytical formula in Appendix 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, ξ ce,pk (ν, 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 ν and what extent of r one needs to consider to count haloes merging into a surviving halo with a particular ν pk .

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 δ pk (R pk ) = δc/D(z), where δ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 Rce describes the merger of a satellite halo of scale Rce into the main peak at scale R pk that corresponds to the same redshift, i.e such that δce(Rce) = δ pk (R pk ) = δ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 R0 and that experienced a merger in its past when its scale was R pk , R pk ⩽ R0.Requiring that, at the time of the merger, the surviving peak's scale (mass) is larger than that of the satellite sets the relation R0 ⩾ R pk ⩾ Rce or, conversely, σ0(R0) ⩽ σ0(R pk ) ⩽ σ0(Rce), where σ0 is the r.m.s. of the density field and is defined in Appendix A. From which it follows that the rarity of the relevant critical events at scale Rce, is within the range The lower bound corresponds to mergers that have been completed at the very last moment, R pk = R0.The upper bound is achieved for mergers of two equal mass objects, R pk = Rce.
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.

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 (ν, R, r) → 0 far from the peak, i.e. when r → ∞.While we cannot implement this condition analytically (as it is non-local), we can measure n ce|pk (ν, 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 n ce|pk given by using formally straightforward analytical calculations of critical event -peak correlations, as described in Appendix 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 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: Therefore, in this approximation, the question of where to truncate the integration over the peak neighbourhood remains and we have For scale-free spectra, introducing the dimensionless ratios u = r/R and w = R/R0 and changing the order of integration, equa-pk 2.Here νce is integrated from 0.681ν pk to ν pk , with power spectrum n = −2 and smoothing ratio R pk /Rce = 0.95.We detail the method used to compute this graph in Appendix C1.Rarer peaks have a stronger excess of events in their vicinity, which mirrors the known Kaiser bias (Kaiser 1984).
tion ( 8) can be written as 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 ν in the range 0.681ν pk ⩽ ν ⩽ ν pk which corresponds to (R/R0) 3 = 1/10 for n = −2 power-law spectrum.Fig. 6 shows that critical events are indeed preferentially clustered at distances ∼ 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 α ≈ 3 in the estimate of equation (4), i.e. qualitatively double the value of RTH.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.

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 et al. 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.

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 to a given peak.We proceed as follows: we generate 430 scale-invariant gaussian random fields of size 256 3 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 % increase in the mass associated with the smoothing scale, and 8 steps give the mass increasing approximately 10-fold.The skeleton is identified at each scale using Dis-Perse (Sousbie et al. 2009) with a persistence threshold of 1 % 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, Figure 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 rarity2 .Critical events are counted in the interval of rarity 0.681 ν pk ⩽ νce ⩽ ν 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.Rce = 1/ 3 √ 10R0 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 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 2 .Instead, in Fig. 7 we find an exclusion zone at small separations, an excess probability at r/R ∼ 2−3, and a cut-off 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 % lower values in Table 1 relative to equation (4) (per ln f ) 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 1.5 pk < 2.0 2.0 pk < 2.5 2.5 pk < 3.0 3.0 pk < 3.5 3.5 pk < Figure 7.The radial number density of topologically linked critical event 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.5R, 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.
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..5, 2.0[ [2.0, 2.5[ [2.5, 3.0[ [3.0, 3.5[ [3.5 As a function of peak rarity ν pk , the mean number of critical events per peak in a constant R slab first increases up to ν pk ≈ 3 before decreasing.It is fairly flat in the range ν pk = 2.5 − 3.5, where most of the physical haloes are, which argues for taking α parameter independent on ν 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 multi-scale 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 ⩾ Mce ⩾ M pk /10 that have δce(Rce) ⩾ δ pk (R pk ) where Rce and R pk are now different.We only retain peaks that exist at scales R pk ⩾ Rmin 3 √ 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/10M 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 1.5 pk < 2.0 2.0 pk < 2.5 2.5 pk < 3.0 3.0 pk < 3.5 3.5 pk < Figure 8. Distribution of the number of major mergers for peaks in different rarity bins, as labeled.The mean number of major mergers increases from 1.8 for low-ν peaks to 2 for high-ν peaks and is graphically represented with vertical lines with tilted caps for readability.
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), α = 2RTH/R = 3.1 is a very good proposition.

Orbital spin parameter of mergers
Let us now define the spin parameter of an event at scale Rce relative to a peak of rarity ν at scale R as Here, we recall that σ−1 is the variance of the gradient of the potential, whose definition is given in Appendix A. ∇ψ and ∇ψ pk are the gravitational accelerations at the locus of the critical event and of the peak respectively, and 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 3 ce |r × (∇ψ − ∇ψ 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 ∝ R 3 pk coming from a distance of R pk with velocity equal to the r.m.s. of the velocity √ 2σ−1 would have at the time of the collapse of the peak.We show on Fig. 9 the distribution of the spin parameter λ.The sample consists of mergers with mass ratios greater or equal to 1:10 around rare peaks.In practice, we select ν pk > 2.5; R pk / 3 √ 10 ⩽ Rce ⩽ R pk and δce ⩾ δ pk .The distribution is roughly log-normal with a mean value of µ λ ≈ 0.048 and a standard deviation of σ λ ≈ 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 et al. 2004;Danovich et al. 2015).This translates two effects.First, the mass ratio (R 3 ce /R 3 pk ) is dominated by the more numerous smaller mergers, which drives λ down.In addition, we find that critical events tend to have a mostly radial acceleration so that their tangential velocity is smaller than the r.m.s. of the velocity field ((∇ψ − ∇ψ pk )/σ−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 et al. 2020).

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 −3 pk , the sample size quickly decreases with smoothing scale; we select here peaks that exist at a scale larger than Rmin 3 √ 100 ≈ 4.6Rmin.The practical consequence is that our sample is only complete over two decades in mass ratio.We also only retain rare peaks that have ν 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 tce respectively, using δ = δ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 ∆ = (t pk − tce)/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 dataset 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 (ν pk ⪆ 3.5) typically Figure 10.Top: distribution of the merger as a function of mass ratio and lookback time.We show the median with 68 % 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.
had 80 % of their mergers in the second half of their life (past 7 Gyr for a 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 (Genel et al. 2009;Fakhouri et al. 2010, figure 9), 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 σ(R)−δ(R) are mapped onto 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.

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); 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 focussed here on the analysis of merger events that bring at least 10 % 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 zero th -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 multi-scale 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 2 major mergers.By tracking peaks in positionsmoothing 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 log-normal distribution with a mean of µ λ = 0.048 and σ λ = 0.5.These properties are remarkably close to the distribution of DM halo spins measured in hydrodynamical simulations (µ λ = 0.038, σ λ = 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 et al. 2021a(Cadiou et al. , 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 multi-scale analysis of these initial conditions.

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 (dis-tance 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 multi-scale 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 et al. 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 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 saddlepoint 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.
1 + ξ ce,pk (ν pk , ν, R, r) = Peak(x) Event(y) ⟨Peak(x)⟩⟨Event(y)⟩ . (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 nonetheless compute it numerically.

APPENDIX C: COVARIANCE MATRICES
The covariance matrix of x, xi, xij is given by 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).

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 Σ, which we may compute symbolically for given separation distance r, smoothing scales R pk and Rce and spectral index ns.We aim to sample points following this distribution and evaluate the integrands of equation ( B2) to obtain the expectancies Peak(x) Event(y) .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 unfavorable 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, xi, y, yi, det Hy 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 ν pk over a range to recover equation ( 11) anyway, we replace δD(y − ν pk ) by ΘH(ν max peak − y)ΘH(y − ν min peak ); furthermore, as det Hy depends on several field variables, computing the conditional distribution is not easily feasible.We therefore replace δD(det Hy) by a 'thick Dirac delta', δ ε D (det Hy) = 1 2ϵ ΘH(ε − det Hy)ΘH(det Hy + ε), for some small ε.The smaller ε is, the thinner the Dirac delta and the more faithful it is to the original integral.However, reducing ε 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 νp νpw dν⟨Peak(x) Event(y)⟩ and ⟨Peak(x)⟩, and results in Fig. 6.Practically, to obtain Fig. 6, we drew 2 26 sample points.We used a smoothing ratio R pk /Rce = 0.95 and a thick Dirac delta of size ε = 10 −3 .

Figure 6 .
Figure6.Integration of the analytical expression of the clustering of critical events around peaks, equation (11), for different values of ν pk .Here νce is integrated from 0.681ν pk to ν pk , with power spectrum n = −2 and smoothing ratio R pk /Rce = 0.95.We detail the method used to compute this graph in Appendix C1.Rarer peaks have a stronger excess of events in their vicinity, which mirrors the known Kaiser bias(Kaiser 1984).

Figure 9 .
Figure9.Merging objects bring in orbital spin; we quantify here its distribution.The distribution resembles a log-normal distribution with parameters µ λ ≈ 0.048 and σ λ ≈ 0.51; the distribution resembles its N -body counterpart.