The local theory of the cosmic skeleton

The local theory of the critical lines of 2D and 3D Gaussian fields that underline the cosmic structures is presented. In the context of cosmological matter distribution the subset of critical lines of the 3D density field serves to delineate the skeleton of the observed filamentary structure at large scales. A stiff approximation used to quantitatively describe the filamentary skeleton shows that the flux of the skeleton lines is related to the average Gaussian curvature of the 1D (2D) sections of the field, much in the same way as the density of the peaks. The distribution of the length of the critical lines with threshold is analyzed in detail, while the extended descriptors of the skeleton - its curvature and its singular points, are introduced and briefly described. Theoretical predictions are compared to measurements of the skeleton in realizations of Gaussian random fields in 2D and 3D. It is found that the stiff approximation predicts accurately the shape of the differential length, allows for analytical insight, and explicit closed form solutions. Finally, it provides a simple classification of the singular points of the critical lines: i) critical points; ii) bifurcation points; iii) slopping plateaux.


INTRODUCTION
The concept of random fields is central to cosmology. Random fields both provide initial conditions for the evolution of the matter distribution in the Universe, and represent how the observed signals manifest themselves in 3D, (e.g., in the galaxy or matter density inhomogeneities that form the Large Scale Structure (LSS)), or on the 2D sky (e.g. for the Cosmic Microwave Background (CMB) temperature and polarization, the convergence or shear in weak lensing maps). In the modern cosmological theories where initial seeds for inhomogeneities observed as cosmic structures have quantum origin, the fields of initial density fluctuations (and velocities) are Gaussian. Subsequent evolution retains Gaussianity for the observables that evolve linearly (CMB, very Large Scale Structure) while developing non-Gaussian signature if non-linear effects are involved (e.g. lensing and LSS at smaller scales).
While comparing the observational data to cosmological theory, in particular in order to estimate parameters of cosmological models, the emphasis is traditionally placed on the statistical descriptors of the random fields. For Gaussian fields the two-point correlation function or the power spectrum provide full statistical information, while non-Gaussian properties may be reflected in multi-point correlations. The understanding and the description of the morphology of structures in our Universe, on the other hand, calls for the studies of the geometry and topology of random fields. This subject has an extensive history from the early description of the one-dimensional radio signal time-streams in 1940's, to the study of the 2D ocean wave patterns in 1960's (Longuet-Higgins 1957) to 3D dimensional fields (Adler 1981) that found the most fruitful application in cosmology (Arnol'd et al. 1981; Bardeen et al. 1986). The most prominent geometrical objects in a typical realization of a random field are rear events -regions of unusually high or low values of the field. The rare events are usually related to the most spectacular observed objects -clusters of galaxies at low z, large protogalaxies at high-z or extensive voids. They are associated with the neighbourhoods of extrema -maxima or minima -making studies of such critical points the first step in understanding typical geometry of a field (Kaiser 1984; Bardeen et al. 1986;Regos & Szalay 1995;Scannapieco et al. 2006). The behaviour of the field in the neighbourhood of a rare peak is highly correlated with the peak properties, which allows to describe not only extrema but the extended peak-patch region (Bond & Myers 1996a) as a point process that involves the field and its successive derivatives. Including the shear flow into consideration gives a compelling application of the geometry of rare events to the description of cluster formation through the peak-patch collapse (Bond & Myers 1996b).
The rare events reflect the organization of the field around them and by and large determine the way the high (low) field regions are interconnected by the bridges of enhanced field values. In application to cosmology, the "Cosmic Web" picture emerges, which 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 supercluster scales (Bond et al. 1996). The study of the connectivity of filamentary structures reveal the role of the remaining type of critical points, the saddle extrema, in establishing, in particular, the percolation properties of the Web (Colombi et al. 2000). The next step naturally involves describing the statistical properties of these filamentary structures (Pogosyan et al. 1998;Schmalzing et al. 1999) and developing techniques for mapping the filaments in the simulation and data. Novikov et al. (2006) presented a 2D algorithm to trace the filaments of a density field while introducing the skeleton as the set of locally defined critical lines emanating from the critical points.  (hereafter SPCNP) extended the local theory and algorithm to three dimensions and provided the foundation for this work while introducing the "stiff" approximation. Recently,  presented an algorithm to map out a fully connected version of the skeleton that is defined according to the global properties as the lines of intersections of the patches (see also Aragón-Calvo et al. (2007), Platen et al. (2007) for alternative algorithms). This approach connects the study of the filamentary structure to the geometrical and topological aspects of the theory of gradient flows (Jost 2008) and returns the focus to the notions of peak and void patches.
This paper presents a consistent local theory for the cosmic skeleton, while focusing on the stiff approximation to compute the differential length of the skeleton as a function of the contrast and modulus of the gradient of the field. It allows us to define precisely how the properties of the skeleton depend analytically on the underlying spectral parameters, and understand what type of line prevails where. The crucial advantage of the local approach to the critical lines is that it allows to cast the statistical treatment of the linear objects as a point process that involves the field and its derivatives, which allows for analytical insight, and explicit closed form solutions. Our purpose is to construct the theory of critical lines of a given field corresponding to an intermediate representation of the field, which is more extended than the knowledge of the critical points.
The organization of the paper is the following. Section 2 classifies the various critical lines in 2D and 3D, connects the average length in a unit volume to the flux of the skeleton lines and, within the stiff approximation, to the average Gaussian curvature of the field in transverse sections. It also discusses the meaning of this approximation. Section 3 calculates the differential length of all sets of critical lines in 2D, while Section 4 investigates the corresponding 3D set of critical lines. More generally, the expression for the differential length of the N dimensional skeleton is sketched in Appendix A. Section 5 introduces the extended descriptors of the skeleton, Section 5.3 describes their singular points, while Section 6 provides the discussion and the summary. Appendix D gives the general method for obtaining in close form the joint distribution of the field and any combination of its derivative tensors in arbitrary dimensions. In particular, it exhibits all the statistical invariants and their dependence on the spectral parameters.

Local definition and classification
The subject of our investigation is a random field, ρ(r), that in a cosmological setting describes, for example, the density of the matter in the Universe, or the projected distribution of Cosmic Microwave light on the celestial sphere. Our focus is on the geometrical properties of the critical lines, that connect extrema of the field mapping out the filamentary ridges and valleys of the field. SPCNP have introduced the definition of the local critical lines as the set of points where the gradient of the density, ∇ρ, is an eigenvector of its Hessian matrix, H, H · ∇ρ = λ∇ρ i.e, the gradient and one of the principal curvature axes are collinear. Formally, this can be specified by a set of equations S ≡ (∇ρ · H) · · (∇ρ) = 0 , where is the fully antisymmetric (Levi-Civita) tensor of rank N . In general S is an antisymmetric N − 2 tensor. In 3D, the function S is vector-valued, S i = P klm ikl (∇mρ)H m k (∇ l ρ). However, zeroes of S determine a set of lines rather than isolated points. Let us consider the behaviour of S function in the neighbourhood of a point r = 0 that satisfies criticality condition S(0) = 0: In our case, under the condition S i = 0, the matrix ∇ k S i by definition possesses the left null-vector, furnished by the density gradient, P i ∇iρ`∇ k S i´= 0; hence the gradients ∇S i are not linearly independent. Consequently, there is a nontrivial solution for the right null-vector P k`∇ k S i´δ r k = 0, which determines the local direction of the line along which the criticality condition is maintained, S(δr) = 0. The critical lines intersect where ∇ k S i admits more than one independent right null-vector.
(i) The skeleton, that has the gradient in the λ1 direction and is limited to the region |λ1| |λ2|, which translates to the condition λ1 + λ2 0. The skeleton has always eigenvalues in the directions transverse to ∇ρ negative, λ3 λ2 0 and corresponds to the filamentary ridges spreading from the maxima in the direction of the slowest descent.
(ii) The anti-skeleton, that has the gradient in the λ3 direction and is restricted to the region |λ3| |λ2|, i.e λ3 + λ2 0. In the directions transverse to ∇ρ the anti-skeleton has always positive curvature λ1 λ2 0. It corresponds to the filamentary valleys spreading from the minima in the direction of the slowest ascent. Anti-skeleton can be viewed as a skeleton of the −ρ field.
The formal classification of the critical lines is summarized in Table 1.

The average flux (length per unit volume) of the critical lines
As the average number density is the fundamental quantity that describes point events, e.g. extrema of a field, conversely the most important characterization of the critical lines or skeleton is their flux, i.e. the number of critical lines intersecting a given oriented surface 1 . This flux is equivalent to the length of the lines per unit volume. Following NCD and SPCNP we shall preferentially use the latter terminology as it highlights that we deal with the first geometrical parameter, the length, of the lines. The subsequent parameters of these linear objects are the curvature, and, in 3D, the torsion.
In this paper we consider ρ to be a homogeneous and isotropic Gaussian random field of zero mean, described by the power spectrum P (k). In the statistical description of the skeleton of the field ρ, several linear scales are involved where σ 2 0 = ρ 2 , σ 2 1 = (∇ρ) 2 , σ 2 2 = (∆ρ) 2 , σ 2 3 = (∇∆ρ) 2 , and generally σ 2 (4) These scales are ordered R0 R * R . . .. The first two have well-known meanings of typical separation between zerocrossing of the field R0 and mean distance between extrema, R * ( Bardeen & al. 1986), and the third one,R is, by analogy, the typical distance between the inflection points. These three are the only ones that are involved in determination of the length of the critical lines. The higher order scaleR appear only in computation of the curvature and the torsion (see Section 5).
Let us define a set of spectral parameters that depend on the shape of the underlying power spectrum. Out of these five scales four dimensionless ratios may be constructed that are intrinsic parameters of the theory , and generally γp,q = σ 2 (p+q)/2 σpσq .
From the geometrical point of view γ specifies how frequently one encounters a maximum between two zero crossings of the field, whileγ describes, on average, how many inflection points are between two extrema. From a statistical perspective, γ's are the cross-correlation coefficients between the field and its derivatives at the same point (see Appendix D).
For Gaussian fields, these parameters can be easily calculated from the power spectrum. All γ's range from zero to one. For reference, for the power-law spectra with index n > −3, smoothed at small scales with a Gaussian window, γ = p (n + 3)/(n + 5), γ = p (n + 5)/(n + 7). Note that cosmologically relevant density power spectra have n > −3 and, thus, while γ can attain low values,γ are always close to unity 2 .
Let us introduce the dimensionless quantities for the field and its derivatives as well as for the functions S i and their gradients ∇S i : Note the specific choice of scaling for ∇S which is convenient in view of the subsequent development of the so-called "stiff" approximation. SPCNP has shown that in terms of these dimensionless quantities, the cumulative length per unit volume of the total set of critical lines below the threshold η is given by where a pair ∇s i , ∇s j can be chosen arbitrarily as long as it is linearly independent. In this equation |∇s i × ∇s j | reflects the inverse characteristic area orthogonal to a critical line per one such line while the two δD-functions account for the critical line condition (1). For the complete set of critical lines, there are no restriction to the region of integration. If one is interested in a particular type of the critical lines, the integration should be restricted to the regions consistent with Table 1. The differential length (per unit volume) is simply given by the derivative of equation (9) with respect to η: while the total length of critical lines is Since for Gaussian field, the derivatives of even order are uncorrelated with the odd orders, the joint distribution function P (x, x k , x kl , x klm ) entering equation (9) is factorized as In P0, the only dependence on the power spectrum of the field is through the parameter γ (c.f. equation (5)) that describes the correlation between the field and its second derivatives. Similarly P1(x k , x klm ) only involvesγ which describes the correlation between the gradient of the field and its third derivatives. Therefore, ∂L/∂η depends only on η,R γ andγ. The integrated length, L may depend only onγ andR since the marginalization of P0(η, x kl ) over η eliminates the dependency over γ.

The "stiff " filament approximation
Let us look at the dependence of the length of the critical lines on characteristic scales of the field in more detail. The R −2 * factor that appeared in equation (10) reflect our choice of dimensionless variables (8) and is suggestive but not yet conclusive since |∇s i × ∇s j | that includes third derivative terms, depends also on the other scale,R. Let us write formally If the third derivatives are important and the first term dominates, then the length scaling L ∝γ −2 R −2 * =R −2 would reflect the mean separation between the inflection points,R. Indeed, by definition the local skeleton is almost straight within a volume that has one inflection point ∼R 3 . A straight segment through such volume has length ∼R, thus the expected length per unit volume is ∼ 1/R 2 . But if the last term dominates statistically, the length per unit volume of the skeleton will scale as L ∝ R −2 * that can be interpreted that the critical lines are almost straight within a large volume volume ∼ R 3 * containing one extremum. This is consistent with observation that the integral term does not depend on the third derivatives, thus inflection points play no role, and any dependence onγ drops out.
Which regime holds can be established by measuring the dependence of the critical lines length in the simulations as a function of smoothing length for different spectral indexes. For the power-law spectra with Gaussian smoothing at the radial scale σ, in 3D, R * = √ 2σ/ √ n + 5, whileR = √ 2σ/ √ n + 7. The measurements in SPCNP found that L ∝ (n + 5.5)σ −2 over the range of spectral indexes relevant to cosmology, which points at the subdominant nature played by the third derivatives. In the "stiff" approximation we omit the third derivative, effectively assuming that the Hessian can be treated as constant Figure 1. The neighbourhood of a local critical line (thick blue line). This is a zoomed section of the wide patch shown in Figure 2. 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 saddle extremal point beyond the left edge of the plot (see Figure 2). 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.
during the evaluation of ∇s. This picture corresponds to a skeleton connecting extrema with relatively straight segments. In the "stiff" approximation, equation (10) becomes The differential length is then only the function of γ times L.
The "stiff" approximation can be looked at from another perspective. By definition at a point on a local critical line, two of the characteristic directions defined for the field, namely, the direction of the gradient, ∇ρ, and one chosen eigen direction of the Hessian, H, must coincide. But the direction of the critical line itself, given by ∇S i × ∇S j , is not, in general, aligned with the gradient of the field. Local critical lines are not the gradient lines, and in this sense they differ from the skeleton lines defined globally as void-patch intersections ). In the "stiff" approximation, however, it is parallel to the gradient. Figure 1 shows the details of the calculations for the high-resolution segment of the 2D field. Thus, the essence of the stiff approximation lies in the assumption that the mismatch between the critical lines and gradient directions is statistically small. As Figure 2, which contains an extended view of the same field, illustrates, this assumption holds particularly well for the primary critical lines which more closely correspond to the intuitive picture of sharp ridges and deep valleys. Indeed, at a primary line the gradient points to the least curved direction, i.e, in some sense in the direction in which the changes of the field properties are the slowest. Therefore one can expect that this is the direction in which the condition of criticality will be maintained, i.e which the critical line itself will follow. Figure 2 shows that the primary lines start to deviate from the gradient flow mostly towards their end points when the curvature of the field along the line becomes comparable in magnitude to the transverse one. Secondary critical lines are much less certain to follow the gradient, sometimes exhibiting a "sliding" behaviour, on occasion almost orthogonal to the gradient, as a loop-like secondary line near the right saddle in Figure 2 exhibits. So the stiff approximation for the secondary lines should be taken with more caution, although we include them for completeness.
The stiff approximation provides a framework to compute the total differential length of the critical lines and the local skeleton almost completely analytically. In the next two sections we will carry this calculation in two and three dimensions and argue that it can straightforwardly be extended in N dimensions (see Appendix A). In what follows we shall omit in the derivation for brevity the 1/R * (in 2D) and 1/R 2 * (in 3D) factors, but keep in mind that all the length quantities below scale accordingly. In section 4.4, after the computational machinery is developed, we return to the role the third derivative may play in description of the critical lines. Figure 2. An example of a generic patch of a 2D field. The underlying isocontours correspond to the density field. The thin gold lines are the gradients lines of the field. The blue lines is the local set of critical lines, given by the solution of S ≡ ( H · ∇ρ) × ∇ρ = 0. Primary lines are shown in solid and the 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 (Jost 2008) at the intersection of a peak patch and a void patch. The primary local lines follow fairly well the gradient lines, noticeably near the extrema, 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 the everywhere 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 trough) 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 of this is the loop on the right hand side. A zoomed view of the area left to the top maximum is shown in Figure 1.

CRITICAL LINES OF 2D FIELDS
Even though the large scale structures of the universe are three dimensional, other important observed data sets involve 2D maps such as the cosmic microwave background or lensing convergence maps. Hence analyzing the local statistical properties of filaments in two dimensions is astrophysically well motivated. The 2D case is also a convenient starting point to introduce the details of the calculations that can be generalized to 3D and higher dimensions. The 2D case affords several simplifications over the 3D case. In 2D, S is a (pseudo) scalar function and its zero level, orthogonal to ∇S, determines the critical lines. The expression for the differential length simplifies to There are just four types of critical lines: two primary, the skeleton and the anti-skeleton, and two corresponding secondary ones. The classification of the 2D critical lines is summarized in Table 2. We shall focus on the most interesting primary lines in the main text, leaving the secondaries to the Appendix B. In Figure 2 the critical lines of different types are shown for an example generic patch of a 2D field.

The differential length of the critical lines of 2D fields
For 2D Gaussian fields, the calculation of the length of the critical lines can be carried almost completely analytically in the stiff approximation.
3.1.1 Direct derivation in the field's frame Let us first proceed in the original coordinate frame. Defining the stiff approximation to ∇s involves only up to second derivatives of the field and equation (15) becomes explicitly where the second derivatives are described using u = −(x11 + x22), w = (x11 − x22)/2 and x12. Let us integrate over x12 using the δD-function, which leads to a substitution x12 → (2x1x2w)/(x 2 1 − x 2 2 ) with the Jacobian |1/(x 2 1 − x 2 2 )|. Then equation (18) where Let us now substitute 3w , noting thatw 2 = (w 2 + x 2 12 ) , to obtain The integration over the first derivatives is now easily performed in the polar coordinates of the x1, x2 plane to give This is the final integral form which can be easily investigated in the u,w plane.

Derivation in the Hessian eigenframe
To generalize the derivation to higher dimensions we note that one can just perform all the calculations in the Hessian eigenframe. We shall denote all quantities evaluated in the eigenframe with tilde, e.g.,x11(= λ1),x22(= λ2),ũ,w,x1,x2. What must be taken into account is that, in general, these quantities are not Gaussian random variables (while the corresponding ones in the fixed frame are), since the transformation from the fixed to eigenframe is non-linear. The Gaussian nature is only preserved forũ,x1,x2. In the Hessian eigenframex12 = 0. From equations (16-17) In equation (15) the averaging is now carried over the distribution of the eigenvalues with the measure π(λ1−λ2) (Doroshkevich 1970) that accounts for eigenvalues being sorted, λ1 λ2: or in terms ofũ,w In the argument of the delta-function in equation (26)w can be zero only at special field points , not at a generic point on a skeleton. So vanishings requires eitherx1 = 0 orx2 = 0 which describes, as expected, that in the Hessian eigenframe one of the component of the gradient vanishes on the critical line. Since we have already chosen the coordinates so that the direction "1" is aligned with the largest eigenvalue and the critical lines can go in both eigen-directions, these two possibilities add up: 4 Note that |w −ũ/2| = |λ1| = |x11| and |ũ/2 +w| = |λ2| = |x22|. That is, the length of the critical lines per unit volume is given by the average absolute value of the Gaussian curvature of the field in the space orthogonal to the skeleton, given that in stiff approximation the direction of the skeleton is assumed to coincide with the gradient of the field. The reason for this is clear -the higher the curvature, the closer the next neighbouring segment of the skeleton can be, thus increasing the flux i.e. the length per unit volume. If we replacew → −w in the second integral, we return to the formula (23) with integration over both positive and negativew. The integrated length of the critical lines is reduced to equations (27) and (28) are the results of the stiff approximation for the threshold dependent differential and the integrated lengths of the critical lines in 2D respectively.

Primary critical lines in 2D: Skeleton and anti-Skeleton.
The local skeleton is the subset of all the critical lines, which includes the parts that appear as the ridges in the field profile, rather than the valleys. This subset is described by the constraints that the skeleton lines should go along the largest eigenvalue λ1 and, in addition, that this direction has the smallest curvature, |λ1| |λ2|. The anti-skeleton is a mirror structure describing the valley of the field and in all the results can be obtained by replacing η → −η in the formulae for the skeleton.
To derive the expression for the skeleton differential length let us return to equation (27). The critical lines with ∇ρ aligned with the largest eigenvalue direction havex2 = 0. Thus, only one term is selected by the δD-function: it is ∝ 2|λ2| = |2w +ũ|. The magnitude restrictions translates intoũ 0, thus This result should not be confused with equation (28), wherew is integrated over full range of negative and positive values and which is strictly equivalent to equation (27), counting critical lines aligned both with the lowest and the largest eigen-directions. Performing the last two integrals one obtains for the differential length in closed form and for the integrated skeleton length 5 Note that modulo the stiff approximation, equation (31) gives a universal, spectral parameter independent, scaling. Figure 3 demonstrates the threshold behaviour of the differential lengths for several values of the spectral parameter γ.
The most important and robust result of our theory is the behaviour of the differential length at high density thresholds It represents a bias similar to the one found in Kaiser (1984) for the clustering of high critical points -maxima. According to the latter, the number density of peaks in regions above high thresholds is higher than on average. Similarly, the length density of critical lines above high threshold is enhanced relative to the mean. From the point of view of measurements, perhaps a Right: ∂L/∂η/P (η) (solid) and its asymptotic behaviour (dashed) in 2D for the same spectral parameter values more interesting quantity than the differential length is the length per unit volume within the regions of high excursions of the field L(> η). In terms of the cumulative length given by equation (9), L(> η) = L − L(η). Its asymptotic behaviour at high η for the skeleton is found by direct integration of equation (32) L The first factor here is the fractional volume occupied by these high excursions of the field. Note that, at large η the differential length divided by the PDF scales like ηγ/R * = η/R0 once the proper scaling with 1/R * is introduced. Hence the differential length as a function of η together with the total length give access to two characteristic scales R0 and R * . See Appendix A for a general proof of this result in N dimensions. The threshold dependence of the statistics of critical lines in the stiff approximation is determined solely by the spectral parameter γ. In the limit γ = 0, when the distribution of the second derivatives of the field ρ is completely independent on the threshold, the length of the skeleton per unit volume within the regions with ρ/σ0 in the interval η, η + dη is just proportional to the fraction of the unit volume that these regions occupy. Completely generally, for any type of critical line, When γ → 1 the trace of the Hessian u becomes uniquely determined by the field level η (recall equation (6)). For over-dense regions with positive η equation (32) is exact for γ = 1, while no skeleton exists in under-dense regions in this limit. Near zero (mean density) threshold the dependence of ∂L skel /∂η is Its details, in particular a step-like cutoff at negative η when γ → 1, are sensitive to the definition of the primary lines. In under-dense regions with large negative densities the skeleton is exponentially suppressed. Starting from equation (30) with η → −η for anti-skeleton, we obtain for the union of both primary critical lines with twice the integrated length This function is now symmetric in η with the skeleton providing the dominant contribution described by equation (32) in over-dense regions of space, and the anti-skeleton dominating the under-dense regions. Near the mean, zero, threshold of the field, both critical lines are present

Secondary critical lines in 2D
Secondary critical lines do not allow for a full analytical treatment and are investigated in Appendix B. They are particularly important near zero threshold, since at this transitional regime the exact behaviour of primary or secondary lines depends significantly on our somewhat arbitrary separation of the critical lines in types. In this paper we are tracking the skeleton -density ridges -as primary lines emanating from the maxima, until the largest eigenvalue ceases to be the shallowest. Alternative definition may, for example, somewhat extend the skeleton at the expense of secondary lines at lower densities as long as all the eigenvalues transverse to the gradient are negative, i.e until λ2 becomes positive. As an advantage, the differential length of the skeleton and the corresponding secondary lines defined this way would not exhibit inflections at low densities that can be seen in Figures 3 and B1 for high γ's. But the downside is that then one looses the ability to describe the primary lines analytically in a closed form. At the high density excursions the properties of the skeleton remain robust with respect to the variations in their exact definition. However the important advantage of the definition of the primary lines adopted in this paper lies deeper. The magnitude of the eigenvalue along the direction transverse to the gradient is connected to the stability of these trajectories near the critical lines and to their possible bifurcations. This is discussed in part in Section 5.3.
Let us summarize the results for the total set of critical lines, primary and secondary combined, which are, of course, universal whatever the definition of the separate types. Summing up the results of this Section with the corresponding ones in Appendix B The full behaviour of the total differential length is presented in Figure 3. One should note the linear asymptotic behaviour at high density levels and the regular quadratic behaviour near zero density threshold 6 . Finally recall that in the section 3 we have omitted almost everywhere a 1/R * factor for the quoted lengths and differential lengths.

CRITICAL LINES OF 3D FIELDS
In three dimensions, we carry the computations directly in the eigenframe of the Hessian, following closely the derivation of Sections 3.1 and 3.2. We present the formalism first for all the critical lines and then narrow our focus to the primary ones.

The length of the critical lines of 3D fields
In 3D, let us use the variablesũ = −(λ1 + λ2 + λ3),w = (λ1 − λ3)/2,ṽ = (2λ2 − λ1 − λ3)/2. In the Hessian eigenframẽ and f ∇s In the eigenvalue space the measure is 2π 2 |(λ1 − λ2)(λ2 − λ3)(λ3 − λ1)| and the eigenvalues are considered sorted. For sorted eigenvalues the choice of the directions has been fixed and the (s 2 , s 3 ), (s 1 , s 3 ) and (s 1 , s 2 ) pairs of surfaces describe different possibilities for the critical line. Those choices add together in the average integrated length. Using the variablew,ṽ the condition of eigenvalues being sorted isw 0, −w ṽ w. Let us consider the critical lines that are the intersections of (s 2 , s 3 ). Their differential length is given by 6 For all γ but γ = 1, for which c 0000 RAS, MNRAS 000, 000-000 The local theory of the cosmic skeleton 11 The skeleton ∂L skel /∂η/P (η) (solid) and its asymptotic behaviour at high density thresholds (dashed) in 3D. The anti-skeleton is described by the curves symmetric with respect to a reflection of η. Right: ∂L/∂η/P (η) its asymptotic behaviour for the total set of critical lines. The spectral parameter values are (from bottom to top at high η) γ = 0.3, 0.6, 0.95.
Integration over δD(s2) and δD(s3) leads to the only possibilityx2 = 0,x3 = 0. That is, the choice of the surface S2 and S3 in the Hessian eigenframe describes the skeleton along which the gradient is aligned with the direction 1, correspondent to the largest eigenvalue, while in the directions 2 and 3 the components of the gradient of the field vanish. Withx2 =x3 = 0 we get a simple expression for while the subsequent integration overx2 andx3 using δD-functions and afterwards overx1 gives Notice again that what the integrand involves the Gaussian curvature in the direction orthogonal to the gradient, which in stiff approximation is the direction of the filament itself. The contributions of the critical lines directed along the second and third eigen-direction is given by similar considerations and are added together when all critical lines are considered. Changing variables one finally obtains ∂L ∂η = 3 3 5 5/2 expˆη 2 /24 while the integrated length is The equations (48) and (49) account for all the critical lines. In Figure 4 (right panel) the results for 3D critical lines are plotted while the discussion of the corresponding asymptotics is given in the Appendix C. We shall now turn our attention to the study of the primary lines and, in particular, the 3D skeleton that delineates the over dense filamentary structure and is of more direct observational interest.

Primary critical lines of 3D fields: Skeleton and Anti-Skeleton
The subset of critical lines identified with the skeleton correspond to the lines with the gradient aligned with the largest eigenvalue λ1 while having λ1 + λ2 0. In equation (48) such lines are described by the first term ∼ |λ2λ3|. The differential 7 One should note the correspondence with the well-known result for the number density of extrema of the field (Bardeen et al. 1986) which is determined by the mean three-dimensional Gaussian curvature |λ 1 λ 2 λ 3 |. . Left: Integration zones in theṽ,ũ plane for the 3D skeleton analysis. Variables are given in units ofw.ṽ varies from −w to +w, whileũ must be greater than 1 2 (ṽ + 3w). In the allowed shaded region 0 > λ 2 λ 3 everywhere. Horizontal lines mark the further subdivision of the integration space if the order of integration is changed according to equation (53). Right: Integration zones in the (ṽ,w) plane afterũ has been mapped to the [0 − ∞] interval. Variables are given in the units ofũ. The lower triangular zone corresponds to the semi-open rectangular band above the red dashed line in the left panel. In this region the integrand is given by the first term of equation (53). It dominates the high η asymptotics.
length of the skeleton is then The integration inṽ-ũ plane is limited to the regionũ > 1 2 (ṽ + 3w), as shown in the left panel of Figure 5. The integrated length of the skeleton is that is, one expect on average one skeleton line crossing a random ≈ (5R * ) 2 surface element. The results of integration of equation (51) are presented in the left panel of Figure 4.

Asymptotic behaviour at γη → ∞
To study high η asymptotes it is useful to change the order of integration to have theũ integral as the outmost one. The inner integration inṽ-w plane is then carried out over the region shown in the right panel of Figure 5.
The last term is exponentially suppressed as η → ∞ while the first one gives The leading quadratic and the next linear terms can be recovered found by replacingũ → γη in the pre-exponential factor and treating the exponent as the δD-function. A more detailed asymptotic study of this Laplace-type integral is required to recover the third-order constant term, that also contributes to the accuracy of the expansion at the level demonstrated in Figure 4. One finds that in the leading order in η the skeleton has the differential length growing as (γη) 2 (see also Appendix A) and involves, as expected, a third of all the critical lines (compare with Appendix C) in the regions of high excursions concentrated around the maxima of the field. However, at intermediated thresholds, the skeleton constitutes more than a half of all critical lines, highlighting enhanced importance of the filamentary dense ridges among other critical lines.

Power series at η → 0 and Hermite expansion
Using two alternative series representations of the shifted Gaussian form that encodes the dependence of the skeleton on the threshold η we obtain either power series or Hermite 9 (Novikov et al. 2006) expansion of the differential length where and These two expansions are similar but distinct. The power-law expansion is suitable for an accurate analysis of the differential length near zero threshold for all γ < 1. On the other hand, the expansion in orthogonal Hermite polynomials is useful as an approximation over an extended range of thresholds. Both series are improper for γ = 1. Although these coefficients can be computed analytically, their expressions are too cumbersome. Instead, we plot several leading ones in Figure 6. Remarkably, the power in Hermite expansion is concentrated in a few low order terms, in particular, k = 0, 1, 2, 3 for the skeleton, with subsequent terms forming a slowly decaying oscillating series. This finding confirms in 3D the conjecture of Novikov et al. (2006). The contribution of the first three most dominant terms, P 2 k=0 B k H k (η) = 0.0462 + 0.0751γη + 0.0464γ 2 (η 2 − 1) has the same structure and remarkably similar coefficients as the high η asymptotics of equation (54) which evaluates to 0.0477 + 0.0852γη + 0.0531γ 2 (η 2 − 1). This explains why the high η asymptotics provides a visually good fit through all thresholds when γ is not too high. At γ → 1, the oscillatory tail of Hermite series provides the correction that reflects the irregular nature of the expansion in this limit.

Primary critical lines of 3D fields: Inter-Skeleton and the overall behaviour
The intermediate primary critical lines are associated with saddle-like regions where the largest eigenvalues in magnitude are λ1 0 and λ3 0, and have opposite signs, and the shallowest direction aligned with the gradient is the second one with −λ1 < λ2 < −λ3. Their appearance reflects the complexity of critical lines in space of more than two dimensions.
The differential length of the intra-skeleton computed in the stiff approximation is presented in Figure 7. The conditions for intermediate lines are prevalent for the regions of the field of moderate values -within 2σ (|η| < 2 of the zero mean for γ > 0.6. Although the occurrence of the intra-skeleton within these regions is never large (∂L/∂η/P (η) is relatively small), the regions corresponding to a near mean density occupy large fractions of the total volume, and as the result the total length of the intermediate skeleton is almost twice that of the skeleton or the anti-skeleton: It constitutes nearly a half of the total length of the primary critical lines At high |η| thresholds, in very dense regions near maxima or under-dense regions near minima of the field, the intermediate skeleton is rare. The total set of the primary critical line is even more than the skeleton dominated by the low order terms in Hermite expansion. Indeed, Figure 6 demonstrates that just the first two terms (odd orders are absent due to symmetry) in Hermite series are dominant, P ∞ k=0 B k H k (η) ≈ L prim tot`1 + 0.340γ 2 (η 2 − 1)´.

Validity of the stiff approximation
Let us consider the opposite to "stiff" regime, when the derivatives of the Hessian dominate the ∇s, Although not natural for cosmology-inspired spectra, such a situation arises when the power spectrum has an extended short wave tail with spectral index 10 n between −9 and −5. Such spectra have smallγ,R R * and there are many inflection points of the field per extremum. Interestingly, this regime also automatically means that the correlation between the gradient and third derivatives of the field is small.
Using the Hessian eigenframe formalism, we can obtain the important results without explicit computation of the differential length. Let us focus on the critical lines corresponding to the first eigenvalue. Equations (43) for S-surfaces gives rise to two δD-functions, δD(2wx1x3)δD((w −ṽ)x1x2) that after integration over the transverse gradient componentsx2 andx3 enforcex2 =x3 = 0, with the Jacobian factor 1/|2w(w −ṽ)x 2 1 |. The length element in this frame obeys 10 In a cosmological framework this takes place when the density field with n < −1 spectrum is smoothed with a top-hat window.
c 0000 RAS, MNRAS 000, 000-000 The local theory of the cosmic skeleton 15 where the last expression defines the ψ(x klm ) function. The differential length is now given by The last integral, with P1 given by equation (D6), is a function ofγ only. The first term shows that, since the integrand prefactor is independent onũ, the differential length does not depend on the threshold η at large η (it does at small η only because of non-trivial integration boundaries dependent on the exact type of critical lines). This is not surprising, since in this limit, there is little link between the skeleton length and the second derivatives, the only ones that are correlated with the field value. Such threshold independent behaviour is not observed in simulations with cosmological spectra, which argues once again for the statistical validity of the "stiff" approximation.

Measurements
In this section we compare the predictions of the local theory in stiff approximation with the measurements of the statistical properties of the critical lines done on realizations of the Gaussian fields with different power spectra. We perform the measurements on critical lines found according to the global definition. The measurements are carried as follows: a set (typically ∼ 100) of scale-invariant Gaussian random field of a N-D maps (typically 1024 2 ) or cubes (typically 256 3 ) is generated with a given power index of n = 0, −1 or −2. The N-D cube is then smoothed via convolution with a Gaussian kernel of width 6 pixels. The spectral parameters, γ,γ etc... are computed through the second moments of the derivative of the smoothed field. The set of critical lines is then extracted as the intersection of the peak patches and void patches (see  for details). In Figure 8 an example realization of the primary critical lines in 3D cube is shown. Since the algorithm produces a set of segments describing those critical lines tagged by the underlying (smoothed) density field, it is straightforward to compute the total and differential length per unit volume of the whole set. The differential length per unit modulus gradient is extracted by tagging the critical lines with this modulus (obtained via Fourier transform differentiation) and proceeding as before. Finally, the curvature of the skeleton is measured by computing the local curvature of a set of adjacent segments via finite difference.
Let us emphasize that these measurements correspond to properties of the global skeleton, whereas the theory developed in this paper is focused on the local skeleton. Hence even more remarkable is the match between the measured and the theoretical differential lengths for all values of γ, that is exhibited in Figure 9. This accuracy should be considered as indicative of the correspondence between the stiff approximation to the local theory and the global set of critical lines.

OTHER STATISTICS AND SPECTRAL PARAMETERS
In the previous sections, the emphasis has been on the differential length of the critical lines as a function of the excursion in density. As argued in  and demonstrated here, it provides means of constraining the shape parameter, γ. Let us now explore other statistics which will allow us to constraint other shape parameters. In particular, let us demonstrate that the differential length as a function of the excursion in the modulus of the gradient of the density and the differential curvature depend on the second shape parameter,γ. Finally, we investigate the number density of singular points on the critical lines.

Differential length versus the gradient modulus
The differential length of the skeleton with respect to the threshold η carries information on the spectral parameter γ thanks to the correlation between the value η and the Hessian of the field. In the stiff approximation the Hessian curvature completely determines the length of the critical lines. For the exact formulation, the length also depends on the third derivatives, that are correlated with the first derivatives via the parameterγ. Thus, measuring length as a function of the modulus of the gradient should carry information onγ and provide an estimate of an impact the third derivatives have on the length statistics of the critical lines.
To demonstrate the dependence of the skeleton length on the gradient of the field in "stiff" approximation let us return to equation (45) which we take integrated over all density thresholds. As before, we perform the integration over the δ-functions that enforces alignment of the gradient with the first eigen-direction, x2 = x3 = 0, however this time we do not integrate over but rather take the differential of the result with respect to x1. Noting that |x1| = X ≡ p x 2 1 + x 2 2 + x 2 3 . we obtain in place of equation (47) where the last integral does not depend on X. Dividing by the integrated length, L = R ∞ 0 ∂L/∂XdX, and generalizing the  result to fields in arbitrary N dimensions we conclude that The exact dependence of the differential lengths will deviate from this form in aγ -dependent way. It is natural to Gaussian random fields with scale invariant spectra smoothed over 7 pixels. Right: value of the fit parameters c k (see equation (67)). Note that c 0 = 1.
parameterize such deviation expanding the true statistics in Hermite series around the stiff approximation This choice of expansion is dictated by the orthogonality of the Hermite polynomials with the weight ∝ exp[− √ N X 2 /2] on the interval X 0. Thus, c0 = 1. If the deviation from the stiff approximation is small, one expects the expansion to be dominated by the n = 0 term, while the subsequent terms should quickly fall in a orderly fashion.
To gain understanding on how the coefficients c 2k (γ) behave withγ, let us consider again the lax situation, opposite to the stiff case, when the third derivatives of the field dominate the length statistics. Our starting point is equation (64) which has the following structure when we consider the differential length with respect to the |x1| = X

∂L ∂X
whereP1 is given by equation (D6) with the dependence on u1 factored out. The difference with the stiff approximation is large even forγ = 0 as the gradient's dependence becomes ∝ X 2 exp[−3/2X 2 ] in place of the stiff scaling ∝ exp[−3/2X 2 ]. Using now this factor as the weight, forγ = 0 we expand the expression in the brackets in generalized Laguerre polynomials. The expansion coefficients are of the formγ 2k exp[−3u 2 1 /2] P k m=0 dmγ 2m H2m( √ 3u1); denoting the result of the integration of the expansion coefficients and all of the residual factors over the third derivatives by Ψ k (γ) we obtain where, again, Ψ0(γ) = 1. With the help of the relation between the Laguerre and Hermite polynomials we can cast equation (69) in the form of equation (67) 1 L

∂L ∂X
The coefficients c 2k (γ) are In particular, in the limitγ → 0 the first two coefficients remain finite and of equally significant magnitude c0 = 1, c2 = √ 2, while all the other ones vanish. Figure 10 and Figure 11 present the measurements of ∂L/∂X/L in 2 and 3D respectively, together with the corresponding coefficients given by equation (67). It is found that these coefficients are significantly smaller in two dimensions, a clear indication that the stiff approximation holds better in 2D. Gaussian random fields with scale invariant spectra smoothed over 7 pixels. Right: value of the fit parameters c k , Note its faster convergence combined with a larger amplitude relative to the 2D case.

Statistics of the curvature of the critical lines
The local curvature, κ, at a point on a curve specified by the tangent vector u = dr/dt is determined by the acceleration of the tangent vectoru ≡ du/dt = u · (∂u/∂r) transverse to the curve direction: Importantly, the curvature does not depend on parameterization t, nor on normalization of the tangent vector u. In the local theory, the tangent vector to a critical line is orthogonal to ∇s i (x k , x kl ) and can be taken to be so the curvature κ is the random quantity which involves the derivatives of the field up to fourth order, The curvature of the critical lines fundamentally reflects the derivatives of the field higher than the second. If they are neglected, the curvature is identically zero. Explicitly, the contributions that do not involve higher derivatives, in 2D (2D) (u ×u) 2 = 4x 2 1 x 2 2 λ 4 1 (λ1 − λ2) 6 λ 4 2 + . . . , vanish when the correspondent critical line conditions x2 = 0 or x2 = x3 = 0 are applied 11 . The integrated curvature over the length of the line, C = R κdL is a useful dimensionless characteristics of the overall extend a line is curved. We have seen that the critical line length in volume dV is dL ∝ |∇s|δD(s)dV and dL ∝ |∇s i × ∇s j |δD(s i )δD(s j )dV in 2D and 3D respectively. Averaging over statistical distribution in regions above threshold η we obtain the mean density of the integrated critical line curvature C = dC/dV where the integration is carried over all the derivatives up to the fourth order. The required joint probability function is given in equations (D19) and (D21). Let us consider 2D case and estimate the curvature by using stiff approximation for the tangent vector u while following its local variation which involve higher derivatives of the underlying field. In the Hessian eigenframe, assuming the skeleton lies along 1, if u is approximated by its stiff counterpart we have: so that (taking into account the measure in the eigenframe and the δD function of S in x2): the last evaluation being done for the primary critical lines. We have measured the mean curvature of the skeleton lines at the threshold η, in simulations of the Gaussian random fields of different spectra, using the global skeleton techniques. Figure (12) displays the 2D (left panel) and 3D (right panel) results in terms of the curvature radius, Rcurv = 1/ κ . The measurements show that for the spectra we consider, the averaged curvature is insensitive to the density threshold for low-to-moderate threshold values showing a plateau in the interval −2 < η < 2. This indicates that in this regime the curvature of the skeleton does not depend on γ, but rather onγ and perhapsγ. It follows that in 3D the critical lines are relatively more wiggly than in 2D. If we use the lower value of R 3D curv as a guidance, it seems the stiff approximation is less accurate in 3D that in 2D, as could be expected. The stiff estimate (81) gives the threshold-averaged mean density of the integrated curvature C stiff and, using equation (37), the mean curvature radius, R stiff curv , as This result captures the qualitative dependence onγ observed in simulations, but is a factor of three smaller in the magnitude of the curvature radius. This shows that the global skeleton, used in the numerical measurement, is notably straighter than the local critical lines, although the dependence of curvature on the spectral parameters seems similar. Note in closing that in 2D, (resp. 3D) the knowledge of the differential length, curvature (resp. length, curvature and torsion) corresponds to an exhaustive global statistical description of the critical lines.

Singular points of critical lines
Let us now ask ourselves the following question: are there any special points along the skeleton? The obvious ones are the extrema of the field itself where critical lines intersect. Beyond this, one can anticipate two other types of singular points. 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. From the point of view of the theory of random fields, the frequency of such points is an interesting venue: indeed we expect that steep power spectrum present relatively more bifurcation points asR, the distance between inflection points (see Sec. 2.2), becomes much shorter than R * , the distance between extrema. In an astrophysical context, the statistical properties of the first type of points, and in particular their clustering properties are of interest for understanding the geometry of galactic infall, which in turn is believed to play an important role in defining the morphological properties of galaxies. The multiplicity of the maxima (i.e. the number of connected skeleton segments) is also of interest in the context of galaxy formation and feedback. In more abstract spaces, such as position-time, identifying bifurcations is important to pin down merging events (see e.g. Hanami (2001) and Appendix A).

Defining the skeleton singular points
Formally a singular condition along the skeleton occurs when at some point the determination of the critical line direction fails. It means that at this point the matrix ∇ k S i of equation (2) has more than one distinct right null-vector, or, equivalently, all M k defined by equation (A7) are zero. The only case when it happens exactly is at the extremal points of the field ∇ρ = 0. There are no other formal singularities on the local critical lines, since when ∇ρ = 0, the requirement M k = 0 sets N relations between the field gradient, second and third derivatives which have vanishing probability to be simultaneously satisfied along a line in a random field.
The failure of the formal definition to identify all the physically interesting situations primarily reflects the inadequacy of the local skeleton construction, which only utilizes locally quadratic approximation to the field, to map the field near the singular points. 12 Figure 13 gives a 2D example. In 2D, the critical lines are zero levels of the scalar S-function, while ∇S = 0 at the extrema of S field. In Figure 13 the region where a ridge splits into two is shown. One expect two critical lines cross there, with three branches following the ridges, and one following the through between two of the split branches. Instead, the locally defined critical lines are not allowed two join at the bifurcation point since the formal condition ∇S = 0 is satisfied just off S = 0 contour, rather they artificially reconnect near the bifurcation point into two non-intersecting segments.
We conjecture that the critical lines experience a qualitative change in behaviour in the vicinity of the points where either the Hessian eigenvalue of the orthogonal to the gradient direction vanish, or becomes equal to the one along the gradient. Namely, if, for definiteness, ∇ρ is taken to be along the first eigen-direction, λ2 = 0, or λ2 = λ1. We call the first case the "sloping plateau" as it designates the entering of a flat region, and the second, tentatively, the "bifurcation" as it designates the places of possible reconnection of critical lines. In particular, at the λ2 = λ1 points most of the transitions from primary to secondary behaviour take place. Remarkably, these special points on the critical lines are recovered by the formal singular condition |M k | = 0 if ∇ k S i is evaluated in the stiff approximation. As given in equation (A12), along the ND critical line defined by x2 = . . . = xN = 0, |M stiff | = x N −1 1 Q i>1 λi (λ1 − λi) = 0 gives rise to three classes of situations: (i) x1 = 0 corresponding to extremal points; (ii) one of λi = 0 corresponding to slopping flattened tubes; and (iii) one of λi = λ1, corresponding to an isotropic bifurcation.
Since it is beyond the scope of this paper to develop the full theory of these special points, we will focus here on their number density for isotropic 2D Gaussian random fields, leaving more detailed investigation to future work.

Number density of the singular points of the 2D critical lines
In 2D, the skeleton's singular points correspond to points where S k ≡ ∇ k S = 0. The number density, nB(η) of singular points below the threshold η is equal to The simplest case of the skeleton singular points ∇S = 0 are, according to equation (8), the extrema of the field itself, x1 = x2 = 0. Indeed when both x1 and x2 vanish |det (∇ k ∇ l s) |δD(s1)δD(s2) = |x kl |δD(x1)δD(x2) , which is exactly the integrand involved in the number density of extrema of the field. The extrema number densities, for reference, are given in 2D by (Longuet-Higgins 1957) The singularity of the extrema from the points of view of the critical line theory is manifest in the fact that at extrema several critical lines intersect. The gradient of S, evaluated in the stiff approximation, in the Hessian eigenframe has the components and involves only second derivatives of the field. Remarkably, within this approximation, there are new singular points that lie on the (local) critical lines. The reason is that among two conditions needed for ∇S to vanish, one is already automatically satisfied by being on a critical line.
To be specific, let us consider the critical line that corresponds to the x2 = 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, x1 = 0, but also in two other cases, namely λ2 = 0 or λ2 = λ1, that we conjectured to be of interest. The first situation, the sloping plateau with a flat transverse gradient, only occurs on secondary critical lines since it implies λ1 > 0, and corresponds tǫ hence The second situation (isotropic Hessian) corresponds tǫ Both G F B (γ) and G I B (γ) are weak functions ofγ of order unity. The mainγ dependence ∝γ −2 reflectsR as the fundamental scale for the singular points.
We note that the number density of the "sloping plateaux" is proportional to the density of the saddle points, hence this type of singular points is predominantly concentrated near mean field values (small η). In contrast, the number density of "bifurcation" points is proportional just to the PDF of the field and, hence, the bifurcation points are as frequent in the regions of high field values as in the low ones. This may provide explanation for the observed insensitivity of the curvature of the skeleton to the threshold, if we conjecture that most of the curvature accumulates near the "bifurcation" points.

CONCLUSION & PERSPECTIVES
The filamentary structure is a dramatic feature of the observed or simulated Cosmic Web. This paper investigated how the set of critical lines of a given field corresponds to an intermediate representation of the field, which is more extended than the knowledge of the critical points, but nevertheless much more compact than the field itself. It introduced the stiff approximation, which states that the tangent vector to the critical lines only involves up to the the second derivative of the fields. Within its framework it has been demonstrated that, for stationary Gaussian random fields, ergodicity allows one to recast the description of the ND critical lines into a point process, which only involve the first spectral parameter, γ, when considering the differential length as a function of the contrast, and the second spectral parameter,γ, when considering it as a function of the modulus of the gradient. The former probability distribution was shown to involve the average flux of the Gaussian curvature of the 1D sections. In turn, these averages can be carried out analytically almost to the last integral in 2D and 3D, and provide simple asymptotics at large and small contrast. The detailed contribution of all types of critical lines as a function of thresholding was described. The main results of this investigation corresponds to equations (27) and (28) for the differential length of the skeleton and the total set of critical lines in 2D and equations (48) and (51) in 3D. Their generalization to N dimensions is given by equation (A18) in Appendix A. Table 3 summarizes the average integrated fluxes (i.e length per unit volume) of the critical lines. For instance in 3D one expect on average one skeleton line crossing a random ≈ (5R * ) 2 surface element. These findings were illustrated on scale free power spectra with spectral parameters which are relevant to cosmology 13 . The prediction of the stiff approximation was checked against measurements for global skeletons ) on realizations of these fields in two and three dimensions and was found to be in good qualitative agreement. The differential curvature of the corresponding lines was also measured (section 5.2) and the corresponding radii were found to be ≈ 8R * and ≈ 2.5R * near η = 0 in two and three dimensions respectively. Hence an access to both the curvature and the length of the skeleton provides the means of constraining two shape parameters, γ andγ. The stiff approximation is also implemented to compute the differential curvature in 2D. Finally (section 5.3), the stiff theory of the singular points of the critical lines was laid out in general, identifying generically three types of points: critical points of the underlying field, bifurcation points and slopping plateaux. Again, the stiff approximation provide means of computing the number density of these points. Appendix D derived the general joint probability of the field and its successive derivative in arbitrary dimensions, which come into play when computing these higher order statistics.
Clearly the formalism developed in this paper will be useful in the context of the upcoming surveys such as the LSST, or the SDSS-3 BAO surveys since it yields access to the shape of the power-spectrum without artifacts related to varying light to mass ratio. For instance,  first applied the corresponding theory to the SDSS-DR4 catalogue in order to constraint the global dark matter content of the universe, since the cosmological parameters are directly a function of the spectral parameter, γ. Its application to CMB related full sky data, such as WMAP or Planck should provide insight into, e.g. the level of non-Gaussianity in these maps (see SPCNP for a discussion). Similarly, upcoming large scale weak lensing surveys could be analyzed in terms of these tools (Pichon et al. 2009).
A natural extension of the theoretical component of this work would be to investigate the properties of the bifurcation points in anisotropic settings and extend beyond the stiff approximation the preliminary results of Section 5.3. This will be the topic of a forthcoming paper. Another natural venue would be to also investigate the statistical properties of, e.g. the peak patch walls (surface, curvature) defined as x3 = 0 in the eigenframe of the Hessian. Eventually, a global theory of the critical manifolds beyond the local approximation should also be developed to provide a framework to study the connectivity of the critical lines.

ACKNOWLEDGMENTS
We thank D. Aubert and K. Benabed for comments and D. Munro for freely distributing his Yorick programming language and opengl interface (available at http://yorick.sourceforge.net/). DP thanks the CNRS (France) for support through a "poste rouge" visiting position during Summer 2007 when this investigation was originated. CP, TS, SP and CG also thank the hospitality of the University of Alberta, and the "Programme National de Cosmologie" for funding. Finally, CP and DP thank the Canadian Institute for Theoretical Astrophysics for hosting the work involved in finalizing this paper. This investigation carried within the framework of the Horizon project, www.projet-horizon.fr.
Local direction of the filament corresponds to the right null-vector δr k of the N − 2 + 1-rank tensor of the derivatives of S A non-trivial solution of this set of C 2 N homogeneous equations generally exists, since the existence of the left null-vector The notion of primary skeleton lines is automatically generalized for N D as the subset of critical lines obeying H · ∇ρ = λ1∇ρ , and λ1 + λ2 0 , where λ1 is the largest and λ2 is the second largest of the sorted eigenvalues. Let us derive the general expression for the statistical average of the flux of the lines arbitrarily defined over the properties of the ND random field by N − 1 equations S i = 0, i = 1 · · · N − 1, where S i (x, x k , x kl , . . .) are functions of the field and it's derivatives. We can shortcut the procedure of flux evaluation by marking each line with one intersection point with a fiducial surface Σ, orthogonal to it, and finding the N-1 number density of the intersection points on the surface Σ = 0. The average number density of the points defined as the intersection of n non degenerate hypersurfaces σ i , . . . , σ N is given by dxd N x k · · · P (x, x k , · · · )δD(σ 1 ) · · · δD(σ n )|det(∇σ 1 , · · · ∇σ n )| .
Let us choose S 1 · · · S N −1 as σ 2 · · · σ n and Σ to be σ N . Expanding the determinant |det(∇S 1 , · · · , ∇S N −1 , ∇Σ)| along its last row we obtain where are the corresponding minors. By design the Σ = 0 surface is to be orthogonal to the line and therefore its normal ∇Σ and the "vector" M ≡ (M k ) are parallel,˛X k M k ∇ k Σ˛= |M| |∇Σ| .
Without loss of generality, we can consider the intersection point to be at r = 0 and take Σ = eΣ · r where eΣ is the unit vector in the local direction of the filament, hence |∇Σ| = 1. The average (N-1)D number density of intersection points on Σ surface that gives us the average flux L is obtained by integrating the volume number density over the coordinate along eΣ, z = eΣ · r, with δD(Σ) in equation (A7) properly counting exactly one intersection per line To apply this general formula to the critical lines one must choose an arbitrary subset of N − 1 linearly independent ∇S i 1 ,i 2 ,...,i N −2 from the set of all C 2 N of them. Note that one can also think of L as the average length of lines per unit volume, which is the interpretation we focus on in the main text.

A2 Stiff critical lines in ND
In the theory of ND critical lines, the N-1 independent functions S i that define the critical condition (A2) acquire the following simple form in the eigenframe of the Hessian of the field Here a is the index of the Hessian eigenvector that the gradient is aligned with, as is obvious from the solution xi = 0, i = a.
In the stiff approximation, the gradients s i k ≡ ∇s i have just two non-zero components, s i a = xiλa(λa − λi) (which vanishes on the critical line) and s i i = xaλi(λa − λi). The vector that determines the direction of the critical line becomes On the critical line, it has just one non-vanishing component which shows that in the stiff approximation we equate the direction of the line with the gradient of the field. Substituting this expression into equation (A9) and integrating over δD(s i ) = δD(xi)/(xa(λa − λi)) we obtain a simple expression for the flux of the critical lines in the stiff approximation i.e the flux of critical lines (or the length per unit volume) is given by the average absolute value of the Gaussian curvature of the field in the space orthogonal to the skeleton. Let us write the probability of measuring the set {λi} as where Qγ is a quadratic form in λi and η which functional form is and Q i<j (λi −λj) is the Jacobian of the transformation to the Hessian eigenframe. Here QN involves polynomial combinations of the eigenvalues of the traceless part of the Hessian (see Appendix D): which can be rearranged explicitly in terms of λs as: It now follows that the differential length of the ND-critical lines is for the stiff approximation: Equation (A18) is the formal generalization of equations (30) and (49). For the ND-skeleton, equation (A18) also holds but the integration region should be restricted to the corresponding condition on the sign of the eigenvalues. Since the argument of Qγ is extremal as a function of η when γη ∼ P 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 (Pichon & Bernardeau 1999). Hence given that Q i<j (λi − λj) is the measure, the only remaining contribution in the integrand comes from˛Q i>1 λi˛∝ (λη) N −1 , and the dominant term at large η is given by where R0 = R * /γ is defined in equation (3).

APPENDIX B: SECONDARY CRITICAL LINES IN 2D
In this Appendix we present a study of asymptotic behaviour of the lengths statistics of secondary critical lines for 2D Gaussian field. Secondary critical lines are the ones that have a gradient of the field aligned with the Hessian eigenvector that corresponds to the largest by magnitude eigenvalue, i.e with the direction of maximum curvature of the field. In 2D, this is the direction of λ2 in the skeleton region, |λ1| < |λ2|, and is the direction of λ1 in the anti-skeleton region. We shall explicitly consider the first type, realizing that the second type is a mirror case with η → −η.
Our starting point is the part of expression (27) that corresponds to the lines where the gradient is aligned with the second eigen-direction, in the region when they are secondary,ũ > 0 The absolute value of the transverse to the gradient curvature 2λ1 = 2w −ũ is evaluated differently forũ 2w andũ > 2w. It is convenient to make the inner integration to be overw, since it can be carried out analytically. The integral splits into two terms where so that finally The integrated length of critical lines L sec is obtained by marginalization over all threshold values η. Performing this integration first  Figure B1. Differential length of the secondary critical lines ∂L/∂η/PDF in 2D for the complete set (solid) and the ones with ∇ρ aligned with λ 2 direction in the skeleton |λ 1 | |λ 2 | region (dashed). Different curves from purple to green correspond to the spectral parameter values γ = 0.3, 0.6, 0.95.
The first three coefficients are If we add all secondary critical lines, the odd power terms of the expansion (B8) cancel, while the even double recovering symmetrical behaviour of the differential length with the threshold. In Figure B1 this behaviour is illustrated. Under our definition of the secondary critical lines, for γ > 1/ √ 2 there is an excess of critical lines near zero threshold. The curvature at η = 0 is positive and diverges in the limit γ → 1 when our series expansion formally fails. This divergence in the second derivative of the differential length is exactly opposite the one the primary lines demonstrate in this limit. We should emphasize, that near η = 0 the behaviour of critical lines of individual type depend significantly on how exactly they are defined.
Direct evaluation of the integrated length gives To study high threshold regime it is advantageous to make theũ integration the outmost one, since it depends on the variable threshold I + IV : II + III : The parenthesis in the last term indicate the substitution that must be performed in the integrand. Right panel in Figure C1 illustrates the integration zones now in the (ṽ,w) plane. Although one can perform theṽ integral analytically and reduce the problem to two-dimensional integration, the resulting expression is too cumbersome. We can obtain useful limits already from unreduced formulae. In particular, at high density threshold, γη → ∞, only the first integral in the term (C4), which containsw,ṽ ∼ 0 neighbourhood, is not exponentially small. Moreover, in the leading order the upper limit of the integral overw can be set to infinity.

APPENDIX D: JOINT DISTRIBUTION OF THE FIELD AND ITS DERIVATIVES FOR A GRF
The joint point distribution functions that are needed for the study of the critical lines in this paper are P0(x, x kl ) and P1(xi, x ijk ), taking into account that for Gaussian random field there is no cross-correlation between odd order derivatives and the field itself or even order derivatives. When considering the curvature of the critical lines, fourth order derivatives, and, thus, more general P0(x, x kl , x klmn ) have to be considered. Some well known results in 2D and 3D are first summarized in section D1. More general results can be obtained by resorting to a general framework which is sketched in section D2 and applied in section D3 for the various cases of interest.
It depends only one a single correlation parameter: γ.
Summary and rescaled forms. We collect all previous results into a normalized form. Using the normalized spectral shape parameters of def. (5) and normalized derivative tensors Xn defined as: Xn = 1 σn ∇ n ρ , i.e. xi 1 ···in = σ −1 n ∂ n ρ ∂ri 1 · · · ∂ri n , one finds that Note that the diagonal entries of Γ are always equal to 1.
Special cases and smaller statistics. Our approach compresses a set of derivative tensors into a set b Γ of symmetric matrices of size m × m , yielding P m (m + 1)/2 invariant scalars. There are two special cases where even smaller invariant sufficient statistics can be found.
First, at angular frequency = 0, the detraced tensors are just scalars so that, for = 0, one has [ b Γ0]pq =˙X q . Therefore b Γ0 actually is a rank-one matrix: b Γ0 = vv † where the entries of vector v are vp = X p . Hence, at the null frequency, we can further compress the m0(m0 + 1)/2 statistics (the non-redundant entries of b Γ0) into m0 scalars (the entries of v). Of course, the = 0 term in the quadratic form also reads: Second, there are several cases of interest where m = 2. This happens for instance at = 0 with derivative orders 0 and 2, at = 1 when considering derivatives of order 1 and 3, at = 2 with derivatives of orders 0,2 and 4, etc. Then, for such an , tr`b Γ Γ −1 ´= tr " »˙a | a¸˙a | bb | a¸˙b | b¸- where a and b are rank-tensors and γ is a scalar. Simple algebra yields that is, a form ubiquitous in this paper. However, an equivalent, more regular form is which has the benefit of stressing that, at such , a sufficient statistic is only made of two invariant scalars, namely a 2 + b 2 and˙a | b¸. In the limit of weak correlation γ → 0, one has, of course, tr`b Γ Γ −1 ´= a 2 + b 2 . An even more symmetric form, which stresses the decorrelation between a + b and a − b is tr`b Γ Γ −1 ´= a + b 2 2(1 − γ) + a − b 2 2(1 + γ) .

D3 Some applications
We now work out these expressions in some cases of interest.
Derivative of orders 0+2 in 3D. The case X = {X0, X2} is the simplest non-trivial case. The theory sketched at sec. D2 applies straightforwardly. In the notations of section D2], we are concerned with frequencies = 0 and = 2 for which, in 3D, w0/g0 = 1, w2/g2 = 15/2 (see table D1). The quadratic form (D15) then reduces to tr`b Γ0Γ −1 0´+ 15 2 tr`b Γ2Γ −1 2´. For = 0, we have here m0 = 2 and we can use the specific form (D17) to work out tr`b Γ0Γ −1 0´w ith [a, b] = [X 0 , X 2 ] = [x, xaa], that is the (normalized) field and the trace of its Hessian. For = 2, we have here m2 = 1: we need only scalars. Following expressions (D15) again, we have b Γ2 = X (2) 2 2 = X 2 2 =x abxab and Γ2 = (−) (2−2)/2 γ2,2 = 1. In summary: The 2D case. The theory applies to isotropic fields in any dimension. We have already provided expressions for the spectral moments (4) and the coefficients w /g of equation (D14). It remains to find detracing coefficients. In the 2D case, the first (ranks 0, . . . , 5) de-traced tensors on R 2 are given by y = y, y i = yi, y ij = yij − 1 2 yaaδij, y ijk = y ijk − 3 4 y aa(i δ jk) , y ijkl = y ijkl − y aa(ij δ kl) + 1 8 y aabb δ (ij δ kl) , y ijklm = y ijklm − 5 4 y aa(ijk δ lm) + 5 16 y aabb(i δ jk δ lm) . (D23) For the correlation between the field and its Hessian, we proceed as above in 3D with w0/g0 = 1 and w2/g2 = 4 given in table D1. Therefore the quadratic form is in agreement with equation (D9). For the case of first and third order derivatives, we read w1/g1 = 2 and w3/g3 = 8 from table D1 so that, similar to equation (D19), one finds with x ijk = x ijk − 3 4 x aa(i δ jk) so that 8x ijk x ijk can be checked to equal Q3 in equation (D12). The d-dimensional case. We outline some results in the d-dimensional case. The de-tracing formulae can be extended to the d-dimensional case but, in this paper, we will content ourselves with the correlations between the field and its Hessian: X = {X0, X2}. Therefore, we need only = 0 and = 2 so that de-tracing remains trivial: the normalized de-traced Hessian given byxij = xij − 1 d δij xaa. Hence for the correlation between the field and its Hessian, we obtain the quadratic form » x xaa which is a straightforward extension of the 3D case of equation (D18). Just recall that γ is now defined in terms of the spectral moments (5) and that de-tracing the Hessian requires a factor 1/d instead of 1/3.