<figure><img src="image_1.png" /></figure>

HAL
open science

# On the connectivity of the cosmic web: theory and implications for cosmology and galaxy formation

Sandrine Codis, Dmitri Pogosyan, Christophe Pichon

► To cite this version:

Sandrine Codis, Dmitri Pogosyan, Christophe Pichon. On the connectivity of the cosmic web: theory and implications for cosmology and galaxy formation. Monthly Notices of the Royal Astronomical Society, 2018, 479 (1), pp.973-993. 10.1093/mnras/sty1643. hal-01768012

HAL Id: hal-01768012

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

Submitted on 9 Nov 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# On the connectivity of the cosmic web: theory and implications for cosmology and galaxy formation

Sandrine Codis,<sup>1,2*</sup> Dmitri Pogosyan<sup>3</sup> and Christophe Pichon<sup>1,4,5</sup>

Accepted 2018 June 19. Received 2018 June 19; in original form 2018 March 30

## ABSTRACT

Cosmic connectivity and multiplicity, i.e. the number of filaments globally or locally connected to a given cluster is a natural probe of the growth of structure and in particular of the nature of dark energy. It is also a critical ingredient driving the assembly history of galaxies as it controls mass and angular momentum accretion. The connectivity of the cosmic web is investigated here via the persistent skeleton. This tool identifies topologically the ridges of the cosmic landscape which allows us to investigate how the nodes of the cosmic web are connected together. When applied to Gaussian random fields corresponding to the high-redshift universe, it is found that on average the nodes are connected to exactly $\kappa = 4$ neighbours in two dimensions and $\sim 6.1$ in three dimensions. Investigating spatial dimensions up to $d = 6$, typical departures from a cubic lattice $\kappa = 2d$ are shown to scale like the power 7/4 of the dimension. These numbers strongly depend on the height of the peaks: the higher the peak the larger the connectivity. Predictions from first principles based on peak theory are shown to reproduce well the connectivity and multiplicity of Gaussian random fields and cosmological simulations. As an illustration, connectivity is quantified in galaxy lensing convergence maps and large dark haloes catalogues. As a function of redshift and scale the mean connectivity decreases in a cosmology-dependent way. As a function of halo mass, it scales like 10/3 times the log of the mass. Implications on galactic scales are discussed.

**Key words:** methods: analytical – methods: numerical – galaxies: formation – dark energy – large-scale structure of Universe.

## 1 INTRODUCTION

Over the course of the last decades, our understanding of the extragalactic universe has significantly evolved: the description of its components has evolved from being (essentially) isolated to being multiply connected both on large scales, cluster and galactic scales. This interplay between large and small scales is driven in part by gravity which tends to couple dynamically different scales in the framework of the so-called concordant cosmological model (de Bernardis 2000). This model predicts a certain shape for the initial conditions, leading to a hierarchical formation scenario, which produces the large-scale structure, the most striking feature in the distribution of matter on megaparsecs scale. This distribution was observed more than 30 yr ago by the first CfA catalogue (de Lapparent, Geller & Huchra 1986) followed by many others such

as the SDSS (Adelman-McCarthy 2008), 2dF (Cole 2005), or more recently DES (Abbott et al. 2016) galaxy redshift surveys.

The ‘Cosmic Web’ picture (Klypin & Shandarin 1993; Bond, Kofman & Pogosyan 1996) was developed to explain the origin of this network: it relates the observed clusters of galaxies, and filaments that link them, to the geometrical properties of the initial density field that are enhanced but not yet destroyed by the still mildly non-linear evolution on those scales. It builds from the ellipsoidal collapse model studied by Lynden-Bell (1964) and Lin, Mestel & Shu (1965), followed by Zel'dovich's work (Zel'Dovich 1970) which related the anisotropic nature of the gravitational collapse to the formation of elongated and flattened structures. The concept of cosmic web emerged from those ideas and was extended in the peak-patch formalism (Bond & Myers 1996): the origin of filaments and nodes lies in the asymmetries of the initial Gaussian random field (GRF hereafter) describing the primordial universe which is later amplified by gravitational collapse. These investigations stressed the key role of the theory of random fields in cosmology and the

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

<footnote n="2"><sup>2</sup>Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada</footnote>

<footnote n="3"><sup>3</sup>Department of Physics, University of Alberta, 412 Avadh Bhatia Physics Laboratory, Edmonton, Alberta, T6G 2J1, Canada</footnote>

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

<footnote n="5"><sup>5</sup>Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote>* E-mail: codis@iap.fr</footnote>
importance of non-local tidal effects in weaving the cosmic web. The high-density peaks define the nodes of the evolving cosmic web and completely determine the filamentary pattern in between. Building upon the cosmological peak theory described in Bardeen et al. (1986), local properties of the filamentary cosmic web can be predicted such as the length of filaments, the surface of walls or their curvatures (Pogosyan et al. 2009), while its topology can be fully characterized (Gott, Dickinson & Melott 1986; Mecke, Buchert & Wagner 1994; Matsubara 1994; Gay, Pichon & Pogosyan 2012), all of these observables carrying complementary cosmological information (Zunckel, Gott & Lunnan 2011; Codis et al. 2013).

While traditionally the emphasis is placed on the statistical descriptors of the underlying random field via the hierarchy of N point correlation functions (Scoccamarro et al. 1998; Slepian et al. 2017), in this paper, we focus on the connectivity of this cosmic network as a mean to understand its morphology and geometry. Indeed, recent ridge extractors (such as the skeleton, Novikov, Colombi & Doré 2006) allow for the definition of the filamentary cosmic web as a *connected* network that continuously link maxima and saddle points of a scalar field together. Hence, it is of interest to try and understand the topological and geometrical properties of the underlying density field through the connectivity and hierarchical relationship that the ridges introduce between the critical points. This can be used to establish, in particular, the percolation properties of the Web (Colombi, Pogosyan & Souradeep 2000). Applied to cosmology, these ridges provide a formal definition of the concept of individual filaments. Considering matter distribution on large scales in the Universe, a natural definition of a *single* filament is the subset of the cosmic web that directly links two haloes together. The transposition of such a definition to the skeleton allows the introduction of useful concepts such as neighbouring relationship between haloes in the cosmic web sense. This has implication on both large and small scales.

From the point of view of constraining cosmological parameters, the redshift evolution of the connectivity of the cosmic web on large scales can be used as a probe of fundamental physics. Indeed, it can robustly estimate the growth of structure and therefore the equation of state of dark energy, as the rate of acceleration of the universe disconnects haloes, and gravitational non-linear evolution induces filament coalescence.

From a smaller scale perspective, the importance of the cosmic web's connectivity is sustained by pan chromatic observations of the environment of galaxies which illustrate sometimes spectacular merging processes, following the pioneer work of e.g. Schweizer (1982) (motivated by theoretical investigations such as Toomre & Toomre 1972). The importance of anisotropic accretion on cluster and dark matter (DM) halo scales (Aubert, Pichon & Colombi 2004; Bailin & Steinmetz 2005; Kang & Wang 2015; Aubert & Pichon 2007) down to central galaxies (Kimm et al. 2011; Welker et al. 2018) is now believed to play a significant role in regulating the shape and spectroscopic properties of galaxies. Indeed it has been claimed (see e.g. Ocvirk, Pichon & Teyssier 2008; Dekel et al. 2009) that the geometry of the cosmic inflow on a galaxy is strongly correlated to its history and nature, as can be seen through the observed (Alpaslan et al. 2016; Poudel et al. 2017; Malavasi et al. 2017; Chen et al. 2017; Kraljic et al. 2018; Laigle et al. 2018), virtually measured (among many others Codis et al. 2012; Metuki et al. 2015; Dubois et al. 2014), and predicted (see for instance Kaiser 1984; Codis, Pichon & Pogosyan 2015; Alonso, Eardley & Peacock 2015; Musso et al. 2018) correlations between cosmic web and galactic properties such as its mass, spin, shape, temperature, or entropy distribution. One of the puzzles of galaxy formation involves understanding

<figure><img src="image_4.png" /><figcaption>**Figure 1.** The connectivity, $\kappa$, of maxima vertices on top of the peak patches for a two-dimensional GRF field. The (blue) numbers represents $\kappa$ for the corresponding node. The colour coding of the skeleton (in shades of red) reflects the underlying density of the field.</figcaption></figure>

how galactic discs reform after minor and intermediate mergers, a process which is undoubtedly controlled by the drifting of filaments through cosmic time and anisotropic gas inflow therein which carries coherent angular momentum (Pichon et al. 2011; Prieto et al. 2015). At high redshifts, those streams – the so-called cold flows – can penetrate haloes deeply into their core and feed the central galaxy. The coplanarity of those flows is predicted by numerical simulations (Danovich et al. 2012), evidenced by the observation of planes of satellite galaxies around their central (e.g. Holmberg 1969; Ibata et al. 2013) and could be explained by predicting cosmic connectivity as a function of the rareness of the nodes and prominence of filaments, which is one of the motivations of this paper.

Section 2 first summarizes our understanding of the statistics of extrema, their spatial correlations, and how Morse theory links them to the actual cosmic web. Section 3 then investigates the connectivity of *d* dimensional GRF both locally, and globally within the context of peak theory and numerical simulations of GRFs. Section 4 focuses next on the cosmic connectivity: first for convergence maps then for the three-dimensional DM distribution and dark halo catalogues in concordant Λ cold dark matter (ΛCDM) simulation. Finally, Section 5 wraps up.

## 2 CRITICAL SETS OF GRF

The concept of random fields is central to cosmology. Random fields provide initial conditions for the evolution of the matter distribution in the Universe. They also describe the distribution of the cosmic microwave background on the celestial sphere and are key to understand the geometrical properties of the large-scale structure as was highlighted in this introduction both in the context of cosmology and galaxy formation. An example of a two-dimensional GRF is displayed in Fig. 1. Each peak (dots) is a node of the skeleton (in shades of red) and belongs to $\kappa$ valleys whose boundaries form the filaments. This number $\kappa$ is defined as the connectivity of the node, namely the number of peaks one peak is connected to. In
this section, we will first recap some properties of extrema<fnref n="1" /> in GRF before turning to the theory of the connectivity of GRF in Section 3.

## 2.1 Extrema of GRF

Let us start by reviewing basic facts on the distribution of critical points in GRF. This field will be generically denoted $\rho$ (with a slight abuse of notation since in the cosmological context it will refer to the density contrast instead of the density field itself) and assumed to have zero mean. The one-point distribution of extrema in 2D has been studied in Longuet-Higgins (1957), while an extensive study of 3D extrema in cosmological settings goes back to Bardeen et al. (1986). In 2D, the set of extrema is composed of maxima, saddle points and minima, identified by the signs of the eigenvalues $\lambda_i$, $i \in \{1, 2\}$, of the (Hessian) matrix of the second derivatives of the field denoted $\rho_{ij}$, $i, j \in \{1, 2\}$, as $(-, +, +, +)$, while in 3D, there are two types of saddle points so that extrema come with signatures $(-, -, +, +)$, $(-, +, +, +)$, $(+, +, +, +)$. In this work, we are primarily interested in the distribution of maxima (which we also often call peaks), that are identified with the nodes of the skeleton, and saddle points of ‘filamentary’ type, $(+)$ in 2D or $(+ -)$ in 3D through which the skeleton bridges connecting two maxima necessarily pass as described in Section 2.2 below.

### 2.1.1 Number density of extrema

Let us first recall the well-known results for the total number density of maxima and ‘filamentary’ saddle points in GRF

$$2D: \bar{n}_{\max} = \frac{1}{8\sqrt{3\pi R_*^2}}, \qquad \bar{n}_{\text{sad}} = \frac{1}{4\sqrt{3\pi R_*^2}}, \qquad (1)$$

$$3D: \bar{n}_{\max} = \frac{29\sqrt{15}-18\sqrt{10}}{1800\pi^2 R_*^3}, \qquad \bar{n}_{\text{sad}} = \frac{29\sqrt{15}+18\sqrt{10}}{1800\pi^2 R_*^3}, \qquad (2)$$

where the characteristic length $R_* \equiv \sigma_1/\sigma_2$ is defined by the ratio of the variances of the field first, $\sigma_1^2 = \langle \nabla\rho \cdot \nabla\rho \rangle$, and second, $\sigma_2^2 = \langle (\Delta\rho)^2 \rangle$, derivatives. A characteristic measure of peak separation is given by the radius of a sphere that on average contains exactly one peak, $R_p = (\pi\bar{n}_{\max})^{-1/2} \approx 3.7 R_*$ and $(\frac{4\pi}{3}\bar{n}_{\max})^{-1/3} \approx 4.2 R_*$, in 2D and 3D, respectively. The universal relations

$$2D: \bar{n}_{\text{sad}}/\bar{n}_{\max} = 2, \quad 3D: \bar{n}_{\text{sad}}/\bar{n}_{\max} \approx 3.055, \qquad (3)$$

are very important for connectivity discussion of GRFs. Hereafter, all extrema heights will be measured in units of the rms $\sigma_0 = \sqrt{\langle \rho^2 \rangle}$ of the field by means of the rareness $\nu = \rho/\sigma_0$, all distances in units of $R_p$, and all number densities in units of $\bar{n}_{\max}$. The latter effectively means that instead of number densities we will be quoting the actual number of extrema inside a sphere of radius $R_p$. We shall use $d$ to designate the dimensionality of space wherever general expressions can be used.

We are often interested in the number density of extrema that satisfy some additional properties (e.g. height). Computing it requires the knowledge of the joint probability distribution function of these properties, the gradient, and the second derivatives of the field, all evaluated at zero gradient. The number density of extrema that satisfy these properties is then obtained by averaging

the determinant of the Hessian $\det(\rho_{ij})$ with this distribution over the region of $\rho_{ij}$ with appropriate signs of eigenvalues for the given type of extrema,

$$n_{\text{ext}}(\text{prop}) = \frac{\int d\rho_{ij} P(\text{prop}, \nabla\rho = 0, \rho_{ij}) |\det(\rho_{ij})| \theta(\{\lambda_i\})}{\bar{n}_{\max}}, \qquad (4)$$

the last condition being formally codified in the Heaviside $\Theta$-like function $\theta(\{\lambda_i\})$ equal, in particular to $\theta_p(\{\lambda_i\}) = \prod_{i=1}^D \Theta(-\lambda_i)$ for peaks and $\theta_s(\{\lambda_i\}) = \Theta(\lambda_1) \prod_{i=2}^D \Theta(-\lambda_i)$ for filamentary saddles.

Perhaps the most important property of an extrema is its height $\nu$. We shall denote as $n_{\text{ext}}(\nu)$ the number of extrema (in the spherical volume of radius $R_p$) exceeding the height $\nu$. Fig. 2 shows the correspondent differential numbers of maxima, $\partial_\nu n_{\max}$, and saddle points $\partial_\nu n_{\text{sad}}$ for GRF, which in 2D are given by analytical expressions (e.g. Bardeen et al. 1986; Gay et al. 2012).

$$\partial_\nu n_{\max} \equiv \frac{dn_{\max}}{d\nu} = \frac{\sqrt{3(1-\gamma^2)}}{\pi} \gamma \nu \exp\left(-\frac{\nu^2}{2(1-\gamma^2)}\right) (5) \\ + \sqrt{\frac{3}{\pi(6-4\gamma^2)}} \exp\left(-\frac{3\nu^2}{6-4\gamma^2}\right) \left[1+\text{Erf}\left(\frac{\gamma\nu}{\sqrt{(1-\gamma^2)(6-4\gamma^2)}}\right)\right] \\ + \sqrt{\frac{3}{2\pi}} \gamma^2 (\nu^2 - 1) \exp\left(-\frac{\nu^2}{2}\right) \left[1+\text{Erf}\left(\frac{\gamma\nu}{\sqrt{2(1-\gamma^2)}}\right)\right],$$

$$\partial_\nu n_{\text{sad}} \equiv \frac{dn_{\text{sad}}}{d\nu} = \sqrt{\frac{6}{\pi(3-2\gamma^2)}} \exp\left(-\frac{3\nu^2}{6-4\gamma^2}\right), \qquad (6)$$

while in 3D are evaluated by numerical integrations. The distributions are governed by the spectral parameter $\gamma \equiv \sigma_1^2/(\sigma_0\sigma_2)$ which reflects the correlation between the normalized<fnref n="2" /> field and trace of its Hessian at the same point, $\gamma = -\langle\rho\Delta\rho\rangle/\sigma_0\sigma_2$. We shall consider the extrema as (relatively) high and rare if their height exceeds one standard deviation from the mean. This threshold is $\gamma$-dependent<fnref n="3" /> but as a guidance in future discussions we can use the values from the $\gamma = 0.58^4$ case in Fig. 2, namely $\nu \gtrsim 2$ for maxima and $\nu \gtrsim 1$ for filamentary saddles.

### 2.1.2 Spatial correlations of GRF extrema

The first information about the relative spatial distribution of extrema is given by their two-point correlation function

$$\xi_{\text{ab}}(r) \equiv \frac{\langle C_a(\mathbf{x}) C_b(\mathbf{x}+\mathbf{r}) \rangle}{\langle C_a(\mathbf{x}) \rangle \langle C_b(\mathbf{x}) \rangle} - 1, \qquad (7)$$

where $a, b$ designate either peak or saddle and $C_a(\mathbf{x}) = |\det(\rho_{ij})| \delta_p(\nabla\rho) \theta_a(\{\lambda_i\})$ is the localized number density of extrema as used in equation (4). Additionally, one can add extra (and possibly different) constraints on the two extrema, for instance specific heights. Here, we briefly describe the features of peak–peak and peak–saddle correlations that are of main importance for connectivity discussion and refer the reader to Appendix A for further mathematical and computational details.

<footnote n="1">¹Note that in this paper, we will use the word ‘extrema’ interchangeably with ‘critical points’ to refer to all points satisfying the condition of zero gradient meaning both maxima, minima, and saddle points.</footnote>

<footnote n="2">²In this paper, ‘normalized’ means rescaled by its variance so that fields of zero mean and unit variance are considered.</footnote>

<footnote n="3">³For example for 2D saddles, the standard deviation from the mean is given by $\sqrt{1-2\gamma^2/3}$.</footnote>

<footnote n="4">⁴$\gamma = 0.58$ corresponds to power-law power spectra with spectral index $n_s = -1$ (respectively $-2$) in 2D (respectively 3D), and to a redshift zero $\Lambda$CDM Universe on scales about 1 Mpc.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Differential number of maxima and filamentary saddle points with height $v$ in the spherical volume of radius $R_p$ for 2D (left) and 3D (right) GRFs with $\gamma = 0.58$. Shaded areas mark the regions of high ‘rare’ maxima and filamentary saddles (above $1\sigma$ from the mean).</figcaption></figure>

For peak–peak correlations, the fundamental result of Kaiser (1984) showed that high-density peaks at large separations are positively correlated in excess of the mean correlation of the field (i.e. positively biased). A more detailed study of Baldauf et al. (2016) demonstrated that the peaks are anticorrelated (i.e. avoid each other) at small separations $R \lesssim R_p$, while an enhancement of the correlation between pairs of peaks at large separations is only pronounced if both peaks are of similar, and high, height $v$. A representation of these results is given in Fig. 3.

We shall focus on the configurations where the central peak is rare (i.e. has high $v_c$) while its neighbours have a distribution of heights $v_n$ above some threshold $v_*$, $v_n > v_*$. Fig. 3 shows that such a peak of GRF has a statistically well-defined exclusion neighbourhood where the probability of finding another peak is suppressed. In this paper, we use the extent of the exclusion region to define the size of a peak-patch around the peak, $R_{pp}$. When we consider only high neighbours above the threshold (as displayed on the right-hand panels of Fig. 3), we find that the high peaks demonstrate a sense of a ‘first layer’ of neighbours, the correlation function having a pronounced maximum therefore showing enhanced probability for neighbours to be at a particular distance. It is natural to identify the peak patch radius $R_{pp}(v_c, v_*)$ as the position where the peak–peak correlation function $\xi_{pk}$ shows a maximum i.e. when $\xi'_{pk}(r = R_{pp}) = 0$. From the point of view of the connectivity, it is these first neighbours that the filament bridges will connect to, and the length of such filaments will typically be of the order of $R_{pp}(v_c, v_*)$. When we study $R_{pp}(v_c, v_*)$ as function of $v_*$ at fixed $v_c$, we find that the distance to the next peak decreases as $v_*$ increases to reach $v_c$, which is the essence of Kaiser bias.<fnref n="5" />

The left-hand panels in Fig. 3 describe the case when we consider neighbouring peaks of all heights. The correlation function demonstrates that the extend of the exclusion region $R_{pp}(v_c)$ increases with the height $v_c$ of the central peak. The definition for the exact boundary of the patch in this case is somewhat ambiguous since the correlation function does not exhibit a well-defined maximum and no clear preferred position of the first neighbours is present. Still, it is clear from our numerical results that the height dependence of exclusion radius for rare peaks, $v_c > 2$, is roughly linear, with slope $\sim 0.2$, i.e. $R_{pp} \approx R_0 + v/5$ (in units of $R_p$), both in 2D and in 3D (in

particular one can track the width of the exclusion zone at half-min level of the correlation $\xi_{pk} = -0.5$). If we calibrate this relation at $v = 2$ by taking as the boundary the radius where the correlation turns over to zero (i.e. when its curvature is maximal), we have $R_{pp} \approx 1.1 + v/5$ in 2D and $R_{pp} \approx 0.9 + v/5$ in 3D. The accuracy of our correlation computations does not warrant higher precision than $\sim 0.03$ in these relations.

Fig. 4 shows the correlations between a central peak and its neighbouring saddle points, $\xi_{pk-sad}$, computed in the same way as the peak–peak correlations (only 3D plots are presented).

We will need these results later in the paper when we will count the filamentary saddle points in the vicinity of a peak as a handle on peak connectivity. For now, we note that the presence of a high peak affects the density of saddles only at $r < R_p$. Within $R_p$, a high peak attracts saddles which are also high (even more so when the heights of peak and saddle get closer), but repels the lower ones. This means that when all saddles are counted (left-hand panel), a very high peak ($v > 2.5 - 3$) predominantly repels the saddles (since most saddles will be notably lower than the peak in this case) from their vicinity, but less prominent peaks ($v \sim 2$) have actually an enhanced density of saddles near them. For $r > R_p$, the density of saddles becomes close to the mean value.

## 2.2 Morse theory and the skeleton picture

In order to compute connectivity of the cosmic web, we must properly define and extract it as a set of contiguous filamentary branches.

### 2.2.1 Ridge extractor algorithms

Over the years, several methods have been developed for this purpose (see, for instance Libeskind et al. 2018, for a comparison of some of those cosmic web classification schemes). Two distinct hypotheses can be broadly identified at the heart of these approaches. One can either use the geometrical information contained in the *local* gradient and the Hessian of the density or potential field (e.g. Novikov et al. 2006; Aragón-Calvo et al. 2007b,a; Hahn et al. 2007a,b; Sousbie et al. 2008a,b; Forero-Romero et al. 2009; Bond, Strauss & Cen 2010b,a), or the topology and connectivity of the density field using the watershed transform (Aragón-Calvo, van de Weygaert & Jones 2010), and/or Morse theory (e.g. Colombi et al. 2000; Sousbie et al. 2008a; Sousbie 2011). In the watershed category, Sousbie, Colombi & Pichon (2009) presented for instance

<footnote n="5">⁵This trend reverses with further increase of $v_* > v_c$ when $R_{pp}$ starts to increase, since the distant, now higher, peak begins to dictate how far its neighbours can be.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Top row: 2D GRF with $n_s = -1$, bottom row: 3D ΛCDM GRF smoothed on 5 Mpc $h^{-1}$. Right-hand panels: peak-peak correlation function when the central peak has a specified height (from 2 to 4 as labelled), while the height of its neighbour is unrestricted. The separation is in units of the mean distance between peaks $R_p$. Left-hand panels: correlation function of a peak with fixed height $v_c = 2.5$ (rare peak) and the second peak with the height exceeding $v_*$, for various thresholds $v_*$ between 0 and 2.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Same as Fig. 3 when the distant critical point is a filament-type saddle point.</figcaption></figure>

a method to compute the full hierarchy of the critical subsets of a given density field.

### 2.2.2 The skeleton of smoothed fields

Formally, the skeleton of a continuous field $\rho$ can be defined in the context of Morse theory (Jost 2008) as the set of critical lines (i.e. field lines which go through critical points) connecting the saddle points and the local maxima of that field while departing from the saddle along the first eigenvector of the curvature tensor (Novikov et al. 2006). Peak (respectively void) patches of the density field can then be defined (see e.g. Sousbie et al. 2009) as the set of points converging to a specific local maximum (respectively minimum), while following field lines in the direction (respectively the opposite direction) of the gradient. One expects filaments to lie at the intersection of such void patches. In particular, note that with this definition, one filament linking two maxima together necessarily

passes through one and only one critical point which has to be a filament-type saddle point. Morse theory formalizes this intuitive construct by segmenting space in respectively the void patches, the walls, the filaments, and the peaks of the cosmic web. The main advantage of watershedding algorithms is to provide a fully connected set of critical lines, but two shortcomings remain: (i) the Morse formalism only truly applies to Morse functions that are smooth and non-degenerate (the Hessian is assumed not to be zero at critical points). When dealing with discrete fields, one needs to smooth the density field before computing its skeleton and therefore introduce a smoothing length. (ii) More dramatically, all watershed algorithms implemented on discretized meshes will overproduce filaments, because the segmentation is carried at finite resolution, where the underlying Morse theory is not satisfied.

These unexpected shortcomings are of prime importance when measuring cosmic connectivity as they correspond to the appearance of bifurcation points that need to be identified and properly
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Illustration of void patches A–D whose intersections define the filaments. This figure can be interpreted in 2D, or in 3D as a cross-section perpendicular to the axis of the blue filament. In the former case, the bifurcation points in red where two critical lines seem to merge away from the critical point in blue is a resolution driven artefact (Pogosyan et al. 2009). In the latter case, bifurcation lines in red are also artefact due to the resolution and would wrongly be identified as filaments (respectively at the intersection of voids A–C–D and B–C–D, see Sousbie 2011, for details).</figcaption></figure>

accounted for. Indeed, when the field is constant over an extended region – a generic situation when cloud-in-cell sampling is applied on a grid – even after smoothing(!), the Morse condition is not satisfied anymore and there can be an ambiguity on the exact location of the filamentary segments. This situation is ubiquitous as void patches cannot all be convex: typically some of them will be squashed by their stronger neighbours, hence their boundaries will be almost tangent until reaching the filament. At finite resolution, the measured intersection line occurs before the real location of the filament. This can be seen on the cross section shown in Fig. 5. In 2D, it will occur on locations different from minima, saddles, and maxima, where an under-resolved critical line splits in two (whereas Morse theory states that this can happen only at critical points). In 3D, it leads to pairs of spurious one-dimensional lines, wrongly identified as filaments by the watershedding algorithm. Those were named bifurcation points and lines respectively by Pogosyan et al. (2009).

### 2.2.3 Discrete tracers and topological persistence

In order to resolve the above-mentioned issues, the Discrete Persistent Structure Extractor algorithm (DISPERSE, Sousbie 2011)<fnref n="6" /> implementing discrete Morse theory (Forman 2002) was developed. This geometric three-dimensional ridge finder allows for a scale- and parameter-free topologically consistent extraction of all different components of the fields. It can be applied either directly to the underlying DTFE density of points<fnref n="7" /> or to a regular mesh, from which the discrete Morse–Smale complex of the density function is computed (by assigning a height to all simplices of

the complex: vertices, edges, faces, volumes, etc.). The grid is partitioned according to the discrete gradient flow of density into an ensemble of critical sets – volumes, surfaces, curves, and points corresponding to the voids, walls, filaments, and clusters within the cosmic web, the so-called ascending 3-2-1- and 0-manifolds of the discrete flow, respectively. The code uses persistence ratio (the relative height of connected critical points as a measure of the significance of their topological connections) to filter out filaments dominated by the noise.

Since the construct satisfies discrete Morse theory (Robins 2000), all critical lines only split at critical points as they should, but multiple critical lines may overlap one another up to some bifurcation point (whose position depends on the sampling and smoothing). Filaments are not counted twice, only two branches emanate from each saddle point. The one minor drawback of the algorithm is that there is as-of-today no complete theory relating persistence to smoothing, so one must rely on calibration when comparing measurements to predictions for random fields smoothed on a given scale (see Appendix C). From a practical perspective, DISPERSE keeps track via a structure stored at each critical points of how many saddle points are connected to it. It also stores the segments where bifurcation occurs.

## 3 CONNECTIVITY OF GRF

Let us now turn to the theory of the connectivity of GRF as measured in GRF realizations or predicted from first principles.

### 3.1 Connectivity of peak patches

#### 3.1.1 Numerical simulations

To get two-dimensional maps, we first generate 100 realizations of a GRF with power-law power spectrum $P(k) \propto k^{n_s}$ on a $2048^2$ grid. We then smooth those maps with a Gaussian kernel on 8 pixels. The skeleton of each map is extracted using DISPERSE with a cut in persistence that depends on the power spectrum and is set so that the number of peaks in the maps matches the expected number of peaks $4096(n_s + 4)/\sqrt{3\pi^2}$ with better than 0.5 per cent accuracy. Adopted persistence cuts are given in Table C1. The spurious low-persistence peaks are due to the sampling of the Gaussian field as illustrated in Fig. C1.

The skeleton obtained with DISPERSE is further smoothed following a prescription described in Sousbie et al. (2009) which ensures that the number and positions of extrema are fixed under the smoothing operation. Thus, smoothing only straightens the skeleton segments and shifts bifurcation points but preserves topology and therefore connectivity.

A similar procedure is adopted in the three-dimensional case for which we generate 20 realizations of a GRF on a $256^3$ grid with power-law power spectrum $P(k) \propto k^{n_s}$ and further smoothed on 4 pixels. We also use the total number of 3D peaks to choose the persistence cut. In this case, the persistence cuts are found to be $p_{min} = \{4.5, 6, 9, 12\}\sigma_0/1000$ for $n_s = \{-3, -2, -1, 0\}$. Note that those runs are generated for only one reference spectral index ($-1$ in 2D and $-2$ in 3D). For other spectral indices, we will use only 10 runs because we restrict ourselves to marginals in those cases and statistics is therefore sufficient. Once the skeletons are computed, the number of saddle points connected to each peak are counted and a histogram is computed in order to get the probability density function (PDF) of the connectivity for different power spectra. We

<footnote n="6">⁶The code DISPERSE is publicly available http://www.iap.fr/users/sousbie/disperse.html</footnote>

<footnote n="7">⁷This paper will restrict itself to regular cubic meshes, but all results apply directly to discrete point-like sample via DTFE (see e.g. Malavasi et al. 2017; Kraljic et al. 2018; Laigle et al. 2018, for examples of implementation on observed surveys).</footnote>
also attach to each node the value of its height so that the joint statistics of nodes' connectivity and height can be investigated.

### 3.1.2 Peak connectivity of 2D GRF

The statistics of connectivity in 2D is shown on Fig. 6. We first notice on the right-hand panel that the mean number of connectors found is four with no detectable dependence on the slope of the power spectrum. On average, cosmic web's nodes therefore follow a square lattice. This result is in agreement with expectations from extrema counts as will be described below. The left-hand panel displays the full PDF of 2D connectivity. Most peaks have a connectivity between 2 and 6, but a long tail pervades with (rare) peaks having a connectivity as large as 10. The shape of the power spectrum also matters with more negative spectrum being more skewed and positive spectrum being more symmetric around the mean.

We anticipate that the number of connectors should also depend on peak height. Hence, the joint PDF, $P(\kappa, v)$, of the peak's number of connections $\kappa$ and height $v$ is also measured. This is the central result of this paper. Fig. 7 displays this joint PDF as sets of $v = \text{const}$ and $\kappa = \text{const}$ slices for a fixed spectral index $n_s = -1$. As a complement, Fig. 8 describes the relation between connectivity and peak height by plotting the conditional average number of connectors at fixed height, $\langle \kappa | v \rangle$, and conversely, the mean height for a given number of connectors, $\langle v | \kappa \rangle$, this time for several spectral indices.

We find that the rarer the peak, the higher the connectivity. This was expected since near a high contrast peak all eigenvalues tend to become equal (Pichon & Bernardeau 1999). Therefore, all incoming directions become possible. For low-density peaks, the number of connectors is found to be close to $\kappa = 3$ rising almost linearly for positive contrasts, with very high peaks ($v \approx 4$) reaching a mean connectivity of about six. Note that $v$ is an absolute density threshold, so that low-density maxima are predominantly peaks inside larger underdense (i.e. void) regions. Thus, we find that peaks in voids have a reduced number of webbing connectors with surrounding structures. We note that the correlations between peak's height and connectivity is spectrum dependent.

### 3.1.3 Peak connectivity of 3D GRF

In three dimensions, we follow the same path as in the 2D case. The mean and full PDF of 3D connectivity is shown in Fig. 9. On average, a 3D peak is connected to $\langle \kappa \rangle \approx 6$ saddle points with no detectable dependence on the spectral index. The mean 3D cosmic web is therefore close to a cubic lattice on average. Neglecting the special case of $n_s = -3$, the full statistics of the number of connectors depends only marginally on the shape of the power spectrum, peaking at a connectivity of five and extending to quite high values of the order of 20 for the rarer objects. Note however that the statistics in our 3D measurements is not high enough and we may have some (small but non vanishing) boundary effects given the small volume of the 3D maps we use in this paper. It is however computationally expensive to significantly increase the volume of each map because of the scalability limitations of DISPERSE.

Fig. 10 then displays the joint PDF of peak's height and connectivity $P(\kappa, v)$ and Fig. 11 shows the corresponding marginals namely the mean connectivity given height $\langle \kappa | v \rangle$ and mean height given connectivity $\langle v | \kappa \rangle$. Similar to the 2D case, it is shown that the mean connectivity increases with peak height from three in underdense regions to more than 10 connectors on average for the most massive

peaks ($v \gtrsim 4$). The interplay between peak's connectivity and height does depend on the slope of the power spectrum: a peak of a given connectivity tends to be systematically higher for high (positive) values of the spectral index and smaller for more red-tilted spectra.

## 3.2 Multiplicity of a peak, branches, and bifurcations

Up to now, we have focused on the connectivity of a peak patch, defined as the number of saddle points at the patch boundary through which the skeleton in the patch is connected to its neighbours. This definition is topological and the connecting skeletons lines never formally intersect away of the critical points. However, from the physical perspective such lines may, and often do, pass parallel to each other over significant distance, thus representing the vicinity of the same dense ridge. Indeed, infinitesimally close to the peak, the field can in general be described as a quadratic surface, thus (unless degenerate) has one leading eigen-direction of the curvature tensor, which corresponds to exactly two locally defined ridges, vicinity of which is tracked by every skeleton line emanating from the peak. Examples of such situation can be seen in Fig. 1.

To model physical overdense filaments, one would like to join close skeleton lines into a single object, until the position that will now be identified with a bifurcation point. Such procedure, and, thus, the position of bifurcation points is smoothing scale dependent but can be made robust, since the divergence of the two nearly parallel skeleton lines is usually exponentially quick. Importantly, we do not count bifurcation points within one smoothing radius from the peak, treating the skeleton branches that diverge each other so immediately as always distinct. It is exactly them that are counted in the multiplicity of the peak. This procedure is incorporated into DISPERSE, where the implemented smoothing of the skeleton preserves the positions of extrema (and thus global connectivity properties), but merges nearby skeleton segments and defines the bifurcation points of such 'physical' skeleton.

Let us now count the actual number of filaments incident onto a given maximum. We call this measure multiplicity of the peak, $\mu$. Clearly, this local (intra patch) measure is in fact equal to the number of connecting saddle points minus the number of bifurcation points within the peak patch

$$ \mu = \kappa - n_{\text{bifurcation}} \qquad (8) $$

Fig. 12 displays the corresponding PDF $P(\mu)$ obtained in our simulations, to be contrasted to Fig. 6 which represents the same distribution for $\kappa$. The distribution of multiplicity is almost symmetric, centred at $\langle \mu \rangle \approx 3$ in 2D and $\langle \mu \rangle \approx 4$ in 3D, and does not extend to numbers as large as for the $\kappa$ distribution. Indeed, we have observed no peaks with multiplicity exceeding 6 in 2D and 11 in 3D. A small dependence with the spectral index is found in this case ($\langle \mu \rangle$ varying from 3.85 to 4.1 for $n_s$ between $-3$ and $0$). Hence, typically, four filaments (three in 2D) locally branch out from a peak and later bifurcate in order to connect to eventually six neighbouring nodes. This mean picture obviously varies from one node of the cosmic web to the other.

In particular, the rareness (or height) of the peak changes the mean multiplicity. To investigate this effect, Fig. 13 shows the joint distribution of peak's multiplicity and height. As expected, the larger the multiplicity, the higher the peak. For instance, 3D peaks of multiplicity $\mu = 2, 3, 4$, and $5$ are expected to have increasing height on average: $\langle v | \mu \rangle \approx 1.07, 1.15, 1.38$, and $1.68 (\pm 0.01)$ in the case where $n_s = -2$ as displayed on that figure. As can be seen in Fig. 14, this relation is almost linear with slope about 0.2 almost
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Left-hand panel: PDF of the number of connectors in 2D evaluated from 10 GRF maps with spectral index from +1 (blue) to -2 (red). Error bars represent the error on the mean. Right-hand panel: corresponding mean connectivity as a function of the spectral index. The mean number of connectors seems to be independent on the power spectrum and equal to four.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Joint PDF of the number of connectors and the peak height estimated from 100 realizations of a 2D GRF with spectral index $n_s = -1$. Slices of the PDF at fixed connectivity and fixed height are, respectively, displayed on the left- and right-hand panels.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Left-hand panel: mean number of connectors as a function of the peak height in 2D evaluated from 10 GRF maps with different spectral indices as labelled. Right-hand panel: same as the left-hand panel for the mean peak height as a function of the number of connectors.</figcaption></figure>

independent from the spectral index and y-intercept varying from ~0.3 for $n_s = -3$ to ~0.8 for $n_s = 0$.

Pogosyan et al. (2009) identified the point process allowing us to build a proxy for the number of bifurcation points which were found to lie in the vicinity of where the skeleton splits. It relied on the degenerate condition of equal eigenvalues of the Hessian (i.e. local

isotropy), so that the next order in the Taylor expansion defines the direction of the splitting branches. Since we expect that the third eigenvalue remains distinct, the bifurcation will generically be co-planar. Appendix D quantifies the geometry of the bifurcation as a function of the relative strength of components of the gradient and that of the third-order derivative of the field.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Left-hand panel: PDF of the number of connectors in 3D evaluated from 10 GRF maps with spectral index from 0 (blue) to -3 (red). Error bars represent the error on the mean. Right-hand panel: corresponding mean connectivity as a function of the spectral index. The mean number of connectors seems to be independent of the power spectrum and equal to four.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Left: the $κ$ = const slices of $P(κ, ν)$ from 20 realizations of 3D $γ = 0.72 \text{ GRF}_{30}^{256}$ field, as functions of contrast $ν$; and right: the $ν$ = const slices of $P(κ, ν)$ for the same ensemble.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Left-hand panel: mean number of connectors as a function of the peak height in 3D evaluated from 10 GRF maps with different spectral indices as labelled. Right-hand panel: same as the left-hand panel for the mean peak height as a function of the number of connectors.</figcaption></figure>

## 3.3 Theory of the peak connectivity

In this section, we develop a theoretical framework to predict GRF connectivity from first principles. This is a challenging question as connectivity is by nature a global property and is difficult to catch from purely local considerations. Yet we will show that the mean connectivity in the field can be easily predicted. We will also develop a simplified theoretical approach that explains the increase

of connectivity with increase of rare peaks height and even predicts reasonably well quantitative aspects of this behaviour.

### 3.3.1 Global analysis

The theoretical estimate of the average number of connections of a peak comes simply from counting the total number of saddles and peaks in the field. Since every saddle point is connected to two and
<figure><img src="image_15.png" /><figcaption>**Figure 12.** The PDF of the multiplicity $\mu$ of peaks in 10 GRF with various spectral indices. This quantity reflects the number of incident skeleton filaments at the node. Left-hand panel: 2D estimate from $10\ 2048^2$ GRF smoothed on 8 pixels, $\langle\mu\rangle \approx 3$ for all spectral slopes. Right-hand panel: same quantity in 3D for $10\ 256^3$ smoothed on 4 pixels, here $\langle\mu\rangle \approx 4$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** PDF of the multiplicity $\mu$ of peaks i.e. their connectivity corrected for the number of bifurcations within the peak patch. Left-hand panel: results obtained from 20 2D $2048^2$GRF with spectral index $n_s = -1$; the corresponding mean is $\langle\mu\rangle = 3$ for all values of $n_s$. Right-hand panel: same quantity in 3D for $20\ 256^3$ GRF. Here, $\langle\mu\rangle = 4$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Mean multiplicity as a function of the peak height from 10 2D $2048^2$GRF (left) and 10 3D $256^3$ GRF (right) with spectral index as labelled.</figcaption></figure>

only two peaks and each peak to $\kappa$ saddles, the mean connectivity of a peak simply reads

$$ \langle \kappa \rangle = \frac{2\bar{n}_{\text{sad}}}{\bar{n}_{\text{max}}} \qquad (9) $$

which for GRF, from equations (1)and (2), is exactly 4 in 2D and $\approx 6.11$ in 3D. This means that in 2D, the connected skeleton on average is similar to a cubic lattice, while this is not exactly the case in 3D and some (small) ‘crystallographic’ defects should appear. However, we have not detected the departure from $\langle\kappa\rangle = 6$ having

limited statistics in our 3D measurements. Note that the connectivity can also be predicted for GRF in higher dimensions, it can be shown to follow

$$ \kappa = 2d + \left( \frac{2d - 4}{7} \right)^{7/4} \qquad (10) $$

at least for dimension $d$ between 2 and 11 (see Appendix E for a quick derivation of this relation).

Equation (9) is general and in particular not restricted to GRF. It can therefore be computed in the weakly non-linear regime using
the predictions for extrema counts given, for instance, in Gay et al. (2012) by means of a Gram–Charlier expansion. In 3D, the result at first order in non-Gaussianity reads

$$ \langle \kappa \rangle = \kappa^G \left( 1 + \sum_{i \ge 1} \kappa^{(i)} \sigma_0^i \right), \qquad (11) $$

with the first two terms given by

$$ \kappa^G = 2 \times \frac{29\sqrt{3} + 18\sqrt{2}}{29\sqrt{3} - 18\sqrt{2}} \approx 6.11, \qquad (12) $$

and

$$ \kappa^{(1)} = \frac{4\sqrt{3}}{35\sqrt{\pi}\sigma_0} \left( 8 \langle J_1^3 \rangle - 10 \langle J_1 J_2 \rangle - 21 \langle q^2 J_1 \rangle \right). \qquad (13) $$

Equation (13) is written in terms of the extended skewness parameters $\langle J_1^3 \rangle / \sigma_0$, $\langle J_1 J_2 \rangle / \sigma_0$, and $\langle q^2 J_1 \rangle / \sigma_0$ where $\sigma_0^2 = \langle |\nabla \rho|^2 \rangle$ is the variance of the density gradient, $\sigma_0^2 = \langle (\Delta \rho)^2 \rangle$ the variance of the density Hessian, $\sigma_1^2 q^2 = |\nabla \rho|^2$ the modulus square of the gradient, $\sigma_2 J_1 = \text{Tr} \rho_{ij}$ the trace of the Hessian matrix, and $\sigma_2 J_2 = 3/2 \text{Tr}(\rho_{ij}^2) - 1/2 (\text{Tr} \rho_{ij})^2$. In the cosmological context, these extended skewness parameters are constant in time at tree order in perturbation theory – similarly to $S_3 = \langle \rho^3 \rangle / \sigma_0^4$ – and depend on the slope of the underlying power spectrum.<fnref n="8" /> Note that Section 4.2 and in particular Fig. 20 will compute equation (13) in the cosmological context. Going to second order in the variance then requires to compute extended kurtosis parameters appearing in the next order term of the Gram–Charlier expansion described in Gay et al. (2012).

Note that as usual, equation (11) can be interpreted as a measure of the temporal and scale evolution of the connectivity given that $\sigma$ is a measure of the amplitude of non-linearities and depends on both time and scale. When conducting a dark energy experiment based on the cosmic evolution of the connectivity, the leading contribution will come from $\sigma_0 = D(z)\sigma_0(z=0)$ in equation (11), where $\sigma_0^2(z=0)$ is the linear variance of the density field at redshift zero, $D(z)$ is the growth rate given by

$$ D(z) = \frac{5\Omega_m H_0^2}{2} H(a) \int_0^a \frac{da'}{a'^3 H^3(a')}, \qquad (14) $$

where $H(a)=H_0\sqrt{\Omega_m/a^3+\Omega_\Lambda\exp[3\int_0^z(1+w(z'))/(1+z')dz']}$, with $\Omega_m$, $\Omega_\Lambda$, and $H_0$, respectively, the DM and dark energy densities and the Hubble constant at redshift zero, $a \equiv 1/(1+z)$ the expansion factor and $w(z) = w_0 + w_1/(1+z)$ the parametrized equation of state of dark energy. Since $\kappa^{(1)}$ is a number which can be predicted from cosmological perturbation theory, the cosmic variation of $\langle \kappa \rangle$ puts constraints on ($w_0, w_1$). In practice, one also needs to control the redshift evolution of the tracer threshold (e.g. the luminosity cut), since the mass function of haloes also evolves with redshift.

In turn, in 2D we can also easily obtain the Gaussian

$$ \kappa^G = 4, \qquad (15) $$

and first non-Gaussian-order contribution to the global connectivity

$$ \kappa^{(1)} = \frac{2}{27} \sqrt{\frac{6}{\pi}} \left( 5 \frac{\langle J_1^3 \rangle}{\sigma_0} - 6 \frac{\langle J_1 J_2 \rangle}{\sigma_0} - 18 \frac{\langle q^2 J_1 \rangle}{\sigma_0} \right), \qquad (16) $$

where $\kappa^{(1)}$ is again independent of $\sigma_0$ at tree order in perturbation theory. In equation (16) the definition of $J_2$ is now changed to $\sigma_2^2 J_2 = 2 \text{Tr}(\rho_{ij}^2) - (\text{Tr} \rho_{ij})^2 = (\lambda_1 - \lambda_2)^2$, while $J_1$ an $q^2$ have the same definition as in the 3D case. Again, the first non-Gaussian correction scales like $\sigma_0$ and is therefore a direct tracer of the growth of structure.

### 3.3.2 Connectivity as a function of peak height

A simple idea to estimate the change of connectivity with the peak height $v_c$ is to count the number of saddles in the volume around the peak through which the filamentary connectors to the neighbouring peaks pass. The main issue is therefore to quantify the typical volume where bridging saddles are. We suggest that this should be the volume of the peak patch where there is a suppressed probability for other peaks to be present. Then, the average saddle count, conditional on the presence of the peak of height $v_c$ at the centre of the volume, and therefore the connectivity of the peak, is

$$ \kappa(\nu) = \bar{n}_{\text{sad}} \int_{0}^{R_{\text{pp}}} d^D r (1 + \xi_{\text{pk-sad}}(r, \nu)), \qquad (17) $$

where $\xi_{\text{pk-sad}}(r, \nu)$ is the correlation between peaks of height $\nu$ and saddles of any height, and $\bar{n}_{\text{sad}}$ is the mean number density of saddle points.

To understand what connectivity equation (17) predicts we refer back to Section 2.1.2 for necessary ingredients, namely the peak-saddle correlation function and the size of peak patch $R_{pp}$. The resulting dependence on $v_c$ is the result of the competition between two effects: on the one hand, the suppression of saddles near the peak as the height of this central peak increases and, on the other hand, the increase of the peak-patch volume over which saddles need to be counted. The first effect is practically contained within $r < R_p$ radius, while the peak-patch radius $R_{pp}$ increases linearly with $v_c$ beyond $R_p$, including the outer region where the density of saddles becomes equal to the average one. Fig. 15 shows how the conditional saddle counts accumulate with the radius. The figure demonstrates that with $R_{pp}$ taken to be $R_{pp} = 1.1 + v_c/5$ in 2D and $R_{pp} = 0.9 + v_c/5$ in 3D as discussed in Section 2.1.2, the volume effect dominates the local (near peak) saddle suppression and on balance the theory predicts a growth of connectivity with peak height. Perhaps surprisingly, our estimate for $R_{pp}$ gives even quantitative values for $\kappa(v_c)$ in the interval $v_c = 1.5-4$ close to the 2D and 3D results measured in Figs 8 and 11, as Fig. 16 attests. As a note of caution, one should not overestimate the quantitative accuracy of this prediction given the precision of our measurements and simplicity of the theoretical model.

## 3.4 Local analysis of peak connectivity, multiplicity, and bifurcations

The complexity of predicting peak connectivity arises from the non-local topological nature of the global skeleton. Remarkably, the study of the number of real physical overdense filaments emanating from a peak lends itself more readily to the rigorous local analysis. A way to proceed is to compute the total number of 2D maxima at the surface of the sphere centred on the central peak and study this number as the function of the radius of the sphere. This is an approximation to the number of filaments crossing the sphere assuming that for high peaks, filaments are sufficiently ‘stiff’.

This approach gives a detailed picture of how many filaments leave the immediate peak volume, where they bifurcate at larger

<footnote n="8">⁸They correspond to isotropic moments of the underlying bispectrum, see Gay et al. (2012).</footnote>
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Expected total number of saddles in spheres of increasing radius $r$ in 2D ($n_s = 0$) and 3D ($n_s = -2$) GRFs around a central peak of height $v$ as labelled. The dotted lines represent the size of the corresponding peak patch $r = R_{pp}$ according to $R_{pp} = 1.1 + v/5$ in 2D and $R_{pp} = 0.9 + v/5$ in 3D, as described in Section 2.1.2.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Mean connectivity of peaks of height $v$ in 2D (blue) and 3D (orange) GRF. In 2D, we consider a spectral index $n_s = 0$, while in 3D, we choose for this figure $n_s = -2$. The error bars display the measured relation shown in Figs 8 and 11. The shaded areas represent the theoretical predictions given by equation (17) given that $R_{pp} \approx 1.1 + v_c/5 \pm 0.03$ in 2D and $R_{pp} \approx 0.9 + v_c/5 \pm 0.03$ in 3D as justified in Section 3.3.2.</figcaption></figure>

radii (as the count increases) and, coupled with an estimate for the size of peak patch from Section 2.1.2, gives an alternative way to predict the global connectivity of the central peak. Being a straightforward conditional maxima count, albeit on 2D spherical sections of the underlying 3D field, this analysis can be formulated via the joint distribution of the field values and its derivatives in two different locations: the central peak and the position on a sphere around it. We refer the reader to Appendix B for the technical details of extrema counts on a curved spherical surface, and only describe the results here in the main text.

Figs 17 and 18 demonstrate the main results for 2D and 3D, respectively. The number of filaments crossing the sphere of radius $r$, $N_{\text{fil}}(r)$, is shown in the left-hand panels of Figs 17 and 18.

We find that the function $N_{\text{fil}}(r)$ is essentially universal, independent of the central peak height and the slope of the power spectrum. At small distances, peaks have an ellipsoidal shape with two filaments (that is to say two branches of a same physical filament aligned with the peak longest semi-axis), $N_{\text{fil}}(0) = 2$. The expected number of filaments then grows with distance to $\sim 3$ and $\sim 4$ in 3D at the scale of the smoothing length. This is what corresponds

to the mean multiplicity $\mu$ of the peaks, since our definition of the bifurcations points does not identify them as such within one smoothing radius from the peak.

At larger $r$, our algorithm detects bifurcations of the filaments that lead to further increase of the number of filaments with radius. At $r > R_p$, the effect of the central peak is negligible and the 2D maxima count simply follows the area of a spherical surface, growing linearly with $r$ in 2D and quadratically as $r^2$ in 3D. The predicted connectivity $\kappa$ is obtained by evaluating $N_{\text{fil}}$ at the peak-patch radius $R_{pp}$, $\langle \kappa | v_c \rangle = N_{\text{fil}}(R_{pp})$. The predicted connectivity does depend on the height because again $R_{pp}$ is a growing function of the peak height. For instance in 2D, using the typical peak-peak separation found in the previous section $R_{pp} = 1.3, 1.5, 1.7, 1.9 R_p$ for $v = 1, 2, 3$, and 4, the peak connectivity is found to go from $\sim 4$ to 6, in full agreement with our measurements in GRFs. A similar agreement is found in 3D where the connectivity ranges from 7 and 13 for peak of heights between $2\sigma$ and $4\sigma$ above the mean. This shows the consistency of the local approach with the other discussed measures of connectivity.

Both the method based on saddle counts in Section 3.3.2 and the one relying on maxima count on surrounding spheres developed here readily allow us to add a constraint on the overdensity of the filaments and therefore to study only physically prominent ones. The information obtained with the two techniques is complementary. Saddle counts reflect the change in the topological properties of the connected web. In this context, eliminating filamentary connections below a given density threshold changes the connectivity $\kappa$ of the web. Conversely, the approach of this section describes the local properties of filaments. It allows us to count them as a function of distance from the peak and track the filaments that are dense near the peak even if their density will drop below ‘detectability’ further away. So it tracks the peak branches (peak multiplicity), not necessarily the connecting segments (peak connectivity).

The right-hand panels of Figs 17 and 18, show the number of prominent filaments (e.g. defined as $v_f > 1$ in 2D and $v_f > 2$ in 3D) leaving the peak as function of the distance $r$. We consider two types of peaks, the very rare ones, $v_c = 4$, and lower and more frequent ones, $v_c = 3$. Focusing on the 3D case, one sees that a high-contrast peak is typically surrounded by three main filaments, while the patch of a less rare peak is dominated by two branches of presumably one dense filament. Notably, the most prominent filaments do not
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Left-hand panel: mean number of maxima on a circle of radius $r$ around a peak of height $v$ as a proxy for the number of filaments crossing that circle. The prediction is for a GRF with $k^0$ power spectrum filtered with a Gaussian kernel. The dashed line represents the case where the central point is random (not necessarily a peak position). Vertical dotted lines displays the peak patch size $R_{pp}(v)$. Right-hand panel: same statistics as in the left-hand panel with restriction that only filaments with density contrast $v_f$ exceeding 1(blue) or 2 (yellow) are counted.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 18.** Left-hand panel: mean number of density maxima on a sphere of radius $r$ around a peak of height $v$, as a function of $r$. The prediction is for a GRF with $\Lambda$CDM power spectrum filtered with a Gaussian kernel on 5 Mpc$^{-1}$. The black dashed line represents the case where the central point is random (not necessarily a peak position). Vertical dotted lines displays the peak patch size $R_{pp}(v)$. Right-hand panel: same statistics as in the left-hand panel with restriction that only filaments with density contrast $v_f$ exceeding 1(blue) or 2 (yellow) are counted.</figcaption></figure>

experience bifurcations (hence the yellow line almost constant on Figs 17 and 18), meaning that in any bifurcation only one child follows the dense parent. This situation pervades far from the central peak and all the way to the neighbouring peaks, therefore forming a low connectivity network. Of course, one expects the density in a filament to decrease as the distance from the peak increases. If it drops below the threshold, the filament is considered terminated, which is then evidenced by the drop of $N_{\text{fil}}(r|v_f > v_{f*}, v_c)$ with distance. Closer is the limit on the filament's density $v_{f*}$ to the peak height $v_c$, less probable it is for the filament to extend far out. Our theory therefore predicts the relation between the heights of the peaks and the density in the filaments that allows to build a fully connected web, i.e. to percolate from peak to peak along the filamentary network. The necessary criterium is to have $N_{\text{fil}}(r|v_f > v_{f*}, v_c > v_{c*}) > 2$ for all $r$. We will explore the applications of this criterium in future works.

# 4 COSMIC CONNECTION

Let us now turn to the connectivity of the evolved cosmic web as measured in DM simulations, first in 2D convergence maps, then in 3D for DM and dark haloes, respectively.

## 4.1 Connectivity of cosmic convergence maps

We use convergence maps for four galaxy source redshifts $z_s = 0.5, 1, 1.5, 2$, and $2.5$ from the MassiveNuS suite of $\Lambda$CDM simulations<fnref n="9" /> and corresponding to the case without massive neutrinos (Liu et al. 2017). The cosmological parameters are set to $h = 0.7, n_s = 0.97, \Omega_b = 0.046, \Omega_m = 0.3$, and $\sigma_8 = 0.8523$ concordant with current observations. These maps cover $12.25 \text{ deg}^2$ with $0.1$ arcmin resolution and were obtained from GADGET-2 runs using the LENSTOOLS PYTHON code (Petri 2016). We refer the reader to Liu et al. (2017) for further details on the simulations and pipeline to generate the convergence maps.

For the sake of simplicity, we select 10 of the 1000 realizations of the MassiveNuS suite and smooth the convergence fields with a Gaussian kernel on 8 pixels corresponding to a full width at half-maximum (FWHM) of 1.9 arcmin. We then measure the connectivity of each two-dimensional peak with DISPERSE. The persistence threshold is chosen so that the number of peaks found by DISPERSE is the same as the thoroughly tested MAP2EXT code

<footnote n="9">⁹Those maps are publicly available at columbia.lensing.org</footnote>
(Colombi et al. 2000; Pogosyan, Pichon & Gay 2011) which is described in Appendix C.

Fig. 19 shows the measured mean connectivity as a function of redshift and the corresponding full distribution on the left- and right-hand panels, respectively. As expected higher galaxy source redshifts are the closest to the Gaussian prediction. At lower redshifts, the field becomes more non-Gaussian and peaks are less connected. A similar picture also appears for the three-dimensional matter distribution as will be shown in the next section.

## 4.2 Connectivity of the 3D matter distribution

Let us now focus on the connectivity of the three-dimensional distribution of matter in the Universe by means of analytical arguments first and numerical simulations when non-linearities become important.

Fig. 20 shows the expected evolution of the connectivity of the three-dimensional distribution of matter in the Universe at first order in perturbation theory. To obtain this result, equation (11) is used with the extended skewness parameters analytically computed at tree order in cosmological perturbation theory, for a Gaussian smoothing and different power-law power spectra. Those leading order skewness-like terms involve hypergeometric functions (Gay et al. 2012).

Eventually, the combination of interest $\tilde{S}_3 = 8 \langle J_1^3 \rangle - 21 \langle q^2 J_1 \rangle + 10 \langle J_1 J_2 \rangle$ reads

$$ \tilde{S}_3 = 8\gamma(s_5^5 \mathcal{F}_5^5 + s_7^5 \mathcal{F}_5^7 + s_9^5 \mathcal{F}_9^9 + s_{11}^5 \mathcal{F}_9^{11}), \quad (18) $$

with the Gauss hypergeometric functions

$$ \mathcal{F}_5^5 = {}_2F_1 \frac{n_s + 5}{2} \frac{n_s + 5}{2} \frac{3}{2} \frac{1}{4}, $$

$$ \mathcal{F}_5^7 = {}_2F_1 \frac{n_s + 5}{2} \frac{n_s + 7}{2} \frac{3}{2} \frac{1}{4}, $$

$$ \mathcal{F}_9^9 = {}_2F_1 \frac{n_s + 9}{2} \frac{n_s + 9}{2} \frac{9}{2} \frac{1}{4}, $$

$$ \mathcal{F}_9^{11} = {}_2F_1 \frac{n_s + 9}{2} \frac{n_s + 11}{2} \frac{11}{2} \frac{1}{4}, $$

and coefficients

$$ s_5^5 = \frac{261212n_s^4 + 2882763n_s^3 + 6981958n_s^2 - 35748720n_s - 197568000}{25725n_s^2(n_s + 2)^2}, $$

$$ s_7^5 = \frac{51861600 + 6744654n_s - 1875041n_s^2 - 633302n_s^3 - 51932n_s^4}{10290n_s^2(n_s + 2)^2}, $$

$$ s_9^9 = - \frac{2392(n_s + 5)(n_s + 7)^2}{3087000}, $$

$$ s_{11}^9 = \frac{1276(n_s + 5)(n_s + 7)^2(n_s + 9)}{21609000}. $$

This result only holds for low values of the variance that is to say early times or large scales. At later times and smaller scales, this Gram-Charlier expansion is proved to converge very slowly breaking down for $\sigma_0 \lesssim 1$ and numerical simulations are needed in order to investigate accurately the time evolution of cosmic connectivity. To do so, we use 18 ΛCDM and 4 CDM simulations of a 50 Mpc$^{-1}$ periodic box with 256$^3$ particles evolved until redshift 0 with GADGET using as a cosmology $\Omega_M = 0.3$, $\Omega_\Lambda = 0.7$, and $\sigma_8 = 0.9$. As an illustration, Fig. 21 displays the skeleton of one of the ΛCDM run. From these persistent skeletons, we then measured the connectivity at various time slices and for a constant comoving smoothing $R \approx 0.8$ Mpc$^{-1}$ corresponding to 4 pixels. Once again,

the persistence cut is chosen at each time-step to match the number of peaks found by the MAP2EXT code.

Fig. 22 shows how connectivity and multiplicity evolve in these two types of simulations. As expected, connectivity decreases as voids grow and filaments merge with cosmic time. The presence of dark energy in the form of a cosmological constant changes the overall shape of this curve which is almost linear in the CDM case, while in a ΛCDM Universe, it is clearly convex which is expected since dark energy tends to slow down the growth of cosmic structure and therefore change the slope at late times. On the right-hand panel, multiplicity is shown to vary much less with time. We observe an almost constant shift between the CDM and ΛCDM Universes (with a slight decay in the ΛCDM case which might again be a manifestation of dark energy changing the slope).

Hence, the disconnection of filaments with cosmic time is a direct consequence of gravitational clustering, and can in fact also be understood ab initio in the framework of excursion set theory – which links the time evolution of cosmic structures to a random walk driven by cosmic variance/smoothing. With smoothing of the initial GRF as a proxy for its upcoming time evolution, one can identify a special smoothing scales corresponding to the coalescence of wall-type saddle and filament-type saddle points, which topologically correspond to the disappearance of a tunnel, or equivalently to two filaments merging into one.¹⁰ From the point of view of the local peak, its connectivity decreases by one when the two filaments merge. The relative change in the disconnection of filaments with cosmology reflects the same process, but is also in part driven by the role played by dark energy which will stretch and disconnect the filamentary structure through the increased expansion of voids via the cosmological constant.

Fig. 23 then shows the full PDF of the global connectivity in the ΛCDM simulations at various redshifts. At high redshift, the field is nearly Gaussian and the measured connectivity has a statistics close to the Gaussian prediction shown in Fig. 9 for a spectral index about $n_s \approx -2$ as expected for a ΛCDM spectrum at one megaparsec scale. Towards lower redshifts, the PDF gets more and more skewed and concentrated, the most likely connectivity is shifted towards lower values (3 instead of 4) and the tail is suppressed. The global connectivity therefore decreases with cosmic time and highly connected nodes become rarer and rarer due to filaments merging.

Given those findings, we anticipate that measuring the redshift evolution of cosmic connectivity could prove to be an interesting probe of dark energy.

## 4.3 Connectivity of dark matter haloes

We make use of the 43 million DM haloes detected at redshift zero in the Horizon 4π N-body simulation (Teyssier et al. 2009). This simulation contains 4096$^3$ DM particles distributed in a 2 h$^{-1}$Gpc periodic box and is characterized by the following ΛCDM cosmology: $\Omega_m = 0.24$, $\Omega_\Lambda = 0.76$, $n = 0.958$, $H_0 = 73$ km s$^{-1}$ Mpc$^{-1}$, and $\sigma_8 = 0.77$ within one standard deviation of WMAP3 results (Spergel et al. 2003). The initial conditions were evolved non-linearly down to redshift zero using the adaptive mesh refinement code RAMSES (Teyssier 2002), on a 4096$^3$ grid. The friend-of-friend algorithm (Huchra & Geller 1982) was then used over 18$^3$ overlapping subsets of the simulation with a linking length

<footnote>$^{10}$This is the analogue of maxima and filament-type saddle points coalescence identified by Hanami (2001) as slopping saddles tracing merger events.</footnote>
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Connectivity of convergence maps for source redshifts between 0.5 and 2.5 as labelled. The left-hand panel shows only the mean connectivity, while the right-hand panel displays the full distribution function.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** Evolution of the mean number of connectors in 3D as a function of the variance σ₀ – i.e. across cosmic time – for different spectral indices as labelled. The early evolution of the 3D connectivity is rather slow. We expect the subsequent highly non-linear evolution to modify this result for higher σ₀.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 21.** The skeleton computed by DISPERSE on DM particles in one of the ΛCDM simulations used in this paper using one particle in 40 at redshift zero.</figcaption></figure>

of 0.2 times the mean interparticle distance to define DM haloes. In this work, we only consider the 43 million haloes with more than 100 particles (the particle mass being 7.7 × 10⁹ M☉). The mass dynamical range of this simulation spans about five decades.

The dark halo catalogue was split into 250 subregions for which the skeleton was computed with DISPERSE using a persistence level of three. For each identified node, the mass of most massive dark halo in the vicinity was assigned to this node (out of five nearest neighbour; the answer shown not to depend strongly with this choice). The connectivity was then sampled as a function of mass and its PDF(κ, M) and mean value are shown on Fig. 24. As expected the connectivity at low mass is lower than for GRF, but grows significantly with the mass of the cluster. Its growth is well fitted by κ(M) ≈ 10/3log (M/10¹¹M☉).

Small peaks that typically host isolated galaxies tend to be fed by a small number of filaments. Typically those objects are accreting matter from two main filaments which can bring coherent angular momentum from large-scale flows. They will naturally be able to form disc-like galaxies with a spin aligned with galactic filaments. On the other hand, rare peaks are multiply connected so that galaxies in this cluster-like environment can hardly acquire angular momentum consistently (Dubois et al. 2012). The more isotropic accretion will favour the formation of elliptical galaxies in dense environments dominated by tidal forces. This relationship between number of filaments and spin was already highlighted in Pichon et al. (2011). We have also shown that three main filaments are expected to dominate, a result which was empirically found in numerical simulation (Prieto et al. 2015; González et al. 2017).

Using DISPERSE, one could also compute the skeleton of subhalo or (virtual or observed) galaxy catalogues, hence their connectivity as a function halo mass or colour and morphology. The result will depend on the chosen persistence threshold, which in turn will depend on the sparsity of the catalogue (hence its limiting mass or magnitude). This will be the topic of future work.

# 5 DISCUSSION AND CONCLUSIONS

This paper studied the connectivity of cosmic random fields in two and three dimensions, quantifying preliminary results presented in Pichon et al. (2010). Starting with generic GRFs, peaks were shown to be connected on average to 4 other peaks in 2D and ~6.1 in 3D (and 8.3, 10.7, and 13.2 in 4D, 5D, and 6D) independently of the shape of the power spectrum. The full statistics of the
<figure><img src="image_25.png" /><figcaption>**Figure 22.** Mean connectivity (left-hand panel) and multiplicity (right-hand panel) of the skeleton as a function of the expansion factor for ΛCDM and CDM simulations as labelled. As expected, the CDM simulation is essentially featureless, whereas the ΛCDM connectivity changes slope when the dark energy expansion kicks in.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 23.** PDF of cosmic connectivity at various redshifts as labelled is the ΛCDM simulations smoothed on a constant comoving length $R \approx 0.8\text{Mpc}^{-1}$.</figcaption></figure>

number of connectors was precisely characterized using DISPERSE: the distribution is skewed, picking around 3 in 2D (respectively 4 in 3D) with a strong tail that can extend up to quite large values 10 in 2D (respectively 20 in 3D). The PDF of the number of connectors was shown to depend on the power spectrum and more importantly on the height of the peak. Rare peaks (with high-density contrast) tend to be connected to numerous neighbours, the rarest ones having more than 7 (15 in 3D) connections on average. Interestingly, the overall number of connected neighbours does not correspond to the local number of filaments sticking out of a peak. Locally, the geometry of a peak is ellipsoidal with two filaments following the axis of minimal curvature. Further away those filaments split and bifurcation points appear. Typically, some of those bifurcations occur so close to the central peak that they are not detectable in a map with finite resolution. The corresponding multiplicity of a peak is given by its connectivity minus the measured number of bifurcations. For GRFs, the mean multiplicity was shown to be 3 in 2D and 4 in 3D. Here, the dependence with the power spectrum and the typical scatter is reduced. However, the dependence with peak height remains strong.

We have also developed a theoretical framework which explains the connectivity properties measured with the persistent skeleton. This theory relies on the statistical properties of saddle points in the vicinity of peaks and allowed us to understand precisely the mean

connectivity but also its dependence on peak height. A precise count of the number of filaments crossing a sphere centred on a peak was used to investigate how bifurcation points appear from the peak to the boundary of its peak patch. These calculations allowed us to study the peak's multiplicity and quantify how many filaments dominate (typically we find three dominant filaments for high peaks in 3D, and two for smaller ones).

Finally, the subsequent evolution of connectivity across cosmic time was also investigated. We developed predictions based on perturbation theory to probe the first stages of structure formation and relied on numerical simulations to confirm the predictions and extend them to a more non-linear regime. As expected, the evolution of cosmic connectivity depends primarily on the growth factor (and on generalized cumulants, see Gay et al. 2012) and therefore on cosmology. Also as expected, the non-linear evolution reduces the number of connections as filaments merge. Compared to alternative probes of dark energy, the connectivity (or the multiplicity) may prove to be a robust estimator given that it is a topological property that can be measured locally.

As an astrophysical application, we focused on galaxy lensing convergence maps and dark halo catalogues. As expected, the connectivity of kappa maps also decreases towards low redshifts, although the evolution is milder than for the three-dimensional fields as the lensing kernel tends to smooth the signal along the line of sight (see also Gouin et al. 2017; Codis et al. 2017, for an investigation of the conditional connectivity around clusters using strong lensing). In the near future, it should be possible to compute the 3D connectivity of HI density maps reconstructed from QSO absorption fluxes from the point spread function, WEAVE surveys, or intensity maps from CHIME, MeerKAT, ASKAP, MWA or HERA, and at some later stage, the E-ELT and SKA, respectively. We also computed the connectivity of dark haloes and found that it scales logarithmic with mass with a scaling going like 10/3.

Cosmic connectivity is not only of interest in the context of cosmology: on smaller scales, it is also paramount to understand galactic assembly. Indeed DM filaments have a baryonic continuation within dark haloes, which connect closely the cosmic environment to the galaxies within dark haloes. Beyond the number of connected filaments to a given galaxy, the mass load, geometry, and torques advected along filaments is also of interest. In this context, one should investigate in more details the small-scale connectivity of galaxies and dark haloes, feeding them with cold gas with
<figure><img src="image_27.png" /><figcaption>**Figure 24.** Left-hand panel: mean connectivity of dark haloes at redshift zero as a function of log mass as labelled for about one million nodes of the cosmic web. The mean cosmic connectivity here is found to be well approximated by the simple linear relation $\kappa(M) \approx 10/3\log(M/10^{11}M_{\odot})$. Right-hand panel: the corresponding PDF $P(\kappa) = P(\kappa|M)$.</figcaption></figure>

stratified angular momentum which allows them to reform stellar discs (Pichon et al. 2011). As a topological quantity, it would be of interest to quantify the (expected) robustness of the multiplicity and connectivity with respect to redshift space distortions, shot noise, photometric errors, etc. As shown in Section 3.3.1, non-linear gravitational coupling decreases the connectivity of dark haloes. This could in principle be quantified by predicting the rate of coalescence of wall-saddle and filament-saddle critical points, which reflect filaments coalescence.

Beyond astrophysics, the present theory of connectivity could prove to be of importance in the context of percolation theory. For instance, the percolation threshold can be explained in terms of the statistical properties of the connectivity of the relevant nodes (keeping track of heights of peaks). At a more abstract level, this work on GRF could be of interest to generically connect their properties to graph theory, since the skeleton provides a mapping from the density field to sets of connected vertices (see e.g. Coutinho et al. 2016). This is left for future work.

## ACKNOWLEDGEMENTS

We thank Stéphane Colombi, Yohan Dubois, Julien Devriendt, Simon Prunet, and Thierry Sousbie for numerous discussions about this project over the last 10 yr! We thank our collaborators and the staff at the CCRT for producing the horizon-4π simulation, D. Munro for freely distributing his YORICK programming language and opengl interface (available at http://yorick.sourceforge.net/), and the community of mathematica stackexchange for help. Many thanks to Stéphane Rouberol for customizing the Horizon cluster, hosted by the Institut d'Astrophysique de Paris, for our purposes, and Eric Pharabod for Fig. 6. This work is partially supported by the Spin(e) grant ANR-13-BS05-0005 of the French Agence Nationale de la Recherche. Special thanks also go to Jia Liu who kindly pointed us to the publicly available Columbia lensing group lensing maps. These maps were obtained with support from NSF grant AST-1210877 and NSF XSEDE allocation AST 140041. We thank New Mexico State University (USA) and Instituto de Astrofísica de Andalucía SCIC (Spain) for hosting the Skies and Universes site for cosmological simulation products. SC and CP thank Lena for hosting multiple visits, DP thanks the ILP for a senior visiting fellowship, while CP and DP thank CITA for hospitality during the completion of the project.

## REFERENCES

* Abbott T. et al., 2016, Phys. Rev. D, 94, 022001

* Adelman-McCarthy e. a., 2008, ApJS, 175, 297

* Alonso D., Eardley E., Peacock J. A., 2015, MNRAS, 447, 2683

* Alpaslan M. et al., 2016, MNRAS, 457, 2287

* Aragón-Calvo M. A., Jones B. J. T., van de Weygaert R., van der Hulst J. M., 2007a, A&A, 474, 315

* Aragón-Calvo M. A., van de Weygaert R., Jones B. J. T., van der Hulst J. M., 2007b, ApJ, 655, L5

* Aragón-Calvo M. A., van de Weygaert R., Jones B. J. T., 2010, MNRAS, 408, 2163

* Aubert D., Pichon C., 2007, MNRAS, 374, 877

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

* Bailin J., Steinmetz M., 2005, ApJ, 627, 647

* Baldauf T., Codis S., Desjacques V., Pichon C., 2016, MNRAS, 456, 3985

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

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

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

* Bond N. A., Strauss M. A., Cen R., 2010a, MNRAS, 406, 1609

* Bond N. A., Strauss M. A., Cen R., 2010b, MNRAS, 409, 156

* Chen Y.-C. et al., 2017, MNRAS, 466, 1880

* Codis S., Pichon C., Devriendt J., Slyz A., Pogosyan D., Dubois Y., Sousbie T., 2012, MNRAS, 427, 3320

* Codis S., Pichon C., Pogosyan D., Bernardeau F., Matsubara T., 2013, MNRAS, 435, 531

* Codis S., Pichon C., Pogosyan D., 2015, MNRAS, 452, 3369

* Codis S., Gavazzi R., Pichon C., Gouin C., 2017, A&A, 605, A80

* Cole e. a., 2005, MNRAS, 362, 505

* Colombi S., Pogosyan D., Souradeep T., 2000, Phys. Rev. Lett., 85, 5515

* Coutinho B. C., Hong S., Albrecht K., Dey A., Barabási A.-L., Torrey P., Vogelsberger M., Hernquist L., 2016, preprint (arXiv:1604.03236)

* Danovich M., Dekel A., Hahn O., Teyssier R., 2012, MNRAS, 422, 1732

* de Bernardis e. a., 2000, Nature, 404, 955

* de Lapparent V., Geller M. J., Huchra J. P., 1986, ApJ, 302, L1

* Dekel A. et al., 2009, Nature, 457, 451

* Dubois Y., Pichon C., Haehnelt M., Kimm T., Slyz A., Devriendt J., Pogosyan D., 2012, MNRAS, 423, 3616

* Dubois Y. et al., 2014, MNRAS, 444, 1453

* Forero-Romero J. E., Hoffman Y., Gottlöber S., Klypin A., Yepes G., 2009, MNRAS, 396, 1815

* Forman R., 2002, Sém. Lothar. Combin., 48, 35

* Gay C., Pichon C., Pogosyan D., 2012, Phys. Rev. D, 85, 023011

* González R. E., Prieto J., Padilla N., Jimenez R., 2017, MNRAS, 464, 4666

* Gott J. R., III, Dickinson M., Melott A. L., 1986, ApJ, 306, 341

* Gouin C., Gavazzi R., Codis S., Pichon C., Peirani S., Dubois Y., 2017, A&A, 605, A27

* Hahn O., Porciani C., Carollo C. M., Dekel A., 2007a, MNRAS, 375, 489
Hahn O., Carollo C. M., Porciani C., Dekel A., 2007b, MNRAS, 381, 41

Hanami H., 2001, MNRAS, 327, 721

Holmberg E., 1969, Arkiv Astron., 5, 305

Huchra J. P., Geller M. J., 1982, ApJ, 257, 423

Ibata R. A. et al., 2013, Nature, 493, 62

Jost J., 2008, Riemannian Geometry and Geometric Analysis, Fifth Edition.
Springer, Berlin

Kac M., 1943, Bull. Am. Math. Soc., 49, 938

Kaiser N., 1984, ApJ, 284, L9

Kang X., Wang P., 2015, ApJ, 813, 6

Kimm T., Devriendt J., Slyz A., Pichon C., Kassin S. A., Dubois Y., 2011, preprint (arXiv:1106.0538)

Klypin A., Shandarin S. F., 1993, ApJ, 413, 48

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

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

Libeskind N. I. et al., 2018, MNRAS, 473, 1195

Lin C. C., Mestel L., Shu F. H., 1965, ApJ, 142, 1431

Liu J., Bird S., Zorrilla Matilla J. M., Hill J. C., Haiman Z., Madhavacheril M. S., Petri A., Spergel D. N., 2017, preprint (arXiv:1711.10524)

Longuet-Higgins M. S., 1957, Phil. Trans. R. Soc. Lond. A, 249, 321

Lynden-Bell D., 1964, ApJ, 139, 1195

Malavasi N. et al., 2017, MNRAS, 465, 3817

Marcos-Caballero A., Fernández-Cobos R., Martínez-González E., Vielva P., 2016, J. Cosmol. Astropart. Phys., 4, 058

Matsubara T., 1994, ApJ, 434, L43

Mecke K. R., Buchert T., Wagner H., 1994, A&A, 288, 697

Metuki O., Libeskind N. I., Hoffman Y., Crain R. A., Theuns T., 2015, MNRAS, 446, 1458

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

Novikov D., Colombi S., Doré O., 2006, MNRAS, 366, 1201

Ocvirk P., Pichon C., Teyssier R., 2008, MNRAS, 390, 1326

Petri A., 2016, Astron. Comput., 17, 73

Pichon C., Bernardeau F., 1999, Astron. Astrophys., 343, 663

Pichon C., Gay C., Pogosyan D., Prunet S., Sousbie T., Colombi S., Slyz A., Devriendt J., 2010, in Alimi J.-M., Fuözfa A., eds, AIP Conf. Ser. Vol. 1241, The Skeleton: Connecting Large Scale Structures to Galaxy Formation. Am. Inst. Phys., New York, p. 1108

Pichon C., Pogosyan D., Kimm T., Slyz A., Devriendt J., Dubois Y., 2011, MNRAS, 418, 2493

Pogosyan D., Pichon C., Gay C., Prunet S., Cardoso J. F., Sousbie T., Colombi S., 2009, MNRAS, 396, 635

Pogosyan D., Pichon C., Gay C., 2011, Phys. Rev. D, 84, 083510

Pogosyan D., Codis S., Pichon C., 2016, in van de Weygaert R., Shandarin S., Saar E., Einasto J., eds, IAU Symp. Vol. 308, The Zeldovich Universe: Genesis and Growth of the Cosmic Web. Kluwer, Dordrecht, p. 61

Poudel A., Heinämäki P., Tempel E., Einasto M., Lietzen H., Nurmi P., 2017, A&A, 597, A86

Prieto J., Jimenez R., Haiman Z., González R. E., 2015, MNRAS, 452, 784

Rice S. O., 1945, Bell System Tech. J., 25, 46

Robins V., 2000, PhD thesis, Department of Applied Mathematics, University of Colorado

Schweizer F., 1982, ApJ, 252, 455

Scoccimarro R., Colombi S., Fry J. N., Frieman J. A., Hivon E., Melott A., 1998, ApJ, 496, 586

Slepian Z. et al., 2017, MNRAS, 468, 1070

Sousbie T., 2011, MNRAS, 414, 350

Sousbie T., Pichon C., Colombi S., Novikov D., Pogosyan D., 2008a, MNRAS, 383, 1655

Sousbie T., Pichon C., Courtois H., Colombi S., Novikov D., 2008b, ApJ, 672, L1

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

Sousbie T., Pichon C., Kawahara H., 2011, MNRAS, 414, 384

Spergel D. N. et al., 2003, ApJS, 148, 175

Teyssier R., 2002, A&A, 385, 337

Teyssier R. et al., 2009, A&A, 497, 335

Toomre A., Toomre J., 1972, ApJ, 178, 623

Welker C., Dubois Y., Pichon C., Devriendt J., Chisari E. N., 2018, AAP, 613, A4

Zel'Dovich Y. B., 1970, Astron. Astrophys., 5, 84

Zunckel C., Gott J. R., Lunnan R., 2011, MNRAS, 412, 1401

## APPENDIX A: EXTREMA CORRELATIONS

Peak theory originates from the so-called Kac–Rice formula (Kac 1943; Rice 1945) and was first used in a cosmological context in the pioneering work of Bardeen et al. (1986). Let us sketch here how the one- and two-point statistics of extrema can be obtained. First, for a random field $\rho$ (e.g. the cosmic density field), we define the moments

$$ \sigma_0^2 = \langle \rho^2 \rangle, \quad \sigma_1^2 = \langle (\nabla \rho)^2 \rangle, \quad \sigma_2^2 = \langle (\Delta \rho)^2 \rangle. \quad (A1) $$

Two characteristic lengths can be built from those moment $R_0 = \sigma_0/\sigma_1$ and $R_* = \sigma_1/\sigma_2$, as well as the spectral parameter

$$ \gamma = \frac{\sigma_1^2}{\sigma_0\sigma_2} = \frac{d \log \sigma(R)}{d \log R}. \quad (A2) $$

For simplicity, let us work with fields having unit variances:

$$ x = \frac{1}{\sigma_0} \rho, \quad x_i = \frac{1}{\sigma_1} \nabla_i \rho, \quad x_{ij} = \frac{1}{\sigma_2} \nabla_i \nabla_j \rho. \quad (A3) $$

The one-point PDF will be denoted $P(\mathbf{X})$ and the joint PDF $P(\mathbf{X}, \mathbf{Y})$, where $\mathbf{X} = \{x, x_i, x_i\}$ and $\mathbf{Y} = \{y, y_{ij}, y_i\}$ of dimension $(d+1)(d+2)/2$ represent the normalized field and its (up to second) derivatives, at two prescribed comoving locations $(\mathbf{r}_x$ and $\mathbf{r}_y$ separated by a distance $r = |\mathbf{r}_x - \mathbf{r}_y|$). In the supposingly Gaussian initial conditions, this joint PDF is a multivariate normal distribution

$$ \mathcal{N}(\mathbf{X}, \mathbf{Y}) = \frac{\exp \left[ -\frac{1}{2} \left( \frac{\mathbf{X}}{\mathbf{Y}} \right)^T \cdot \mathbf{C}^{-1} \cdot \left( \frac{\mathbf{X}}{\mathbf{Y}} \right) \right]}{\det|\mathbf{C}|^{1/2} (2\pi)^{(d+1)(d+2)/2}}, \quad (A4) $$

where $\mathbf{C}_0 \equiv \langle \mathbf{X} \cdot \mathbf{X}^T \rangle$ and $\mathbf{C}_\gamma \equiv \langle \mathbf{X} \cdot \mathbf{Y}^T \rangle$ are the diagonal and off-diagonal components of the covariance matrix

$$ \mathbf{C} = \begin{pmatrix} \mathbf{C}_0 & \mathbf{C}_\gamma \\ \mathbf{C}_\gamma^T & \mathbf{C}_0 \end{pmatrix}. \quad (A5) $$

All these quantities solely depend on the separation $r$ because of homogeneity and isotropy. Equation (A4) is sufficient to compute the expectation of any quantity involving the field, its first and second derivatives. In particular, the two-point correlation $\xi_{\text{crit}}(r, \nu)$ of critical points at threshold $\nu$ separated by $r$ is given by

$$ 1 + \xi_{\text{crit}}(r, \nu) = \frac{\langle \rho_{\text{crit}}(\mathbf{X}) \rho_{\text{crit}}(\mathbf{Y}) \rangle}{\langle \rho_{\text{crit}}(\mathbf{X}) \rangle^2}, \quad (A6) $$

with the ‘localized’ density of critical point

$$ \rho_{\text{crit}}(\mathbf{X}) = \frac{1}{R_*^d} |\det(x_{ij})| \delta_D(x_i) \delta_D(x - \nu). \quad (A7) $$

This density is formally zero unless the condition for a critical point is satisfied. In particular,

$$ \langle \rho_{\text{crit}}(\mathbf{X}) \rangle = \left( \frac{\sigma_2}{\sigma_1} \right)^d \int d\mathbf{X} \det(x_{ij}) \delta_D(x_i) \delta_D(x - \nu) P(\mathbf{X}) \\ \equiv \partial_\nu n_{\text{crit}}, $$

which appears in the denominator of equation (A6), is the average number density of critical points at threshold $\nu$, while

$$ \langle \rho_{\text{crit}}(\mathbf{X}) \rho_{\text{crit}}(\mathbf{Y}) \rangle = \int d\mathbf{X} \int d\mathbf{Y} P(\mathbf{X}, \mathbf{Y}) \rho_{\text{crit}}(\mathbf{X}) \rho_{\text{crit}}(\mathbf{Y}) $$
is the cross-correlation. If one wants to restrict the signature of the critical points (--- for peaks, -- + for filament-type saddles, - + + for wall-type saddles, and ++ + for minima), an additional constraint on the sign of the second derivatives is required. Then, the peak two-point correlation function reads

$$1 + \xi_{pp}(r, v) = \frac{\langle \rho_{pk}(\mathbf{X}) \rho_{pk}(\mathbf{Y}) \rangle}{\langle \rho_{pk}(\mathbf{X}) \rangle^2}. \quad (A8)$$

where the localized peak number density $\rho_{pk}(\mathbf{X})$,

$$\rho_{pk}(\mathbf{X}) = \frac{1}{R_{\star}^{d}} |\det(x_{ij})| \delta_{D}(x_{i}) \Theta_{H}(-\lambda_{i}) \delta_{D}(x - v), \quad (A9)$$

implements the peak condition. For $d > 1$, $\delta_d(x_i)$ is understood as $\prod_{i \le d} \delta_D(x_i)$, while $\Theta_H(-\lambda_i)$ stands for $\prod_{l \le d} \Theta_H(-\lambda_l)$, with $\{\lambda_l\}_l$ the eigenvalues of the Hessian. Because of these restrictions on the signature, the integral typically is not analytical. In dimension d, we define the conditional probability that $x_{ij}$ and $y_{ij}$ satisfy the PDF, subject to the condition that $x_i = y_i = 0$ and $x = y = v$, and rely on Monte Carlo methods in MATHEMATICA in order to evaluate numerically equation (A8). Namely, we draw random numbers of dimension $d(d+1)$ from the conditional probability that $x_{ij}$ and $y_{ij}$ satisfy the PDF, subject to the condition that $x_i = 0$ and $x = y = v$. For each draw $^{(k)}$ if $\lambda_l(x_{ij}^{(k)}) < 0$ and $\lambda_l(y_{ij}^{(k)}) < 0$ ($l \le d$), we keep the sample and evaluate $\det(x_{ij}^{(k)})\det(y_{ij}^{(k)})$ and otherwise we drop it; eventually,

$$\langle \rho_{pk}(\mathbf{X}) \rho_{pk}(\mathbf{Y}) \rangle \approx \frac{P(x=y=v, x_i=y_i=0)}{N} \sum_{k \in S} \det(x_{ij}^{(k)}) \det(y_{ij}^{(k)}),$$

where N is the total number of draws, and S is the subset of the indices of draws satisfying the constraints on the eigenvalues. The same procedure can be applied to evaluate the denominator $\langle \rho_{pk}(\mathbf{X}) \rangle = \bar{n}_{pk}(v)$. Equation (A8) then yields an estimation of $\xi_{pp}(r, v)$. This algorithm is embarrassingly parallel and can be easily generalized, for instance, to the computation of the correlation function $\xi_{pk}(r, >v)$ of peaks above a given threshold in density, the correlation between peaks and saddle points and to any dimension d. In practice, it is fairly efficient as the draw is customized to the shape of the underlying Gaussian PDF. Obviously, if correlation functions above a given threshold are considered, the required number of draws is larger and increases with the value of the threshold (as the event $x > v$ becomes rarer).

## APPENDIX B: FILAMENT CROSSINGS

We propose here to describe how to compute peak counts on the surface of a sphere around a central peak. The case of curved manifold for one-point statistics was addressed in Marcos-Caballero et al. (2016) and Pogosyan, Codis & Pichon (2016). Here, we develop a formalism to deal with two-point statistics when spherical sections of higher dimensional spaces are considered.

### B1 One-point statistics on the sphere

Following Codis et al. (2013), we use Cartesian coordinates with indices 1–3. Here, 3 will refer to the direction of the separation r between the central peak and the current point. The unit variance field under consideration is again denoted x and its unit variance first and second derivatives $x_i$ and $x_{ij}$ for i, j between 1 and 3. We want to find the function $B_{\text{fil}}(x(0), x_i(0), x_{ij}(0), x(r), x_i(r), x_{ij}(r))$ such that the mean number of 2D peaks at a distance r from a central peak reads $\langle B_{\text{fil}} \rangle$ where all variables follow a known Gaussian distribution. To

do so, we need to express the gradient, trace and determinant of the Hessian at the surface of the sphere in terms of the 3D field variables.

The surface gradient can be defined as

$$\nabla_{\perp} x = \nabla x - \frac{\partial x}{\partial r}. \quad (B1)$$

This operator acts over the unit sphere and is trivially independent of the frame orientation. In our frame, the surface gradient has coordinates $(x_1, x_2, 0)$.

The surface Laplacian (often called the Beltrami–Laplace operator) can also be defined in a frame-orientation independent way

$$r^2 \Delta_{\perp} x = r^2 \Delta x - r \frac{\partial^2(rx)}{\partial r^2}, \quad (B2)$$

which in our frame reads $I_1 \equiv \Delta_{\perp} x = x_{11} + x_{22} - 2x_3/r$.

The only remaining quantity to compute is now the determinant on the sphere. It can be shown that the Hessian matrix of covariant derivatives reads

$$H_{2D} = \begin{pmatrix} x_{11} - x_3/r & x_{12} \\ x_{12} & x_{22} - x_3/r \end{pmatrix}, \quad (B3)$$

so that the determinant simply reads $I_2 \equiv x_{11}x_{22} - x_{12}^2 - x_3/r(x_{11} + x_{22}) + x_3^2/r^2$. This surface determinant $I_2$ has therefore mean $\langle I_2 \rangle = \sigma_1^2/(3r^2)$ and the variance of the surface Laplacian $I_1$ has a curvature correction given by $\langle I_1^2 \rangle = 8/15\sigma_2^2 + 4/(3r^2)\sigma_1^2$.

In particular, the genus (i.e. signed critical points) can be computed $\chi = 4\pi r^2 \times \sigma_1^2/(3r^2) \times 3/(2\pi\sigma_1^2) = 2$ as expected. Note that the factor $3/(2\pi\sigma_1^2)$ comes from the zero gradient condition which reads $P(x_1 = 0, x_2 = 0) = (1/\sqrt{2\pi\sigma_1^2/3})^2$ as the variance of the gradient along each direction is $\sigma_1^2/3$.

We also recover that $J_2 = I_1^2 - 4I_2 = (x_{11} - x_{22})^2 + 4x_{12}^2$ is independent of the curvature.

### B2 Two-point statistics

We now consider the vector containing the field and its first and second derivatives at the origin and on a sphere at a distance r, $X = \{y, y_i, y_{ij}, x, x_i, x_{ij}\}$. The genus on the sphere given a peak constraint at the centre can then be computed as

$$\chi_{\text{pk}} = \frac{\langle I_2 \delta_D(x_i^{2D}) \times \det y_{ij} \delta_D(y_i) \delta_D(y - v) \mathcal{B}(\lambda_3 < 0) \rangle}{\langle \det y_{ij} \delta_D(y_i) \delta_D(y - v) \mathcal{B}(\lambda_3 < 0) \rangle},$$

where $\lambda_1 < \lambda_2 < \lambda_3$ are the eigenvalues of the Hessian matrix $y_{ij} \delta_D(y_i) = \delta_D(y_1)\delta_D(y_2)\delta_D(y_3)$ and $\delta_D(x_i^{2D}) = \delta_D(x_1)\delta_D(x_2)$, and is found to be exactly 2 as expected.

As a proxy for the number of skeleton branches around a peak, we can now compute the mean number of 2D peaks on the surface of the sphere around a central peak. For that purpose, a condition on the sign of the eigenvalues has to be implemented and leads to

$$N_{\text{fil}} = \frac{\langle I_2 \delta_D(x_i^{2D}) \mathcal{B}(\mu_2 < 0) \det y_{ij} \delta_D(y_i) \delta_D(y - v) \mathcal{B}(\lambda_3 < 0) \rangle}{\langle \det y_{ij} \delta_D(y_i) \delta_D(y - v) \mathcal{B}(\lambda_3 < 0) \rangle},$$

where $\mu_1 < \mu_2$ are the eigenvalues of the Hessian matrix $H_{2D}$. As displayed in Fig. 18, the number of peaks is two in the zero separation limit (while it would be one without zero gradient constraint at the centre of the sphere) as the local geometry is ellipsoidal and it does not change with the height of that central peak. As separation grows, bifurcation points occurs and the number of filaments increases.
## APPENDIX C: PERSISTENCE CUTS

To get the persistent skeleton of our two- and three-dimensional GRF realizations, we use the code DISPERSE. In addition to its ability to work with sampled data sets while assuming nothing about its geometry or homogeneity, this code can select structures on the basis of their significance via the notion of persistence ratio, a measure of the strength of the topological connection between individual critical points. This persistence threshold is expressed in terms of the typical rms of the noise. To set this threshold, in practice we use a persistence cut $p_{min}$ so that the number of peaks in the maps matches the expected number of peaks (which is known exactly for GRFs) with better than 0.5 per cent accuracy. This cut has to depend on the shape of the power spectrum. Indeed, if a constant persistence cut is chosen, the left-hand panel of Fig. C1 shows that the error on the total number of peaks in the map is significantly non-zero and is a strong function of the spectral index. Instead, we therefore select a cut different for each spectral index. In 2D, adopted persistence cuts are given in Table C1 and display on the right-hand panel of Fig. C1. The relation seems to be almost linear with the spectral index and well fitted by $p_{min} = 0.0031(1 + n_s/2.8)\sigma_0$.

The same procedure is followed in 3D. In that case, the persistence cuts are found to be $p_{min} = \{4.5, 6, 9, 12\}\sigma_0/1000$ for $n_s = \{-3, -2, -1, 0\}$.

For non-GRFs, the analytical prediction is not known in general. Instead we compare the measurements of DISPERSE (Sousbie, Pichon & Kawahara 2011) and subsequently adapt the persistence cut to measurements from the MAP2EXT code (Colombi et al. 2000; Pogosyan et al. 2011). Here, for every pixel a segment of quadratic surface is fit in the tangent plane based on the field values at the pixel of origin and its neighbours. The position of the extremum of this quadratic surface, its height and its Hessian are computed. The extremum is counted into the tally of the type determined by

<table><caption>**Table C1.** Cut in persistence $p_{min}$ (in units of $\sigma_0/1000$) for different values of the spectral index $n_s$ chosen so that the number of peaks is not different from the expected value by more than 10. This mapping is well fitted by the linear relation $p_{min} = 0.0031(1 + n_s/2.8)\sigma_0$.</caption><tr><td>$n_s$</td><td>1</td><td>0.5</td><td>0.</td><td>-0.5</td><td>-1</td><td>-1.5</td><td>-2</td></tr><tr><td>$p_{min}$</td><td>4.2</td><td>3.7</td><td>3</td><td>2.4</td><td>2.1</td><td>1.5</td><td>0.9</td></tr></table>

<figure><img src="image_28.png" /><figcaption>**Figure C1.** Left-hand panel: error on the number of 2D peaks measured compared to expected as a function of the spectral index when a fixed cut of $p_{min} = 0.003\sigma_0$ is chosen for the persistence. Right-hand panel: cut in persistence (in units of $\sigma_0/1000$) for different values of the spectral index $n_s$ chosen so that the number of peaks is not different from the expected value by more than 1 per cent. This mapping is well fitted by the linear relation $p_{min} = 0.0031(1 + n_s/2.8)\sigma_0$.</figcaption></figure>

its Hessian (two negative eigenvalues for peaks) if its position falls within the original pixel. Several additional checks are performed to preclude registering extrema in the neighbouring pixels and minimize missing extrema due to jumps in the fit parameters as region shifts to the next pixel. This procedure performs with better than 1 per cent accuracy when the map is smoothed with a Gaussian filter whose FWHM exceeds 6 pixels.

## APPENDIX D: GEOMETRY OF CONNECTION

Let us briefly investigate the geometry of the field in the vicinity of bifurcation points. As defined by Pogosyan et al. (2009) bifurcation points will locally correspond to two eigenvectors of the Hessian having the same eigenvalue.

### D1 Bifurcation in 2D

Let us expand the field, $x$, in the Hessian eigenframe to third order in the vicinity of $x_2 = 0$, $x_{11} = x_{22}$ along the direction $\delta\mathbf{r} = \epsilon(\cos\theta, \sin\theta)$; its variation, $\Delta x$ will obey (having set $\epsilon$ to one)

$$6\Delta x = 6x_1 \cos(\theta) + 3x_{22} + x_{111} \cos^3(\theta) + 3 \sin(\theta) x_{112} \cos^2(\theta) + 3 \sin^2(\theta) x_{122} \cos(\theta) + \sin^3(\theta) x_{222} + o(x_{ijk}).$$

Calling $y = \cos^2\theta$, the extrema of $\Delta x$ with respect to $\theta$ on the $\epsilon$-circle satisfy the algebraic equation

$$y \left( \frac{(y(3x_{112} - x_{222}) - 2x_{112} + x_{222})^2}{(2x_1 + y(x_{111} - 3x_{122}) + x_{122})^2} + 1 \right) = 1, \quad (\text{D1})$$

or geometrically

$$\left(\frac{1}{2} \delta \mathbf{r} \cdot \nabla \nabla \rho \cdot \delta \mathbf{r} + \nabla \rho\right) \times \delta \mathbf{r} = \mathbf{0}, \quad (\text{D2})$$

which is a cubic in $y$; hence at most six solutions are possible for $\theta$ which correspond to three maxima and three minima. Note that if the gradient dominates equation (D2) relatively to the third derivative tensor, the solutions $\pm\delta\mathbf{r}$ should be globally in the direction of $\nabla\rho$ (the dipole), whereas if it is subdominant (e.g. for the first bifurcation near the peak), $\delta\mathbf{r}$ should follow the 'eigendirections' of $\nabla\nabla\nabla\rho$ (the octupole), possibly weakly skewed by the local gradient. The inspection of equation (D2) shows indeed that the stronger the gradient $x_1$, the weaker the deviation of the bifurcation
from the gradient's direction ($y = 1$). If the octupole has special symmetries ($x_{111} = 3x_{122}$ and $x_{222} = 3x_{112}$) then only one solution in $y$ remains: $y_0 = (x_{112}^2/(2x_1 + x_{122})^2 + 1)^{-1}$. Note also that if $x_{222} = x_{112} = 0$ only $\cos \theta = \pm 1$ is a solution: this corresponds to a straight bifurcation since the corresponding curvature is also zero.

## D2 Bifurcation in 3D

Let us expand the field, $x$ to third order in the vicinity of $x_2 = x_3 = 0$, $x_{11} = x_{22}$ along the infinitesimal direction $\delta\mathbf{r} = \epsilon(\cos \theta \cos \phi, \cos \theta \sin \phi, \sin \theta)$. The extrema of $\Delta x$ on the infinitesimal sphere of radius $\epsilon$ obey

$$ \frac{\partial \Delta x}{\partial \delta \mathbf{r}} \times \delta \mathbf{r} = \left( \frac{1}{2} \delta \mathbf{r} \cdot \nabla \nabla \nabla \rho \cdot \delta \mathbf{r} + \nabla \nabla \rho \cdot \delta \mathbf{r} + \nabla \rho \right) \times \delta \mathbf{r} = \mathbf{0}. $$

Again the critical directions correspond to a mixture between the eigenvector of the anisotropic part of the Hessian, the gradient, and the 'eigenvectors' of $\nabla \nabla \nabla \rho$, depending on the relative strength of the three components. Unless the anisotropic part of the Hessian is null the bifurcation will be coplanar, and can be investigated with the 2D formalism presented above.

## APPENDIX E: ND CONNECTIVITY

In order to compute the connectivity in arbitrary dimensions, it is best to move to the eigenframe of the Hessian and rely on the joint statistics of its eigenvalues. From Pogosyan et al. (2009), the probability of measuring the set of $d$ (ordered) eigenvalues of the $d$-dimensional Hessian $\{\lambda_i\}$ and density $\nu$ obeys

$$ \prod_{i \le d} d \lambda_i \prod_{i < j} (\lambda_j - \lambda_i) \exp \left( -\frac{1}{2} Q_{\nu}(\nu, \{\lambda_i\}) \right), \quad (E1) $$

where $Q_{\nu}$ is a quadratic form in $\lambda_i$ and $\nu$ given by

$$ Q_{\nu}(\nu, \{\lambda_i\}) = \nu^2 + \frac{(\sum_i \lambda_i + \gamma \nu)^2}{(1 - \gamma^2)} + Q_d(\{\lambda_i\}), \quad (E2) $$

with

$$ Q_d(\{\lambda_i\}) = (d + 2) \left[ \frac{1}{2}(d-1) \sum_i \lambda_i^2 - \sum_{i \neq j} \lambda_i \lambda_j \right]. \quad (E3) $$

It now follows that the extrema number counts in dimension $d$ read:

$$ \partial_{\nu} n^d = \langle \Theta_H(-\lambda_i) | \prod_i \lambda_i | \rangle, \quad (E4) $$

where this expectation is computed using equation (E5). From equation (E2), integration over $\nu$ yields the marginal probability of $\{\lambda_i\}$:

$$ \prod_{i \le d} d \lambda_i \prod_{i < j} (\lambda_j - \lambda_i) \exp \left( -\frac{1}{2} Q_d(\{\lambda_i\}) - \frac{1}{2} \left( \sum_i \lambda_i \right)^2 \right). \quad (E5) $$

Finally, the $d$ dimensional connectivity is simply given by the ratio

$$ \kappa = \frac{2 \langle \Theta_H(-\{\lambda_i\}_{i \le d}) \Theta_H(\lambda_d) | \prod_i \lambda_i | \rangle}{\langle \Theta_H(-\{\lambda_i\}_{i \le d}) | \prod_i \lambda_i | \rangle}, \quad (E6) $$

where this expectation is now computed using equation (E5). Implementing equation (E6) in dimension up to 10 yields

Table E1 which gives a few values of $\kappa$ as a function of dimension. The global connectivity of a GRF is shown to slightly depart from a cubic lattice in dimensions greater than 3. These defects are shown to affect the connectivity even more so as the dimension increases with a scaling proportional to the power 7/4 of the dimension, so that a good fit to the global connectivity is given by

$$ \kappa = 2d + \left( \frac{2d - 4}{7} \right)^{7/4}. \quad (E7) $$

We checked that integration over the field variables, $\mathbf{X} = (x, x_{ij}, x_i)$ also yield the same numbers.

<table><caption>**Table E1.** $\kappa_G$ as a function of dimension $d$ for GRFs. The scaling is well fitted by $\kappa = 2d + ((2d-4)/7)^{7/4}$.</caption><thead><tr><th>d</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th></tr></thead><tbody><tr><th>κ</th><td>4</td><td>6.11</td><td>8.35</td><td>10.73</td><td>13.23</td><td>15.85</td><td>18.7</td><td>21.4</td><td>24.4</td><td>27.4</td></tr></tbody></table>

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