{'fraction_non_alphanumeric': 0.033593976252534026, 'fraction_numerical': 0.00974997586639637, 'mean_word_length': 4.5309663641217295, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Skeleton: Connecting Large Scale Structures to Galaxy Formation

We report on two quantitative, morphological estimators of the filamentary structure of the Cosmic Web, the so-called global and local skeletons. The first, based on a global study of the matter density gradient flow, allows us to study the connectivity between a density peak and its surroundings, with direct relevance to the anisotropic accretion via cold flows on galactic halos. From the second, based on a local constraint equation involving the derivatives of the field, we can derive predictions for powerful statistics, such as the differential length and the relative saddle to extrema counts of the Cosmic web as a function of density threshold (with application to percolation of structures and connectivity), as well as a theoretical framework to study their cosmic evolution through the onset of gravity-induced non-linearities.

Over the course of the last decades, our understanding of the extragalactic universe has undergone a paradigm shift: the description of its structures has evolved from being (mostly) isolated to being multiply connected both on large scales, cluster scales and galactic scales. This interplay between large and small scales is driven in part by the scale invariance of gravity which tends to couple dynamically different processes, but also by a the strong theoretical prejudice associated with the so-called concordant cosmological model [1]. This model predicts a certain shape for the initial conditions, leading to a hierarchical formation scenario, which predicts the formation of the so-called Cosmic Web, the most striking feature of matter distribution on megaparsecs scales in the Universe. This distribution was confirmed, observationally, more than twenty years ago by the first CfA catalog [2], and later by the SDSS [3] and 2dFGRS [4] catalogs. On these scales, the "Cosmic Web" picture [5] 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 [6]. The analysis of the connectivity of this filamentary structure is critical to map the very large scale distribution of our universe to establish, in particular, the percolation properties of the Web [8].
On intermediate scales, the paradigm shift is sustained by panchromatic observations of the environment of galaxies which illustrate sometimes spectacular merging processes, following the pioneer work of e.g. [9] (motivated by theoretical investigations such as [10]). The importance of anisotropic accretion on cluster and dark matter halo scales [11,12,13] is now believed to play a crucial role in regulating the shape and spectroscopic properties of galaxies. Indeed it has been claimed (see e.g. [14,15]) that the geometry of the cosmic inflow on a galaxy (its mass, temperature and entropy distribution, the connectivity of the local filaments network, etc.) is strongly correlated to its history and nature. Specifically, simulations suggest that cold streams violently feed high redshift young galaxies with a vast amount of fresh gas, resulting in very efficient star formation. One of the puzzles of galaxy formation involves understanding how galactic disks reform after minor and intermediate mergers, a process which is undoubtedly controlled by anisotropic gas inflow.
Recently, [16] presented a method to compute the full hierarchy of the critical subsets of a given density field. This approach connects the study of the filamentary structure to the geometrical and topological aspects of the theory of gradient flows [17]. In this paper, we focus on the connectivity of the corresponding network. Specifically, since galaxy formation seems to be geometrically regulated by the accretion of cold gas from the LSS, how can we use the skeleton to characterize this anisotropic accretion, and predict its evolution through perturbation theory ?
Let us first qualitatively introduce the two operating skeleton extraction algorithms, and present our findings

The skeleton: algorithms
We may qualitatively define the skeleton as the 3D analog of ridges in a mountainous landscape. It consists of the special lines which connect saddles and peaks together. Mathematically these are critical solutions of the gradient flow (ṙ = ∇ρ) between critical points. Here we want to construct such 3D ridges to trace the filaments of the large scale structures. Two venues have been explored recently: the so called local skeleton [19,20], which defines a (degenerate) point-like process corresponding to the zeros of a set of functions; the idea is that on a ridge the gradient should be along the direction of least curvature (see figure 1, left panel), which translates into the set of 3 equations: where ∇ρ and ∇∇ρ are respectively the gradient and the Hessian of the field, ρ and λ 1 > λ 2 > λ 3 its eigenvalues, with the condition λ 1 + λ 2 ≤ 0, which ensures that the critical line is a ridge. Alternative conditions on the eigenvalues of the Hessian may be applied to pick up the whole set of critical lines, see figure 3. Following [19], this local skeleton and its properties provide an alternative description to classical approaches of galaxy clustering (powerspectrum, bispectrum etc.) and attempt to achieve data compression via a mathematical description of the morphology of the cosmic web. An alternative, global definition [19,16] is to define it as the border between N void-patches, where N is the dimensionality and void-patches are the attraction basins of the (anti-)gradient of the field (see figure 1, right panel in 2D, and figure 2 in 3D). The actual implemented algorithm is based on a watershed technique and uses a probability propagation scheme to improve the quality of the segmentation (see also [18]). It can be applied within spaces of arbitrary dimensions and geometry. The corresponding recursive segmentation yields the network of the primary critical lines of the field: the fully connected skeleton that continuously link maxima and saddle-points of a scalar field together. Both constructions are of interest: the local formulation provides means of conducting detailed statistical investigation of the corresponding degenerate point-like process, while the global skeleton yields a totally connected network of lines. This critical set of lines is a compact description of the geometry of the field, richer than the knowledge of the critical points alone. For the purpose of this review, it also allows us to explore the connectivity of peaks, both from the point of view of the number of connections to other peaks, and of the number of skeleton branches leaving a given maximum. Conversely, the local skeleton formulation allows us to extend the BBKS theory [21] of peaks in the context of critical lines to investigate their cosmic evolution, see below. Figure 3 compares the two sets of lines in the neighbourhood of a given peak.

The skeleton: connectivity
Let us now make use of this algorithm to explore the connectivity of the corresponding network. A set of two dimensionnal Gaussian random fields (GRF) is produced, and for each of them, the peak patch and the skeleton was generated and its connectivity computed, following the second prescription described in [16]: we chose here to smooth and label the skeleton having fixed the field extrema, (i.e. not the bifurcation points of the skeleton, see Figure 4 (top left panel)). By fixing the extrema of the field, one ensures that the skeleton subsets that link these extrema are treated independently.
Let us first focus on the degree of the peakpatch, defined as the number of saddle points within a given patch which connects the skeleton of one patch to its neighbours. Hence the connectivity count reflects the number of maxima a given peak is connected to (or in the language of graph theory, the most likely degree of the vertices:n = n ). Figure 4 (top right panel) displays this PDF as a function of n; in particular it allow us to compute n , the mean n which is found to be independent of γ, the shape parameter of the underlying GRF. The bottom left panel shows the corresponding normalized PDF(n, η) versus the contrast η. As expected the number of connections increases with contrast. Indeed, near the maximum at high contrast all eigen values tend to become equal [22]. Therefore all incoming directions become possible.
Let us now take the other (astrophysically relevant) perspective and count the actual number of filaments incident onto a given maximum. This local (intra patch) degree is in fact equal to the number of saddle points minus the number of bifurcation within the peakpatch: n max = n saddle − n bifurcation . The bottom right panel of figure 4 shows the corresponding PDF(n max ). Note that this distribution is almost symmetrical, centered at n max = 3. Figure 5 (left panel) presents the same distribution in 3D, which is very skewed and presents a sharp mode at 3. Its cosmic evolution, measured in CDM dark matter simulations with and without dark energy is qualitatively shown on the right panel. Our purpose in the rest of this paper is to explain qualitatively this cosmic trend by 1) deriving the statistics of bifurcation points within each patch and 2) predicting the cosmic evolution of saddle points and peaks.  (1). Primary lines are shown in solid and secondary lines are dashed. The green lines correspond to global critical lines: the skeleton and the anti-skeleton, which delineate a special bundle of gradient lines at resp. the intersection of a peak-patch and a void-patch [17]. The primary local lines follow fairly well the gradient lines, noticeably near the saddle points, where the "stiff" approximation holds best. In contrast, the approximation worsens for the secondary critical lines. The main distinction between the global and local skeletons is that the global one follows everywhere the smooth gradient line that uniquely connects a maximum to a saddle, at the cost of deviating from being exactly on the ridge (see how in the vicinity of the minimum at the bottom, the right green line does not follow the valley). The local skeleton tries to delineate the ridges as far from extrema as possible, but then the lines that follow this local procedure from different extrema do not meet and have to rather suddenly reconnect. A particularly striking example is the loop on the right hand side. Bottom left medallion: the neighbourhood of a local critical line (thick blue line). Thin lines are isocontours of the field. Three sample points are investigated in detail. The signature, orientation and the magnitude of the local Hessian are represented by the golden shapes. Near the maximum on the right edge, the signature of the eigenvalues of the Hessian is (-1,-1), which is shown by ellipses oriented according to eigen-directions with longer semi-axis along the direction of the least curvature. At the leftmost point the eigenvalue signature is "saddle-like", (1,-1), which is represented by a pair of hyperbolae, also oriented with respect to eigen-directions. By definition, on the critical line the gradient of the field ∇ρ, shown by red arrows, is aligned with one of the eigen-directions (i.e the axis of the ellipse or hyperbola in the graph). The light cyan arrows are the tangent vectors to the critical line ∝ ε · ∇S, while stiff approximation to them would be parallel to the gradient. The direction of the critical line is close to the gradient when it follows the ridge near the maximum, but slides at an angle in the "saddle-like" region, before joining the corresponding saddle extremal point. Note that the gradient line that takes us to the same saddle as a segment of the global skeleton (dashed green) does not follow the ridge too closely in this instance.

The skeleton: statistics
As illustrated first in 3D on Figure 7, the critical condition given by equation (1) defines 3 isosurfaces corresponding to its components, and the local skeleton direction is simply given by the cross products of two of its normals. Hence Top left panel: A typical (gray coded) 2D peak-patch of a Gaussian random field together with the degree of each vertex of the skeleton graph; for instance, within the red star the mean degree is 4, while the corrected degree is 3 within the blue star, as one bifurcation point is present; top right panel: the corresponding PDF of this degree, for different values of the shape parameter γ as labeled; the mean degree is 4 and corresponds to the theoretical expectation of 2n saddle /n max ; bottom left panel: the (log) differential PDF of the degree as a function of contrast: as expected, the denser peaks are more connected; interestingly, the distribution seems Gaussian at fixed n even though its marginal isn't; bottom right panel: the corrected PDF of the number of branches connecting onto a given peak, which involves substracting the number of bifurcation points within the patch to the number of saddle points on the edge of the patch. the differential length (per unit volume) is simply given by the statistical expectation where ergodicity allowed us to replace volume average by ensemble average over the statistical distribution, P(η, x k , x kl , x klm ), of the successive (a-dimensional) derivatives, σ 1 x k ≡ ∇ k ρ, σ 2 x kl ≡ ∇ k ∇ l ρ, σ 3 x klm ≡ ∇ m ∇ l ∇ k ρ of the field, ρ. The scale factor 1/R * 2 = σ 2 2 /σ 2 1 in equation (2) follows dimensionally (given σ 2 2 σ 1 ∇s i ≡ ∇S i , σ 0 η ≡ ρ). Here the variances, σ 2 p , obey σ 2 p = 2π N/2 /Γ[N/2] ∞ 0 k 2p P(k)k N−1 dk in N dimensions. (see [20] for a formal derivation of equation (2)). Left panel: the 3D corrected PDF of the number of branches connecting onto a given peak for a GRF: the mode is at 3; insert: the typical neighbourhood entropy map of a massive galaxy in the Marenostrum simulation from [15]; right panel: the qualitative cosmic evolution of the mean degree in dark matter simulation with and without dark energy; as expected this degree decreases with the expansion factor as connectors get washed out by gravitationnal clustering; the rate of decrease seems sensitive to the moment when dark energy kicks in for ΛCDM models (around a ∼ 0.5). Note that the absolute amplitude of this mean degree is not shown here as we expect the global skeleton to over estimate the mean connectivity in 3D (see main text).
More generally, for ND critical lines, the N-1 independent functions S i that define the critical condition (1) acquire the following reduced form in the eigenframe of the Hessian of the field: s i = x 1 x i (λ 1 − λ i ) = 0. The measurements in [23] found that, over the range of spectral indexes relevant to cosmology, the third derivatives of the field x klm play a subdominant nature. In the so called "stiff" approximation we therefore omit the third derivative, effectively assuming that the Hessian can be treated as constant during the evaluation of ∇s. This picture corresponds to a skeleton connecting extrema with relatively straight segments. Let us first assume that the underlying field is Gaussian. In the stiff approximation, the gradients s i k ≡ ∇s i have just two non-zero components, s i 1 = x i λ 1 (λ 1 − λ i ) (which vanishes on the critical line) and s i i = x 1 λ i (λ 1 − λ i ), which shows that in this approximation we equate the direction of the line with the gradient of the field. Substituting this expression into equation (2) and integrating over δ D (s i ) = δ D (x i )/(x 1 (λ 1 − λ i )) we obtain a simple expression for the differential length of the ND-critical lines: where the shape parameter γ = σ 2 1 /(σ 0 σ 2 ) describes the correlation between the field and its second derivatives, Q γ is a quadratic form in λ i and η which functional form is Equation (3) simply states that the stiff differential length is the expectation of the product of the "smaller" eigenvalues, which is quite reminiscent of the classical extremal count result (the expectation of the product of all eigenvalues). It also qualitatively makes sense, as the larger the curvarture orthogonal to the skeleton, the more skeleton segment one may pack per unit volume. Since the argument of Q γ is extremal as a function of η when γη ∼ ∑ i λ i , the largest contribution at large γη in the integral should arise when λ i ∝ γη since near the maximum at high contrast all eigen values are equal [22]. Hence given that ∏ i< j (λ i − λ j ) is the measure, the only remaining contribution in the integrand comes from |∏ i>1 λ i | ∝ (λ η) N−1 , and the dominant term at large η is given by Note that in 2D the differential length, equation (3) can be reduced to a particularly simple closed form and for the integrated skeleton length In other words, one expects to find one segment of skeleton per linear section of ≈ (4.2R * ). Similarly, in 3D one segment of skeleton is found per surface section of ≈ (4.65R 2 * ). The match of the detailed PDF with the corresponding measurements are shown in Figure 6 in 2 and 3D. From the point of view of cosmology, the skeleton is invariant w.r.t. any monotonic bias, and traces the denser regions of the field. As shown in equations (5) and (7) its statistical description yields a measure of both R * and R 0 , hence on the shape of the underlying powerspectrum on the corresponding smoothing scale over which the skeleton was computed. An implementation of this estimate on the SDSS catalog was carried by [24] and provided constraints on Ω M .

Bifurcation counts
In [25] we conjectured that critical lines experience a qualitative change in behaviour in the vicinity of the points where either a Hessian eigenvalue orthogonal to the gradient direction vanish, or becomes equal to the one along the gradient. The first type corresponds to points where the curvature transverse to the direction of the critical line vanishes along at least one axis: typically, in 2D, they mark regions where a crest becomes a trough, or vanish into a plateau. The second type correspond to points where the critical lines would split, even though the field does not go through an extremum: a bifurcation of the lines occurs along the slope; the occasional skier or mountaineer will be familiar with a crest line splitting in two, even though the gradient of the field has not vanished. Namely, if, for definiteness, ∇ρ is taken to be along the first eigen-direction, λ 2 = 0, or λ 2 = λ 1 . We called the first case the "sloping plateau" as it designates the entering of a flat region, and the second, the "bifurcation" as it designates the places of possible reconnection of critical lines. Remarkably, these special points on the critical lines are recovered by the To be more specific, in 2D, the skeleton's singular points correspond to points where S k ≡ ∇ k S = 0. The number density, n B (η) of singular points below the threshold η is equal to n B (η) = η>x dxd 2 x k d 3 x kl d 4 x klm d 5 x klmn P(x, x k , x kl , · · · )|det (∇ k ∇ l s) |δ D (s 1 )δ D (s 2 ) .
The gradient of S, evaluated in the stiff approximation, in the Hessian eigenframe has the components s stiff 1 = x 2 λ 1 (λ 1 − λ 2 ) , and s stiff 2 = x 1 λ 2 (λ 1 − λ 2 ) , and involves only second derivatives of the field. Let us consider the critical line that corresponds to the x 2 = 0 condition in the Hessian eigenframe. Then s stiff 1 vanishes everywhere along this line. The requirement s stiff 2 = 0 has a solution at the extremal points, x 1 = 0, but also in two other cases, namely λ 2 = 0 or λ 2 = λ 1 , that we conjectured to be of interest. The second situation (isotropic Hessian) has a number density given by We note that the number density of "bifurcation" points is proportional just to the PDF of the field and, consequently the bifurcation points are as frequent in the regions of high field values as in the low ones. Finally, in 3D, we expect the finite resolution bifurcation branch (connecting maxima to bifurcation points) to become bifurcation plates; in turn the boundary of these plates will be counted as critical lines by our voidpatch algorithm; we therefore expect that the algorithm numerically over estimate the degree of peakpatches. Note nonetheless that the 3D counterpart of equation (9) should allow us to correct for the number of bifurcations within each patch and compute the statistics of the number of branches connected to a given maximum.

Departure from Gaussianity
While the Gaussian limit provides the fundamental starting point in the study of random fields [26,27,21], non-Gaussian features of the CMB and the large scale structure (LSS) fields are also of great interest. CMB inherits high