# What makes a cosmic filament?
the dynamical origin and identity of filaments I. fundamentals in 2D

Job Feldbrugge,1 Rien van de Weygaert2
1University of Edinburgh, Higgs Centre for Theoretical Physics, James Clerk
Maxwell Building, Edinburgh EH9 3FD, UK
2Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV
Groningen, The Netherlands E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Cosmic filaments are the main transport channels of matter in the Megaparsec
universe, and represent the most prominent structural feature in the matter
and galaxy distribution.

Here we describe and define the physical and dynamical nature of cosmic
filaments. It is based on the realization that the complex spatial pattern and
connectivity of the cosmic web are already visible in the primordial random
Gaussian density field, in the spatial pattern of the primordial tidal and
deformation eigenvalue field. The filaments and other structural features in
the cosmic web emerging from this are multistream features and structural
singularities in phase-space. The caustic skeleton formalism allows a fully
analytical classification, identification, and treatment of the nonlinear
cosmic web. The caustic conditions yield the mathematical specification of
weblike structures in terms of the primordial deformation tensor eigenvalue
and eigenvector fields, in which filaments are identified – in 2D – with the
so-called cusp caustics. These are centered around points that are maximally
stretched as a result of the tidal force field.

The resulting mathematical conditions represent a complete characterization of
filaments in terms of their formation history, dynamics, and orientation. We
illustrate the workings of the formalism on the basis of a set of constrained
$N$-body simulations of protofilament realizations. These realizations are
analyzed in terms of spatial structure, density profiles, and multistream
structure and compared to simpler density or potential field saddle point
specifications.

The presented formalism, and its 3D generalization, will facilitate the mining
of the rich cosmological information contained in the observed weblike galaxy
distribution, and be of key significance for the analysis of cosmological
surveys such as SDSS, DESI, and Euclid.

###### keywords:

large-scale structure of Universe – cosmology: theory – dark matter

††pubyear: 2024††pagerange: What makes a cosmic filament?  the dynamical
origin and identity of filaments I. fundamentals in 2D–A

##  1 Introduction

The present study aims to infer and present a fundamental physical definition
for the identity of cosmic filaments. On the basis of the analytical fully
non-linear Caustic Skeleton model of the cosmic web (Feldbrugge et al., 2018),
the outline of its filamentary spine is inferred from the eigenvalues of the
primordial tidal and deformation field. It yields an unequivocal, solid
physical definition of the principal constituent of the cosmic web, the
pervasive web-like network that constitutes the system of cosmic arteries
along which matter is transported from low-density voids toward massive
compact nodes.

On scales of a few to hundreds of Megaparsec dark matter, gas, and galaxies
have organized into a complex spatial network, defining an intricate spatial
pattern of elongated filaments and flattened walls surrounding near-empty void
regions and connecting at cluster nodes (Bond et al., 1996; van de Weygaert &
Bond, 2008). By far the most prominent feature in the cosmic web are the
filaments (Giovanelli et al., 1986; Colberg et al., 2005; Cautun et al., 2014;
Galárraga-Espinosa et al., 2024). The pervasive interconnecting skeleton of
filamentary ridges defines the spine of the cosmic web (Aragón-Calvo et al.,
2010a, b). More than 50% of matter, gas, and galaxies have aggregated in the
filamentary arteries of the cosmic web (Cautun et al., 2014; Ganeshaiah Veena,
2020).

The prominence of filaments is particularly striking in large-scale galaxy
surveys. Their existence and importance was already noticed by early galaxy
redshift surveys (Joeveer & Einasto, 1978; Giovanelli et al., 1986; de
Lapparent et al., 1986; Shectman et al., 1996). Subsequent generations of wide
redshift surveys have increasingly enforced the realization that cosmic
filaments form a fundamental building block of the universe’s infrastructure.
The 2dFGRS survey (Colless et al., 2003), the Sloan Digital Sky Survey SDSS
(Tegmark et al., 2004) and the 2MASS redshift survey (Huchra et al., 2012) all
emphasized this observation, culminating in the recently published astounding
vista offered by the DESI survey,

In recent years, numerous studies have addressed the influence of the
filamentary environment on the properties, formation, and evolution of
galaxies. This ranges from the direct impact on the generation of rotation and
spin orientation of halos and galaxies (Porciani et al., 2002; Aragón-Calvo et
al., 2007b; Hahn et al., 2007; Aragón Calvo, 2007; Schäfer, 2009; Hahn et al.,
2009, 2010; Tempel & Libeskind, 2013; Codis et al., 2012, 2015, 2018b;
Ganeshaiah Veena et al., 2018, 2019; López et al., 2021; López, 2024) to the
question in how far the mass, stellar content and star formation activity are
influenced (Cautun et al., 2014; Kraljic et al., 2018; Kuutma et al., 2017;
Laigle et al., 2018; Hellwing et al., 2021; Parente et al., 2023). Also, while
representing the dominant morphological feature of the spatial organization of
matter in the Universe, they function as the main transport channels on
Megaparsec scales (Courtois et al., 2013; Pomarède et al., 2017; Tully et al.,
2023). As mass streams out of lower-density regions to higher-density areas,
filaments funnel matter, gas, and galaxies toward the compact high-density
nodes of the cosmic web, the location of massive clusters. Within this
context, it is of considerable importance that filaments also dominate the
gravitational and tidal force field on Megaparsec scales (Kugul, 2020), the
same large-scale force field that even induces the rotation of filaments (Wang
et al., 2021; Xia et al., 2021). Moreover, various subtle gas dynamical
processes involved in the anisotropic flow into and along filaments may have a
decisive influence on the star formation activity and history of galaxies
residing in or migrating towards them (Aragon Calvo et al., 2019; Lu et al.,
2024).

The filamentary network covers a wide spectrum of scales. In terms of mass
density, we find filaments over at least five orders of magnitude (Cautun et
al., 2014). Most outstanding are the heavy and thick filamentary extensions of
the cluster nodes, the outer edges of the most prominent arteries in the
cosmic web. In the local Universe, the Pisces-Perseus chain is a beautiful
specimen of such features (Giovanelli et al., 1986). A representative sample
of filaments, identified in the SDSS survey, is the filament catalogue of
Tempel et al. (2014). Most filaments are more moderate-density arms, often
branching out of the major filaments. Towards low-density voidlike regions,
this culminates in a branching forest of tenuous tendrils (Cautun et al.,
2014; Alpaslan et al., 2014; Shivashankar et al., 2016; Jaber et al., 2024;
Aragon-Calvo, 2024). This multiscale pattern is a direct manifestation of the
hierarchical fashion in which structure has assembled itself over the
evolution of the Universe, in which structure gradually builds up through the
merging of smaller-scale features that emerged earlier. By implication, at
earlier cosmic eons we find a similar web-like arrangement on smaller scales
as we find currently on Megaparsec scales. In the context of the formation of
galaxies, we see that galaxies and galaxy halos formed at the node of an
intricate network of filamentary extensions, with substantial repercussions
for the outcome and process of galaxy formation (see e.g. Springel et al.,
2008; Bi et al., 2024). On a more general note, we may observe that the
complex multiscale cosmic network resembles that of other intricate ridgelike
networks in nature, not unlike that of nerve and artery systems (Neyrinck et
al., 2018; Vazza & Feletti, 2020).

###  1.1 Filament identification and classification

While the cosmic web represents the largest structure that emerged in the
Universe, its complex multiscale, morphological, and topological properties
have formed a challenge for its classification and analysis. Furthering our
insight into the role of filaments in the formation of structure (Icke, 1973;
Eisenstein et al., 1997; Cautun et al., 2014; Kugul, 2020; Wang et al., 2024),
their influence on the formation and evolution of galaxies, and in the
sensitivity of their structure and dynamics to global cosmological factors
(see e.g. Bonnaire et al., 2022) makes it necessary to establish an objective
physical definition of what constitutes a filament. Hence, understanding out
of which primordial configurations in the initial Gaussian random field these
structures emerged, and how different global cosmological conditions –
including that of the nature of dark matter, dark energy, and the mass of
neutrinos – impact their evolution and structure is of key importance. Yet, a
physically motivated, unequivocal understanding and agreement on what defines
a filament, or any of the other structural components of the cosmic web, has
as yet still to be established.

To this adds the complication that we do not yet have a complete picture of
how different galaxy populations may yield a biased view of the spatial
structures and connections in the cosmic matter distribution. It may imply
that filaments and voids detected in galaxy surveys may not always be
representative for those in the underlying dark matter distribution. A few
recent studies do indicate subtle yet significant structural and topological
differences (Zakharova et al., 2023; Bermejo et al., 2024).

Driven by these considerations, over the past decades, a range of methods and
formalisms have been forwarded for the detection and classification of
filaments. Libeskind et al. (2018) provide a review and comparison of more
than a dozen formalisms. One may distinguish at least five classes of filament
finders. Geometric filament finders are usually based on the Hessian of the
density or gravitational potential at each location. It includes the Twerb and
Vweb formalisms (Hahn et al., 2007; Forero-Romero et al., 2009; Hoffman et
al., 2012). The most sophisticated ones explicitly take into account the
multiscale nature of the mass distribution, of which MMF/Nexus is a particular
example (Aragón-Calvo et al., 2007a; Cautun et al., 2013). Topological methods
address the connections between structural singularities in the mass
distribution. They are amongst the most widely used formalisms in the studies
of the cosmic web, in particular, that of Disperse (Sousbie, 2011; Sousbie et
al., 2011). Other representatives are Spineweb (Aragón-Calvo et al., 2010b)
and Felix (Shivashankar et al., 2016). Possible biases that may affect the
various means of filament detection have been investigated in a recent study
by Dhawalikar & Paranjape (2024).

In addition to the geometric and topological formalisms, several alternative
methods have played a substantial role in the study of the cosmic web. Bisous
is a well-known stochastic formalism, involving Bayesian exploration based on
stochastic geometric modeling of filaments (Tempel et al., 2014). It forms the
basis for the widely used filament catalogue extracted from the SDSS survey
(Tempel et al., 2014). More recent developments often incorporate machine
learning codes (see e.g. Canducci et al., 2022). Perhaps the oldest
representatives for a systematic analysis of filamentary patterns are graph-
based methods. The Minimal Spanning Tree (MST) is a prime example and is
figuring prominently in the cosmic web analysis of the GAMA survey (Alpaslan
et al., 2014). There is even a class of cosmic web identifiers that exploit
the resemblance of the cosmic web to biological branching networks. The Monte
Carlo Physarum Machine, inspired by the growth of Physarum polycephalum ’slime
mold’, has been successfully applied to the structural analysis of the cosmic
web in both simulations and observations (Elek et al., 2021, 2022; Wilde et
al., 2023).

Possibly the most profound techniques for classification of the cosmic web are
those emanating from the analysis of the 6D phase-space structure of the
cosmic mass distribution (Shandarin, 2011; Shandarin et al., 2012; Neyrinck,
2012; Abel et al., 2012). Restricted to situations in which the initial
conditions are known, they yield the identification of the matter streams
constituting the migration of mass in the buildup of structure. It allows the
definition of an objective physical criterion for what constitutes the various
structural elements of the cosmic web. The present study pursues the ideal of
understanding the cosmic web in terms of the gravitational physics underlying
its formation and development by elaborating on the phase-space evolution of
the cosmic matter distribution and inferring the implied emergence, branching,
and connectivity of the cosmic web and its structural components.

###  1.2 Tides, Deformation and the Caustic Skeleton

The cosmic web represents a key phase in the dynamical buildup of structure in
the Universe. It emerges when the original long phase of linear evolution of
the primordial density and velocity field turns into a more advanced non-
linear stage involving contraction and collapse of mass inhomogeneities. As
such, it marks the transition from the primordial (Gaussian) random field to
highly non-linear structures that have fully collapsed into halos and
galaxies. The formation and evolution of the characteristic anisotropic
structures, i.e., the filaments and walls, are the result of anisotropic
deformations. It is the anisotropy of the force field and the resulting
deformation of the matter distribution which are at the heart of the emergence
of the web-like structure of the mildly non-linear mass distribution, which
has been recognized and accurately described in the mildly non-linear stage by
the Zel’dovich formalism (Zel’dovich, 1970). Hence, by implication, the
gravitational tidal force field, induced by the inhomogeneous mass
distribution, is the key to modeling and understanding the structure of the
cosmic web.

The seminal role of gravitational tidal force fields in shaping the
anisotropic wall-like and filamentary structures in the cosmic web has been
recognized for many years (Zel’dovich, 1970; van de Weygaert & Bertschinger,
1996; Bond et al., 1996; Hahn et al., 2007; van de Weygaert & Bond, 2008;
Paranjape et al., 2018; Paranjape, 2021). In recent years the dominant
position of the tidal force field in laying out the spatial structure of the
cosmic web has led to the formulation of an analytical model for the fully
non-linear evolution of the cosmic web, the caustic skeleton (Arnol’d, 1982;
Hidding et al., 2014; Feldbrugge et al., 2018; Feldbrugge & van de Weygaert,
2023; Feldbrugge et al., 2023). It is based on the realization that the
evolution of the cosmic web can be understood in detail in terms of the
singularities and caustics that are arising in the matter distribution as a
result of the structure of the corresponding multi-stream flow field. The
geometric pattern of the large-scale cosmic matter distribution is to be
understood in terms of the folding of the dark matter sheet in phase-space
(Shandarin, 2011; Shandarin et al., 2012; Neyrinck, 2012; Abel et al., 2012).
A hierarchy of non-linear structures emerges after shell crossing in the
corresponding multi-stream regions. The mathematical framework of catastrophe
theory (Thom, 1975; Zeeman, 1977; Poston & Stewart, 1996; Arnold, 1984)
defines and classifies the stable folding configurations, whose geometrically
assemble into the connecting framework of the caustic skeleton (Arnol’d, 1982;
Hidding et al., 2014; Feldbrugge et al., 2018). The resulting phase-space
based caustic skeleton description of the evolving web-like pattern in the
cosmic matter distribution is codified in terms of a complete set of caustic
conditions. The study by Feldbrugge et al. (2018) demonstrated that a full
understanding of the cosmic web structure is obtained through the spatial
characteristics of the eigenvalue and eigenvector fields of the cosmic tidal
force field. Underlining this is the realization that the embryonic outline of
the cosmic web, in particular its filamentary network, can already be seen in
the primordial tidal eigenvalue field (Feldbrugge et al., 2023; Wilding, 2022)
(see figures 1 and 2).

In the current study, we focus on the nature of filaments within the caustic
skeleton framework. Of key importance is the observation that it implies that
the proper dynamical identity of filaments is dictated by the distribution and
connectivity of the relevant singularities in the tidal eigenvalue field.
Instrumental is the realization that this concerns the complex phase-space
singularities that are inferred from the Lagrangian description of the
evolving mass distribution (Feldbrugge et al., 2018; Feldbrugge & van de
Weygaert, 2023).

As such it differs fundamentally from studies – and identification methods –
that base themselves on the visual impression of filaments as bridges between
maxima in the mass distribution (see e.g Colberg et al., 2005), identifying
filaments with the saddle points and corresponding integral lines in the
corresponding cosmic density field (Pogosyan et al., 2009b, a; Sousbie, 2011;
Sousbie et al., 2011; Codis et al., 2018a). However, the implicit
identification of morphological and structural features on the basis of
density contour levels is probably based on an unwarranted assumption (see
e.g. Cautun et al., 2014, figure 13): filaments span no less than five orders
of magnitude in density, ranging from the massive prominent filamentary
extensions of clusters (Kuchner et al., 2020, 2021) to the tenuous and small
intravoid tendrils (Park & Lee, 2009; Shivashankar et al., 2016; Jaber et al.,
2024). This renders any density threshold-based identification
unrepresentative. This may also concern more sophisticated topological
identifications on the basis of density field filtrations (Sousbie et al.,
2011), which was demonstrated by the study of Shivashankar et al. (2016) and
which may substantially affect the detection of filaments in the galaxy
distribution (see e.g. Kuchner et al., 2020, 2021).

Of even more fundamental importance is the issue of the nature of the
filamentary network, and the question of whether filaments are exclusively
Bridges or whether there is a substantial population of Branches that fork off
from the main arteries, gradually fading into lower density areas without
revealing any connection to a higher density tip. In particular, in low-
density void regions we tend to find a system of tenuous tendrils, emanating
from major filamentary arteries at the boundary of voids, which bear a
striking resemblance to a network of branches (see e.g. Rieder et al., 2013;
Cautun et al., 2014; Alpaslan et al., 2014; Shivashankar et al., 2016; Jaber
et al., 2024; Aragon-Calvo, 2024). This may actually be a direct manifestation
of the multiscale nature of the filamentary web and its hierarchical
evolution. The migration of mass towards high-density regions along the small-
scale filaments culminated in the aggregation of mass in surrounding large-
scale filaments and the evacuation of mass from the low-density tips of the
small-scale filaments, turning these into tendrils. The direct implication is
that we do need to take care that any identification of filaments takes into
account the fact that filaments are not necessarily Bridges, but may also be
Branches.

Based on the caustic skeleton formalism, we seek to identify the objective and
dynamical identity of filaments in the cosmic matter distribution. In previous
studies, we found that filaments can be identified – in the three-dimensional
context – with two distinct singular features, the $A_{4}$ swallowtail
singularity, and the $D_{4}$ elliptic/hyperbolic umbilic. Visual inspection of
$N$-body simulations reveals that the latter are the high-density filamentary
extensions from the cluster nodes in the cosmic matter distribution. In the
present study, we elaborate on this by addressing the question of the precise
nature of the identity and outline of the filamentary network of the cosmic
web. The first aspect is to establish the formal criteria of the filamentary
spine of the cosmic web, on the basis of the related caustic singularities and
their connections. Subsequently, we study the gravitationally contracted
matter distribution around the filamentary spine, in Eulerian space, on the
basis of the recently developed non-linear constrained formalism (Feldbrugge
et al., 2023).

###  1.3 Filamentary explorations

The ultimate intention of the presented formalism is to establish a physically
based definition of filaments that allows us (1) to explore in a targeted
laboratory approach the properties of filaments as a function of different
primordial conditions, different spatial scales and at different cosmic
epochs, (2) to identify the rich cosmological information contained in the
intricate and complex filamentary network of the cosmic web and (3) to
investigate the formation, evolution and properties of dark halos and galaxies
as a function of filamentary environment.

With respect to the latter aspect, the present study provides a physically
based template for systematic studies of halo and galaxy formation as a
function of large scale environment. In this sense, it forms an an elaboration
and extension of earlier systematic studies, ranging from the study of disk
galaxy orientation in and around a zoom-in simulation of a cosmic filament
(Hahn et al., 2010), of (galaxy) halos forming in a large cosmic void (Rieder
et al., 2013) to explicit and systematic techniques involving the splicing and
tuning of large-scale and small-scale fluctuations in the primordial Gaussian
density field. Following up on a similar ideas employed in the peak-patch
formalism of Bond & Myers (1996), the principal example of this is the MIP
ensemble simulation formalism forwarded by Aragon-Calvo (2016) and its
elaborations (e.g. Cadiou et al., 2021).

In the present study, we limit ourselves to the two-dimensional situation, in
order to establish the key characteristics of the phase-space-based
filamentary network in the cosmic web. The considerably more complex three-
dimensional situation, involving a larger number of caustic singularities, is
the subject of the sequel to this study.

###  1.4 Outline

Towards this end, the present paper on the nature of filaments is organized as
follows. In section 2 we summarize the essential elements of the theory of the
caustic skeleton and its definition in terms of the deformation tensor
eigenvalue and eigenvector field. In section 3 we discuss the initial
conditions for structure formation and consider the caustic skeleton in scale
space to focus attention on structures at a specific length scale. In section
4 we identify the location of the ridge and center of filaments. Subsequently,
in section 5, we assess the statistical properties of the filament population
focusing on the filament formation time and number density. This leads to the
treatment in section 6 of the relation and connectivity of filaments within
the context of the cosmic web. In section 7, we compare the caustic definition
of cosmic filaments with different suggestions for the identity of filaments
using $N$-body simulations. Finally, in section 8 we outline the
generalization of the current two-dimensional analysis to the full three-
dimensional situation.

##  2 The caustic skeleton

The Caustic Skeleton is an analytical model that describes the emergence and
fully non-linear buildup of the cosmic web. The formalism is based on the
assessment of the evolving cosmic mass distribution in 6D phase-space. It
focuses on the folding of the dark matter sheet in phase-space and the
emergence and formation of multi-stream regions as the mass distribution
evolves from its pristine Gaussian configuration towards an ever more complex
non-linear arrangement of mass. For a detailed exposition of caustic skeleton
theory and its mathematical nuances see Arnol’d et al. (1982); Arnol’d (1982);
Arnold (1984); Hidding et al. (2014); Feldbrugge et al. (2018); Feldbrugge &
van de Weygaert (2023); Feldbrugge et al. (2023).

In the current study we consider a heuristic two-dimensional model of large-
scale structure formation in the context of the Caustic Skeleton model, in
order to elucidate the mathematical structure of filaments in the cosmic web.
We consider an Einstein-deSitter (EdS) cosmology. In an EdS Universe, the
growing mode of (linear) structure growth is proportional to the scale factor,

$b_{+}(t)\,\propto\,a(t)\,,$ (1)

so that it will suffice to parametrize structure formation in terms of the
growing mode, ranging from $b_{+}(0)=0$ to $b_{+}(t_{0})=1$ at the current
epoch.

We start our investigation by discussing the theoretical concepts that
underpin the caustic skeleton model. Fundamental is the description of the
cosmic mass distribution in a Lagrangian context, formulated in terms of
Lagrangian fluid dynamics. The folding of phase-space sheets leads to families
of Lagrangian singularities, whose identity is fully determined by the
eigenvalues and eigenvectors of the deformation tensor (Feldbrugge et al.,
2018, 2023). The expressions for the specific two-dimensional context of the
present study are specified in a separate subsection. Within the context of
cosmic structure formation, the first-order Lagrangian scheme is of
fundamental importance, the Zel’dovich approximation (Zel’dovich, 1970). It
plays a central role in the work described in this study and allows a fully
analytical treatment of the Lagrangian singularities. The description of these
also provides the notation used in the present study.

###  2.1 Lagrangian fluid dynamics

The caustic skeleton describes the cosmic web in terms of the stable
singularities of the dark matter sheet in Lagrangian fluid dynamics.
Lagrangian fluid dynamics describes the evolution of the cold dark matter
fluid in terms of the Lagrangian map $\bm{x}_{t}:L\to E$, mapping initial
positions in Lagrangian space $L$ to final positions in Eulerian space $E$ at
time $t$, i.e., a mass element starting at the initial position $\bm{q}\in L$
moves in time $t$ to the position

$\displaystyle\bm{x}_{t}(\bm{q})=\bm{q}+\bm{s}_{t}(\bm{q})\,,$ (2)

with he displacement map $\bm{s}_{t}$ describing its relative motion. The
associated density field is a derived property of the fluid following from the
conservation of mass of the mass elements,

$\displaystyle\rho(\bm{x})$
$\displaystyle=\sum_{\bm{q}\in\bm{x}_{t}^{-1}(\bm{x})}\frac{\bar{\rho}}{|\det(\nabla\bm{x}_{t}(\bm{q}))|}$
$\displaystyle=\sum_{\bm{q}\in\bm{x}_{t}^{-1}(\bm{x})}\frac{\bar{\rho}}{|\det(I+\nabla\bm{s}_{t}(\bm{q}))|}\,,$
(3)

where $\det(\nabla\bm{x}_{t}(\bm{q}))$ is the relative signed volume of a mass
element at time $t$, with the mean density $\bar{\rho}$, and the preimage

$\displaystyle\bm{x}_{t}^{-1}(\bm{x})=\\{\bm{q}\in
L\,|\,\bm{x}_{t}(\bm{q})=\bm{x}\\}\,,$ (4)

consisting of the Lagrangian points reaching $\bm{x}$ in time $t$. The number
of streams $n$ (defined as the number of elements in the preimage
$\bm{x}_{t}^{-1}(\bm{x})$) partitions Eulerian space into $n$-stream
regions111The number of streams is an odd natural number by the topology of
folding.. At the boundaries of a multi-stream region, the deformation tensor
$\nabla\bm{x}_{t}$ is singular and the density spikes to infinity. This spike
marks the onset of the non-linear gravitational evolution of the mass element
and is known as the phenomenon of shell-crossing. More formally, this is known
as the fold caustic. Subsequently, we see the emergence of a family of ever
more complex and intricate Lagrangian singularities. Their definition,
identification with structural features of the cosmic web and mathematical
expressions form the core of the Caustic Skeleton model (Feldbrugge et al.,
2018).

###  2.2 Lagrangian singularities

The formal definition and expressions for the Lagrangian singularities can be
fully specified in terms of the eigenvalue and eigenvector fields of the
deformation tensor $\nabla\bm{s}_{t}$. The eigenvalue and eigenvector fields
of the deformation tensor $\nabla\bm{s}_{t}$ are defined by the eigen equation

$\displaystyle[\nabla\bm{s}_{t}(\bm{q})]\bm{v}_{t,i}(\bm{q})=\mu_{t,i}(\bm{q})\bm{v}_{t,i}(\bm{q})\,,$
(5)

with the ordering $\mu_{t,1}(\bm{q})\leq\mu_{t,2}(\bm{q})$. On the basis of
the eigenvalues and eigenvectors we can identify two families of
singularities, the $A$ singularities that depend on only one eigenvalue, and
the umbilic D singularities that are defined by two eigenvalues. For the two-
dimensional situation, we deal with two eigenvalues $\mu_{1}$ and $\mu_{2}$,
and identify three A singularities – the fold caustic $A_{2}$, the cusp
caustic $A_{3}$ and the swallowtail caustic $A_{4}$ – and one $D$ singularity,
the umbilic $D_{4}$ singularity.

The fold caustic $A_{2}$ in Lagrangian space corresponding to the eigenvalue
field $\mu_{t,i}$ (with $i=1,2$) consists of the manifold

$\displaystyle A_{2}^{i}=\\{\bm{q}\in L\,|\,1+\mu_{t_{0},i}(\bm{q})=0\\}\,,$
(6)

with $t_{0}$ the current time. The density spikes at the fold manifold.
Conversely, the fold manifold bounds the mass elements having undergone shell-
crossing along the $\bm{v}_{t_{0},i}$ direction.

Interestingly, the structural complexity and classification of caustics do not
stop at the fold caustic. Caustic theory provides a refined description of the
intricate folding of the dark matter sheet in phase-space (Arnol’d, 1982;
Zeeman, 1977; Thom, 1975). The next-order caustic in this hierarchy of
singularity is the cusp caustic $A_{3}$,

$\displaystyle A_{3}^{i}=\\{\bm{q}\in L\,|\,$ $\displaystyle
1+\mu_{t,i}(\bm{q})=0\text{ and
}\bm{v}_{t,i}(\bm{q})\cdot\nabla\mu_{t,i}(\bm{q})=0$ $\displaystyle\text{ for
some }t\leq t_{0}\\}\,.$ (7)

In addition to $\nabla\bm{x}_{t}$ being singular due to the condition
$1+\mu_{t,i}=0$, the derivative of the eigenvalue field $\mu_{t,i}$ in the
direction of the corresponding eigenvector $\bm{v}_{i,t}$ vanishes. The cusp
manifold marks the mass elements that form a kink in the boundary between the
different $n$-stream regions. Within the context of the formation of the
cosmic web, (in the 2D context) it corresponds to filaments.

In addition to the cusp caustic, a two-dimensional Lagrangian fluid
generically also contains the higher order $A_{4}$ swallowtail caustics. The
swallowtail caustic consists of the points for which

$\displaystyle A_{4}^{i}=\\{\bm{q}\in L\,|\,$ $\displaystyle
1+\mu_{t,i}(\bm{q})=0\,,\bm{v}_{t,i}(\bm{q})\cdot\nabla\mu_{t,i}(\bm{q})=0\,,$
$\displaystyle\text{ and
}\bm{v}_{t,i}(\bm{q})\cdot\nabla(\bm{v}_{t,i}(\bm{q})\cdot\nabla\mu_{t,i}(\bm{q}))=0$
$\displaystyle\text{ for some }t\leq t_{0}\\}\,.$ (8)

(a) $\lambda_{1}$

(b) $\lambda_{2}$

(c) Zel’dovich approximation

(d) $N$-body simulation

Figure 1: The caustic skeleton in Lagrangian (upper panels) and Eulerian space
(lower panels). The upper left and upper right panels show the Lagrangian
skeleton in relation to the first ($\lambda_{1}$) and second eigenvalue field
($\lambda_{2}$). The lower left and right panel show the skeleton in Eulerian
space, based on the Zel’dovich approximation and an $N$-body simulation.

Finally, the umbilic caustic $D_{4}$ is defined by the points for which the
two eigenvalue fields coincide

$\displaystyle D_{4}=\\{\bm{q}\in L\,|\,$ $\displaystyle
1+\mu_{t,1}(\bm{q})=0\text{ and }1+\mu_{t,2}(\bm{q})=0\,.$
$\displaystyle\text{ for some }t\leq t_{0}\\}\,,$ (9)

Around the umbilic caustic, the density spikes more strongly. This follows
from the fact that the matrix $\nabla\bm{x}_{t}$ is singular along two
directions.

Singularity | Singularity | Feature in the
---|---|---
class | name | 2D cosmic web
$A_{2}$ | fold | collapsed region
$A_{3}$ | cusp | filament
$A_{4}$ | swallowtail | cluster or knot
$D_{4}$ | hyperbolic/elliptic | cluster or knot
Table 1: The identification of caustic singularities with cosmic web
components, for (for the $2$-dimensional situation)

###  2.3 Caustic skeleton & Cosmic Web

The embryonic caustic skeleton in Lagrangian space can be mapped onto Eulerian
space by means of the Lagrangian map. The resulting varieties
$\bm{x}_{t_{0}}(A_{2}^{i}),$ $\bm{x}_{t_{0}}(A_{3}^{i}),$
$\bm{x}_{t_{0}}(A_{4}^{i}),$ and $\bm{x}_{t_{0}}(D_{4})$ delineate the various
structural components of the cosmic web. Table 1 lists the identification of
the caustic singularities with the well-known structural features of the
cosmic web. In other words, the Lagrangian maps capture the geometry of the
cosmic web, and their evolution their formation history.

###  2.4 The Zel’dovich approximation

The focus of the present study is on the exploration and presentation of the
mathematical structure of filaments in the cosmic web. To this end, we center
the analysis on a fully analytical heuristic model of cosmic structure
formation. To optimize transparency and insight into the implied structure and
dynamics of filaments, we here explore the two-dimensional situation and
relegate the full three-dimensional situation to the follow-up paper.

The Zel’dovich approximation provides a fully analytical description of the
displacement of mass elements and hence entails a fully analytical expression
for the deformation tensor. Within the context of the caustic skeleton
formalism, it therefore allows a fully analytical treatment. The Zel’dovich
approximation describes the evolution of mass elements in terms of their
ballistic motion, in which time evolution and spatial structure are separated.
The mass elements follow linear trajectories determined by the primordial
gravitational potential. This follows from the first-order approximation of
the displacement $\bm{s}_{t}(\bm{q})$ of the mass element with Lagrangian
coordinate $\bm{q}$, yielding

$\displaystyle\bm{s}_{t}(\bm{q})=-b_{+}(t)\nabla\Psi(\bm{q})\,.$ (10)

The time evolution is represented by the growing mode $b_{+}(t)$, while the
spatial structure is encapsulated by the displacement potential $\Psi$,
expressed in terms of the current total energy density and the linearly
extrapolated gravitational potential $\phi_{0}$,

$\displaystyle\Psi(\bm{q})\,=\,\frac{2}{3\Omega_{0}H_{0}^{2}}\phi_{0}(\bm{q})\,,$
(11)

contains the geometric information of the evolution of the fluid.

The Zel’dovich approximation accurately describes single-stream regions from
early times up to surprisingly advanced phases of quasi-linear structure
formation (Zel’dovich, 1970; Shandarin & Zel’dovich, 1989; Hidding et al.,
2014). For our purpose, it represents an ideal context, given the focus on the
first shell-crossing events in the caustic skeleton formalism. It fails at
late times in multi-stream regions when the gravitational backreaction of the
mass elements becomes important. This follows directly from the observation
that the Zel’dovich approximation describes a ballistic motion, and hence
unable to capture gravitational interactions at later stages of structure
formation. In particular of relevance in the present context of phase-space
evolution, is that the Zel’dovich approximation does not capture the
turnaround typically occurring during phase mixing.

For conciseness, we define the eigenvalues $\lambda_{i}(\bm{q})$ of the
Hessian of the displacement potential $\Psi(\bm{q})$

$\displaystyle[\mathcal{H}\Psi]\bm{v}_{i}=\lambda_{i}\bm{v}_{i}\,,$ (12)

with the ordering $\lambda_{1}(\bm{q})\geq\lambda_{2}(\bm{q})$. The
corresponding (non-linear) density field in the Zel’dovich approximation
assumes the form (see eq. (3)),

$\displaystyle\rho(\bm{x})=\sum_{\bm{q}\in\bm{x}_{t}^{-1}(\bm{x})}\frac{\bar{\rho}}{|1-b_{+}(t)\lambda_{1}(\bm{q})||1-b_{+}(t)\lambda_{2}(\bm{q})|}\,.$
(13)

The caustic conditions, presented in section 2.2, in the Zel’dovich
approximation assume the following form for the caustic manifolds
corresponding to the first eigenvalue field (in which we have dropped the
eigenvalue subscript),

$\displaystyle A_{2}=\\{\bm{q}\,|\,\lambda_{1}(\bm{q})=1/b_{+}(t_{0})\\}\,,$
(14) $\displaystyle A_{3}=\\{\bm{q}\,|\,\lambda_{1}(\bm{q})\geq
1/b_{+}(t_{0})\,,\bm{v}_{1}(\bm{q})\cdot\nabla\lambda_{1}(\bm{q})=0\\}\,,$
(15) $\displaystyle A_{4}$
$\displaystyle=\\{\bm{q}\,|\,\lambda_{1}(\bm{q})\geq
1/b_{+}(t_{0})\,,\bm{v}_{1}(\bm{q})\cdot\nabla\lambda_{1}(\bm{q})=0\,,$
$\displaystyle\phantom{=\\{\bm{q}\,|\,\,}\bm{v}_{1}(\bm{q})\cdot\nabla(\bm{v}_{1}(\bm{q})\cdot\nabla\lambda_{1}(\bm{q}))=0\\}\,,$
(16) $\displaystyle
D_{4}=\\{\bm{q}\,|\,\lambda_{1}(\bm{q})=\lambda_{2}(\bm{q})\geq
1/b_{+}(t_{0})\\}\,.$ (17)

An illustration of a realization of the caustic skeleton of the two-
dimensional cosmic web is shown in figure 1. The two upper panels show the
outline of the caustic skeleton in Lagrangian space. The two lower panels show
the skeleton in Eulerian space, both for the Zel’dovich approximation as well
as for the fully non-linear outline inferred from the corresponding $N$-body
simulation. Each of the various caustic features is indicated using different
colours, red for the fold caustics $A_{2}$, blue for cusp caustics $A_{3}$
that delineate the filamentary spine of the cosmic web.

##  3 the Multiscale Skeleton

The cosmic web has evolved gravitationally from a field of primordial density,
gravity, and potential perturbations. In the present study, we assume the
primordial density fluctuations to be a realization of a statistically
homogeneous and isotropic stationary Gaussian random field.

As a result, the emerging cosmic web is an intricate multiscale structure,
with structure and features over a vast range of scales. In order to focus
attention on features around a specific length scale, we invoke the concept of
the Scale Space of the Caustic Skeleton.

###  3.1 Primordial Density and Potential field

A Gaussian random field is fully specified by its second-order moment. The
power spectrum fully quantifies the statistical properties of the Fourier
modes $\hat{\delta}(\bm{k})$ of the density field,

$\displaystyle\hat{\delta}(\bm{k})=\int\delta(\bm{q})e^{i\bm{k}\cdot\bm{x}}\mathrm{d}\bm{q}\,,$
(18)

in terms of their covariance,

$\displaystyle\langle\hat{\delta}(\bm{k}_{1})\hat{\delta}^{*}(\bm{k}_{2})\rangle=(2\pi)^{2}\delta_{D}^{(2)}(\bm{k}_{1}-\bm{k}_{2})P_{\delta}(\|\bm{k}_{1}\|)\,,$
(19)

in which $\delta_{D}^{(2)}$ is the two-dimensional Dirac delta function.

For the purpose of the present study, we use a heuristic power spectrum that
consists of a scale-free power spectrum in combination with a Gaussian cutoff
at scale $R_{s}$ at high frequencies,

$\displaystyle P_{\delta}(k)=\mathcal{N_{\delta}}k^{n}e^{-R_{s}^{2}k^{2}}\,.$
(20)

The scale-free part of the spectrum is characterized by the spectral index
$n$, the cutoff part by the Gaussian smoothing scale $R_{s}$, while
$\mathcal{N_{\delta}}$ is the normalization factor. Notice that the the
Gaussian spectral cutoff is meant to represent an intrinsic physical spectral
aspect, not unlike that seen in warm dark matter or neutrino-dominated
cosmologies and to some extent not unlike the more gradual high-frequency
downturn seen in ($\Lambda$)CDM cosmologies.

By implication, also the displacement potential $\Psi$ of the Zel’dovich
approximation is a Gaussian random field. Given the Poisson equation,

$\nabla^{2}\Psi\propto\delta\,,$ (21)

it follows that the power spectrum

$\displaystyle P_{\Psi}(k)=\mathcal{N}_{\Psi}k^{n-4}e^{-R_{s}^{2}k^{2}}\,.$
(22)

In the present study, we take the primordial gravitational potential and
displacement potential to follow the Harrison-Zel’dovich scale-invariant
spectrum on large scales (Harrison, 1970; Zel’dovich, 1972; Mukhanov &
Chibisov, 1981; Guth & Pi, 1982; Starobinsky, 1982; Bardeen et al., 1983), for
which the spectral index $n=2$. On small scales, we invoke a cutoff filter
scale $R_{s}$ to smooth the fields. In our calculations, we normalize the
power spectrum on the basis of the amplitude of the density perturbations,

$\langle\delta^{2}(\bm{q})\rangle=\alpha^{2}\,,$ (23)

which sets the normalization $\mathcal{N_{\delta}}$ to

$\displaystyle\mathcal{N_{\delta}}=\frac{\alpha^{2}4\pi
R_{s}^{n+2}}{\Gamma\left(1+\frac{n}{2}\right)}\,,$ (24)

with the gamma function $\Gamma$. For our calculations, we set the density
field amplitude to $\alpha^{2}=10$. While this sets the time scale of the
simulations, the employed cutoff length scale $R_{s}=1$ sets the spatial scale
of the simulations.

(a) $\Psi$

(b) $\delta$

(c) $\lambda_{1}$

(d) $\lambda_{2}$

Figure 2: The displacement potential (upper left), with the associated
primordial density perturbation (upper right), and associated eigenvalue
fields (lower left and lower right). The primordial eigenvalue fields display
the progenitors of the late-time cosmic web.

###  3.2 Scale Space Cosmic Web

While the filamentary structure of the multiscale cosmic web is that of a
complex and intricate aggregate of elongated bridges and branches, in the
present study we chose to zoom in on features of a specific length scale
$\sigma$, and the implied corresponding mass scale.

To direct attention to features on a scale $\sigma$, we chose to follow the
strategy of studying the formation of the caustic skeleton from
correspondingly smoothed primordial conditions. A similar strategy has been
followed in various studies on the multiscale evolution and structure of the
cosmic web by (e.g. Aragón-Calvo et al., 2010a; Aragon-Calvo & Szalay, 2013;
Jaber et al., 2024; Aragon-Calvo, 2024). While the non-linear gravitational
evolution of cosmic structure generates its structure over a spectrum of
scales, the cosmic web finds itself in the transition from linear to non-
linear scales and as such retains a considerable level of memory of its
primordial spatial structure. This is in particular so for the spatial pattern
of the tidal/deformation field that outlines the cosmic web (see Feldbrugge et
al., 2023).

To accomplish this Scale Space approach, we smooth the displacement field on a
Gaussian scale $\sigma$. In a sense, this linearly takes into account the
nonlocal influences on the emergence of the caustic skeleton: the skeleton is
not only sensitive to the displacement at one single point but emerges through
the influence of the neighborhood on the mass element.

To this end, we convolve the displacement potential $\Psi$,

$\displaystyle\Psi_{\sigma}(\bm{q})=\int\Psi(\bm{q}-\bm{q}^{\prime})W_{\sigma}(\bm{q}^{\prime})\,\mathrm{d}\bm{q}^{\prime}\,,$
(25)

with a Gaussian kernel222An alternative choice frequently used in late-time
cosmology is to invoke a spherical (or rather circular in the two-dimensional
context) tophat kernel  $\displaystyle
W_{\sigma}(\bm{q})=\frac{\Theta_{H}(\|\bm{q}\|-\sigma)}{\pi\sigma^{2}}\,,$
(26) with the Heaviside step function $\Theta_{H}$. While we do not expect the
results to significantly alter for similar alternative smoothing kernels, here
we only consider Gaussian smoothing. on scale $\sigma$,

$\displaystyle
W_{\sigma}(\bm{q})=\frac{1}{2\pi\sigma^{2}}e^{-\frac{\bm{q}^{2}}{2\sigma^{2}}}\,.$
(27)

The resulting scale-space power spectrum of the smoothed displacement
potential $\Psi_{\sigma}$ assumes the form

$\displaystyle
P_{\Psi_{\sigma}}(k)=\mathcal{N}k^{n-4}e^{-(R_{s}^{2}+\sigma^{2})k^{2}}\,.$
(28)

We emphasize the difference between the smoothing scales $R_{s}$ and $\sigma$:
the spectral smoothing scale $R_{s}$ is a physically intrinsic aspect of the
fluctuation spectrum, while $\sigma$ is a user-imposed smoothing scale meant
to study the multiscale structure of the cosmic web.

Within the context of the scale dependence of the filamentary web, and our
corresponding sampling schemes, an important role is played by the generalized
moments,

$\displaystyle\sigma_{i}^{2}$
$\displaystyle=\frac{1}{(2\pi)^{2}}\int\|\bm{k}\|^{2i}P_{\Psi_{\sigma}}(\bm{k})\mathrm{d}\bm{k}$
(29)
$\displaystyle=\frac{1}{2\pi}\int_{0}^{\infty}k^{2i+1}P_{\Psi_{\sigma}}(k)\mathrm{d}k$
$\displaystyle=\frac{\alpha^{2}R_{s}^{n+2}}{\left(R_{s}^{2}+\sigma^{2}\right)^{i+\frac{n}{2}-1}}\frac{\Gamma\left(i+\frac{n}{2}-1\right)}{\Gamma\left(\frac{n}{2}+1\right)}\,.$

These represent convenient parameterizations of the statistical properties of
Gaussian random fields,

As described below, these generalized moments play an important role in the
properties of filaments and our sampling schemes. For the particular case of a
2D Harrison-Zel’dovich primordial spectrum, with $n=2$, we have

$\displaystyle\sigma_{2}=\frac{\alpha R_{s}^{2}}{R_{s}^{2}+\sigma^{2}}\,.$
(30)

###  3.3 Scaling considerations: space vs. time

A scaling in real space $\bm{x}\mapsto L\bm{x}$ for some $L>0$ leads to a
scaling $\bm{k}\mapsto\bm{k}/L$ in Fourier space. Under this transformation,
the power spectrum scales as

$\displaystyle
P_{\Psi_{\sigma}}(k/L)=L^{4-n}\mathcal{N}k^{n-4}e^{-(R_{s}/L)^{2}k^{2}}\,.$
(31)

Rescaling in position space $\bm{x}\mapsto L\bm{x}$ thus corresponds to a
Gaussian smoothing with the length scale

$\displaystyle\sigma=\frac{\sqrt{1-L^{2}}}{L}R_{s}\,,$ (32)

and a scaling of the amplitude of the Gaussian random field

$\displaystyle\alpha\mapsto L^{2-\frac{n}{2}}\alpha\,.$ (33)

The latter can be interpreted as a rescaling of the growing mode,

$\displaystyle b_{+}(t)\mapsto L^{2-\frac{n}{2}}b_{+}(t)\,.$ (34)

For these simple power law power spectra, it thus suffices to study the
properties of the caustic skeleton at a fixed smoothing scale $\sigma$ as a
function of time $t$, or conversely for a fixed time $t$ as a function of the
scale $\sigma$. For more general power spectra this degeneracy is generically
broken.

##  4 Caustic Filaments

Within the context of the caustic skeleton formalism, the Zel’dovich
approximation enables us to infer specific analytical expressions for
filaments and relate their configuration directly to the Gaussian primordial
initial conditions. Given the fact that the cosmic web is a non-linear complex
geometric structure, this is a remarkable observation. It emanates from the
realization that the spatial outline of the cosmic web may already be
recognized in the primordial Universe when assessing the primordial
tidal/deformation eigenvalue field (see fig. 2 and e.g. Feldbrugge et al.
(2023)).

As a result of gravitationally driven evolution of the primordial mass
distribution, the embryonic pattern of the cosmic web is shaped into a complex
multiscale weblike network. It includes major arteries suspended between
clusters at the nodes of the network (top panel, figure 3), arteries tapering
off while stretching into low density void regions (central panel, figure 3
and configurations of tendrils that branch out from more massive peers (lower
panel, figure 3).

Figure 3: Top: a sketch of a filament (blue) connecting two clusters (red).
Matter flows into the clusters along the filaments (black arrows). The
stretching centers of the filaments, where due to the shearing migration flow
along the filament mass elements get maximally stretched, are indicated by the
blue points. The left panel illustrates a filament between two clusters.
Centre: The right filament branching off via cusp caustics. Bottom: multiscale
filament network and filament branching.

###  4.1 Filament Caustic

Filaments of the two-dimensional cosmic web correspond to the $A_{3}$ cusp
manifold of the first eigenvalue field $\lambda_{1}$. In the Zel’dovich
approximation, in Lagrangian space, these assume the form

$\displaystyle A_{3}=\\{\bm{q}\,|\,$ $\displaystyle\lambda_{1}(\bm{q})\geq
1/b_{+}(t_{0})\,,\bm{v}_{1}(\bm{q})\cdot\nabla\lambda_{1}(\bm{q})=0\\}\,,$
(35)

with an orientation in Lagrangian space given by the normal

$\displaystyle\bm{n}=\nabla(\bm{v}_{1}(\bm{q})\cdot\nabla\lambda(\bm{q}))\,.$
(36)

The fact that the caustic skeleton model allows the determination of $A_{3}$
filaments solely based on the initial eigenvalue and eigenvector field
represents a remarkable simplification concerning a fully non-linear
gravitational assessment. Amongst others, it facilitates a statistical
assessment of filament properties, by its direct relation to the Gaussianity
of the primordial field (see section 5).

The caustic equation (35) identifies all Lagrangian points ${\bm{q}}$ that
belong to the $A_{3}$ singularity. Its spatial outline in the evolving cosmic
matter distribution is obtained through the Lagrangian map towards Eulerian
space, which translates each mass element $\bm{q}_{c}$ of the Lagrangian
singularity to the Eulerian position $\bm{x}_{t}(\bm{q}_{c})$ to which it has
migrated at time $t$. A few complications arise in the pure translation of the
complete $A_{3}$ caustic in Lagrangian space to Eulerian space. This makes it
necessary to distinguish between the main ridge of the resulting filament, the
trunk of the filament, and the parts of the caustic that form its connection
to the higher-density regions.

Firstly, as the trunks of the filaments stretch while the connecting regions
compress during the evolution of the cosmic web, a general point on the cusp
manifold in Lagrangian space is biased to end up in a cluster-like region. In
other words, a general point on the cusp manifold is not the most
representative point on this manifold of mass elements ending up in the
filamentary structures at the current epoch.

Secondly, in the practical implementation within the context of the Zel’dovich
approximation which we use for the analytical description in the present
study, we are not able to follow the development of phase-space and the
resulting caustic structure in the high-density contracting region333Note that
more accurate higher order descriptions or the use of $N$-body simulation
inferred deformation values would allow a more accurate caustic description in
the evolving phase-space in and around high-density regions.. In these
regions, mass elements are part of an ever more complex wrapping phase-space
structure (see e.g. Colombi, 2015), which the straight ballistic motion of
mass elements implied by the Zel’dovich approximation entirely fails to
capture.

The combination of both aspects outlined above makes it necessary to augment
the straightforward projection of the mass elements along the $A_{3}$ caustic
from Lagrangian to Eulerian space with an additional criterion on the
stretching of mass elements along the trunk of the filament. To this end, we
introduce a local criterion for the central stretching region of a filament,
which is outlined in detail below.

###  4.2 Filament Trunk & Stretching Center

To infer a complete criterion for the main ridge of filaments within the
context of caustic skeleton theory we invoke a straightforward dynamical
criterion concerning the effect of the large-scale tidal force field on the
mass elements in a filament and the induced migration flows along the
filament. To this end, we base ourselves on the realization that filaments are
formed and shaped by the large-scale tidal force field. These induce
characteristic shearing migration flows towards and along the emerging
filaments.

As a result of the tidally induced migration streams, mass elements in the
main ridge get stretched. By contrast, those at the higher density regions
experience contraction and find themselves in a region in which the phase-
space structure quickly evolves towards an ever more complex configuration. It
is the stretching part of the caustic as defining the cosmic filament which we
would recognize visually as such: the trunk of the filament. Following the
notion of the induced shearing flow along the filaments’s ridge, we see a
variation of flow velocity along the filament, increasing towards the high-
density nodes towards which mass is migrating (see fig. 3).

On the basis of this, we employ the natural dynamical definition for filaments
as the regions along which the matter distribution is being stretched along
its spine (see fig. 3 for a sketch of this configuration). We may note that
there are also regions along a filament that are not stretching, such as near,
around, and in the cluster nodes at the tips of a filament. Instead, as mass
elements approach such high-density nodes, subsequently, pass through them and
accrete on them, they may get compressed. Hence, we are principally only
interested in the main stretching part of the ridge a filament, the trunk of a
filament.

Following this definition, we may identify a characteristic dynamical center
of a filament, the stretching center of a filament – or, shortly, the center
of a filament – as the location along the filament where the flow corresponds
to a maximum stretching of the spatial mass distribution along the filament’s
ridge. Within our study, these filament centers are a key handle for allowing
a mathematical analysis of the distribution and statistical properties of
filaments. This is particularly true when invoking the Zel’dovich
approximation, which allows a fully analytical formalization in terms of
implied elementary algebraic conditions.

###  4.3 the Filament Stretching Center: Conditions

To infer the formal conditions to specify the location of the trunk and
stretching center of a filament, we invoke and combine several constraints on
the eigenvalue and eigenvector fields:

####  4.3.1 Two Eigenvalues

To develop a definition for the stretching center of a filament, we first note
that clusters – in 2D – are typically identified with the regions in
Lagrangian space that undergo shell-crossing along at least two directions. In
the Zel’dovich approximation, this implies that in a cluster region, both
eigenvalues $\lambda_{1}$ and $\lambda_{2}$ exceed $1/b_{+}(t_{0})$. On the
basis of this realization, we may already conclude that a full filament
criterion should involve both the first and second eigenvalue fields.

####  4.3.2 Tangent Space and Stretching Condition

The Taylor series of the Lagrangian map in the Lagrangian point $\bm{q}_{0}$,

$\displaystyle\bm{x}_{t}(\bm{q}_{1})-\bm{x}_{t}(\bm{q}_{0})=\nabla\bm{x}_{t}(\bm{q}_{0})(\bm{q}_{1}-\bm{q}_{0})+\dots$
(37)

demonstrates that a vector in the tangent space of Lagrangian space $\bm{T}\in
T_{\bm{q}}L$ in the point $\bm{q}\in L$ is mapped to the vector
$\nabla\bm{x}_{t}(\bm{q})\bm{T}$ in the tangent space in Eulerian space
$T_{\bm{x}_{t}(\bm{q})}E$.

In other words, as $\bm{q}_{1}\to\bm{q}_{0}$ the vector

$T=\bm{q}_{1}-\bm{q}_{0}$ (38)

is mapped to the vector $\nabla\bm{x}_{t}\bm{T}$ in Eulerian space.
Restricting the map to the points $\bm{q}_{c}$ located on the filamentary cusp
caustic $A_{3}$, the space $\bm{T}$ is identified with the unit tangent vector
of the cusp manifold. This is identified by the condition

$\displaystyle\bm{T}\cdot\bm{n}=0\,,$ (39)

with the normal of the cusp manifold $\bm{n}$. We will assume a unit tangent
vector, i.e., $\|\bm{T}\|=1.$ The key observation is that when the norm
$\|[\nabla\bm{x}_{t}(\bm{q}_{c})]\bm{T}\|$ exceeds unity, the filament
stretches along the tangent direction, i.e., the filament stretches along its
ridge. On the other hand, when the norm is smaller than one, the filamentary
region contracts at the time $t$.

On the basis of this relation we are therefore led to the formal definition of
the stretching center of a filament: the stretching center of the filament is
the location on the cusp manifold of maximum and enduring dilation.
Analytically, this point of maximum stretch can be inferred analytically
within the framework of the Zel’dovich approximation.

####  4.3.3 Tangent – Second Eigenvalue alignment

Figure 4: A comparison of the distribution of the angle between $\bm{T}$ and
$\bm{v}_{2}$ (in degrees) for general points on the cups line (black) and
points on the cusp line for which $\lambda_{2}<0$ (red) sampled in Lagrangian
space for $b_{+}(t_{c})=0.5$.

A complicating factor for the practical implementation of the stretching
conditions on the tangent space $\bm{T}$ is that they are algebraically and
computationally very hard to evaluate. Though geometrically elegant, the
maximum stretch condition on the tangent vector $\bm{T}$ yields a complicated
algebraic identity that cannot be straightforwardly implemented using linear
constrained theory.

Instead, we follow a simplifying approximation based on the observation of the
alignment between the second eigenvector field $\bm{v}_{2}(\bm{q}_{c})$ and
the tangent vector $\bm{T}$. As demonstrated in our earlier study (Feldbrugge
& van de Weygaert, 2023), on average these are aligned. This strong
correlation between the orientation of the tangent vector and the second
eigenvector is shown in the diagram in figure 4, for general points on the
cusp curve $A_{3}$ (black) and those with negative second eigenvalue field
$\lambda_{2}$ (red). It suggests us to follow the simplifying approximation of
a perfect alignment,

$\bm{T}=\bm{v}_{2}(\bm{q}_{c})\,.$ (40)

This translates into the identity

$\displaystyle\|\nabla\bm{x}_{t}(\bm{q}_{c})\bm{T}\|$
$\displaystyle=\|(I-b_{+}(t)\mathcal{H}\Psi(\bm{q}_{c}))\bm{v}_{2}(\bm{q}_{c})\|$
(41) $\displaystyle=|1-b_{+}(t)\lambda_{2}(\bm{q}_{c})|\,,$ (42)

with the unit two-by-two matrix $I$.

For positive eigenvalue $\lambda_{2}$, the line element contracts and
undergoes shell crossing at

$1-b_{+}(t)\lambda_{2}=0\,.$ (43)

Subsequently, with increasing growing mode the mass element expands.

On the other hand, when the eigenvalue $\lambda_{2}$ is negative, following
the Zel’dovich approximation the filament continues to stretch. We may
immediately observe that the location – or, rather, mass element – where the
filaments get maximally stretched corresponds to the mass element whose
eigenvalue $\lambda_{2}$ along the cusp manifold $A_{3}$ represents a local
minimum.

####  4.3.4 Minimum Second Eigenvalue

Following the relations and approximations on the filament stretching
conditions, we arrive at the set of conditions identifying the stretching
center of filaments. They are identified with the points on cusp manifolds for
which the second eigenvalue obeys the conditions:

  * 1.

the second eigenvalue field is negative,

$\displaystyle\ \ \ \ \ \ \ \ \lambda_{2}(\bm{q}_{c})<0\,,$ (44)

  * 2.

has a vanishing directional first derivative

$\displaystyle\ \ \ \ \ \ \ \
\bm{v}_{2}(\bm{q}_{c})\cdot\nabla\lambda_{2}(\bm{q}_{c})$ $\displaystyle=0\,,$
(45)

  * 3.

and has a positive second directional derivative

$\displaystyle\ \ \ \ \ \ \ \
\bm{v}_{2}(\bm{q}_{c})^{T}[\mathcal{H}\lambda_{2}(\bm{q}_{c})]\bm{v}_{2}(\bm{q}_{c})$
$\displaystyle>0\,.$ (46)

Note that the conditions $\lambda_{1}(\bm{q}_{c})=1/b_{+}(t_{c})$ and
$\lambda_{2}(\bm{q}_{c})<0$ indicate that – in the Zel’dovich approximation –
at time $t_{c}$ the mass elements only undergone a single shell-crossing
event. Non-linear gravitational collapse may include a set of further shell
crossings.

###  4.4 Filament Centers – Umbilic and Swallowtail Nodes

The stretching condition of a negative second eigenvalue, $\lambda_{2}<0$,
ensures that the central region of the filament is far away from umbilic
clusters. This may be directly inferred from the observation that umbilics are
specified by the simultaneous requirement – at some time $t_{0}$ – on the
first and second eigenvalue,

$\displaystyle 1-b_{+}(t)\lambda_{1}=0\,,1-b_{+}(t)\lambda_{2}=0\,,$

for some $t<t_{0}$. Evidently, these points are substantially removed from the
minimum stretching criterion listed above.

A similar observation can be made for swallowtail clusters. That these are
well separated from the filament’s stretching center can be inferred from the
strong alignment between tangent $\bm{T}$ and second eigenvector $\bm{v}_{2}$
guarantees the stretching center to be well separated from the swallowtail
clusters. In the swallowtail caustic the Lagrangian map satisfies the
condition

$\displaystyle\bm{v}_{1}(\bm{q}_{c})\cdot\nabla(\bm{v}_{1}(\bm{q}_{c})\cdot\nabla\lambda_{1}(\bm{q}_{c}))=\bm{v}_{1}(\bm{q}_{c})\cdot\bm{n}=0\,.$
(47)

This condition implies the alignment of the first eigenvector field
$\bm{v}_{1}$ with the tangent vector $\bm{T}$ field. In combination with an
approximation of $\bm{T}$ being parallel to $\bm{v}_{2}$ leads to the
conclusion that the stretching center of a filament is maximally distanced
from the swallowtail caustics $A_{4}$.

Figure 5: The caustic skeleton (blue) and the mirror skeleton (red) with the
corresponding cusp lines in two $N$-body simulations in Eulerian space.

###  4.5 Stretching Filaments and Mirror Crossings

In addition to the formal definition of the trunk and stretching center of a
filament, there is also an alternative, interesting, insightful, and perhaps
even surprising interpretation of the proposed condition for the stretching
center of a filament.

To arrive at this interpretation, we notice that the conditions

$\displaystyle\bm{v}_{1}\cdot\nabla\lambda_{1}$ $\displaystyle=0\,,$ (48)
$\displaystyle\bm{v}_{2}\cdot\nabla\lambda_{2}$ $\displaystyle=0\,,$ (49)

display a neat symmetry between the first and the second eigenvalue fields
$\lambda_{1}$ and $\lambda_{2}$. The key observation is that when we take the
negative of the displacement field,

$\Psi\mapsto-\Psi\,,$ (50)

also the eigenvalue fields flip sign,

$\displaystyle\lambda_{1}\mapsto-\lambda_{1}\,,$
$\displaystyle\lambda_{2}\mapsto-\lambda_{2}\,.$ (51)

Given that this entails the initial Gaussian conditions that flipped sign, on
the basis of these flipped eigenvalue fields we may infer the mirror caustic
skeleton. It defines the mirror cosmic web. A telling illustration of a
corresponding cosmic web realization and its mirror cosmic web is shown in
figure 5.

Inspecting the corresponding eigenvalue conditions (discussed above), we may
directly see that the central stretching regions of filaments correspond to
the intersection of filaments in the caustic skeleton and mirror caustic
skeleton. In a sense, the cosmic web and its mirror web may be seen as dual
patterns, in which the nodes of the cosmic web correspond to voids in the
mirror web, and the edges and web and the mirror web cross roughly at the
midpoints. In this interpretation, clusters and filaments have the same role
as that of vertices and edges in a dual (geometric) tessellation. This
realization resembles the duality between the Voronoi tessellation and its
dual Delaunay tessellation (Okabe et al., 2000; Icke & van de Weygaert, 1987;
van de Weygaert, 1994). Recently, Pontzen et al. (2016); Shim et al. (2021);
Stopyra et al. (2021); Shim et al. (2023); Stopyra et al. (2024) used a
similar realization to identify voids in the cosmic web with clusters in the
mirror web.

##  5 Filament Statistics

The caustic conditions enable a statistical analysis of the walls, filaments
and clusters in the present-day cosmic web using a combination of stochastic
geometry and Gaussian random field theory Rice (1944, 1945); Longuet-Higgins
(1957); Adler (1981); Bardeen et al. (1986); Adler & Taylor (2009). This may
ultimately shed light on several fundamental questions of structure formation
and enable the exploitation of the substantial amount of cosmological
information stored in the intricate and complex structure of the cosmic web.
After all, it would allow us to obtain insight into the dependence of filament
and cosmic web properties on the global cosmology.

Some of the most prominent key issues that are in need of a firm quantitative
statistical assessment, in terms of distribution functions and statistical
moments, are:

  * 1)

at what time do filaments form?

  * 2)

what is the number density of filaments?

  * 3)

what is the mean curvature, areas, and lengths of filaments?

  * 4)

what is the topology of the cosmic web over cosmic time?

  * 5)

how are these filament properties influenced by the power spectrum of the
primordial density perturbations?

In this paper, we study the statistical properties of filaments, extending the
work presented in Feldbrugge et al. (2016, 2019, 2023). We limit the analysis
to the simplest statistics, exploring how the caustic conditions yield solid
analytical predictions.

A full statistical analysis that extends these results to three-dimensional
fluids is the subject of an upcoming paper. In the 3D study, the extended
results include a more detailed analysis including the mean area and curvature
of walls, the mean length and curvature of filaments, the merging history of
filaments, and the implied number density of cluster nodes.

###  5.1 Potential derivatives: Definitions

The caustic skeleton describes the geometry of the cosmic web in terms of the
eigenvalues and eigenvectors of the Hessian $T_{ij}$ of the displacement
potential $\Psi(\bm{q})$,

$\displaystyle T_{ij}(\bm{q}_{c})\,=\,\frac{\partial}{\partial
q_{i}}\frac{\partial}{\partial q_{j}}\Psi_{\sigma}(\bm{q}_{c})\,,$ (52)

The eigenvalue $\lambda_{i}$ fields and eigenvector $\bm{v}_{i}$ fields
represent fields of non-linear functions of the second-order derivatives
$T_{ij}$ of the displacement potential. Also, by Leibniz’s rule, the gradients
of the eigenvalue fields, $\nabla\lambda_{i}$, are non-linear functions of the
second and third-order derivatives at a point. In addition, the normal of the
cusp filament

$\bm{n}=\nabla(\bm{v}_{i}\cdot\nabla\lambda_{i})\,,$ (53)

includes fourth-order derivatives of the displacement potential. In the
following, we denote the higher order partial derivatives of the smoothed
displacement potential $\Psi_{\sigma}$ at the point $\bm{q}_{c}$ as

$\displaystyle T_{ij\dots k}=\frac{\partial}{\partial
q_{i}}\frac{\partial}{\partial q_{j}}\dots\frac{\partial}{\partial
q_{k}}\Psi_{\sigma}(\bm{q}_{c})\,.$ (54)

###  5.2 Potential Derivatives: Gaussian Statistics

With the displacement potential representing a Gaussian random field, also
linear functionals of the potential are Gaussian random fields. This includes
derivatives of the potential $T_{ij\dots k}$ as well as convolutions of the
potential field. Moreover, any linear combination of Gaussian variates has a
Gaussian functional distribution. The probability distribution function
$p(\bm{Y})$ for such a set of Gaussian variates $\bm{Y}\in\mathbb{R}^{n}$
follows the multi-dimensional Gaussian distribution

$\displaystyle
p(\bm{Y})\,=\,\frac{\exp\left[-\frac{1}{2}\Delta\bm{Y}^{T}M^{-1}\Delta\bm{Y}\right]}{[(2\pi)^{n}\det
M]^{1/2}}\,.$ (55)

in which $\Delta\bm{Y}$ specifies the deviation from the mean,

$\Delta\bm{Y}=\bm{Y}-\langle\bm{Y}\rangle\,,$ (56)

and $M$ is the covariance matrix

$\displaystyle
M=\text{cov}(\bm{Y},\bm{Y})=\langle\Delta\bm{Y}^{T}\Delta\bm{Y}\rangle\,.$
(57)

For example, the second order derivatives $T_{ij}$,

$\bm{Y}_{2}=(T_{11},T_{12},T_{22})\,,$ (58)

are normally distributed random variables with zero mean, and covariance
matrix $M_{2,2}$,

$\displaystyle
M_{2,2}=\left\langle\Delta\bm{Y}_{2}^{T}\Delta\bm{Y}_{2}\right\rangle=\frac{\sigma_{2}^{2}}{8}\begin{pmatrix}3&0&1\\\
0&1&0\\\ 1&0&3\end{pmatrix}.$ (59)

where the $\sigma_{2}$’s is the second order spectral moment of the primordial
displacement potential $\Psi_{\sigma}(\bm{q})$ (see eq. (29)). Likewise, the
third order derivatives $\bm{T}_{ijk}$,

$\bm{Y}_{3}=(T_{111},T_{112},T_{122},T_{222})\,,$ (60)

comprise a set of Gaussian random variables with covariance matrix $M_{3,3}$,

$\displaystyle
M_{3,3}=\left\langle\Delta\bm{Y}_{3}^{T}\Delta\bm{Y}_{3}\right\rangle=\frac{\sigma_{3}^{2}}{16}\begin{pmatrix}5&0&1&0\\\
0&1&0&1\\\ 1&0&1&0\\\ 0&1&0&5\end{pmatrix}.$ (61)

Following this further for the fourth order derivatives $T_{ijkl}$

$\bm{Y}_{4}=(T_{1111},T_{1112},T_{1122},T_{1222},T_{2222})\,,$ (62)

the corresponding Gaussian distribution is specified by the covariance matrix
$M_{4,4}$,

$\displaystyle
M_{4,4}=\left\langle\Delta\bm{Y}_{4}^{T}\Delta\bm{Y}_{4}\right\rangle=\frac{\sigma_{4}^{2}}{128}\begin{pmatrix}35&0&5&0&3\\\
0&5&0&3&0\\\ 5&0&3&0&5\\\ 0&3&0&5&0\\\ 3&0&5&0&35\\\ \end{pmatrix}.$ (63)

In these Gaussian distributions, the influence of the power spectrum of the
primordial fluctuations enters via the spectral moments
$\sigma_{2},\sigma_{3},$ and $\sigma_{4}$ defined by (29).

In all, we may consider the Gaussian distribution for a full set of second-,
third-, and fourth-order derivatives of the displacement potential

$Y\,=\,(T_{11},\dots,T_{2222})\,,$ (64)

which corresponds to a Gaussian distribution with a vanishing mean and
covariance matrix

$\displaystyle M=\begin{pmatrix}M_{2,2}&0&M_{2,4}\\\ 0&M_{3,3}&0\\\
M_{2,4}^{T}&0&M_{4,4}\end{pmatrix},$ (65)

in which $M_{i,j}$ are the covariance matrices of the $i$-th and $j$-th order
derivatives of the displacement potential, which includes $M_{2,2}$, $M_{3,3}$
an $M_{4,4}$ specified above in equations (59), (61), (63) and with the
additional matrix $M_{2,4}$,

$\displaystyle M_{2,4}$
$\displaystyle=-\frac{\sigma_{3}^{2}}{16}\left(\begin{array}[]{ccccc}5&0&1&0&1\\\
0&1&0&1&0\\\ 1&0&1&0&5\\\ \end{array}\right),$ (69)

With respect to the Gaussian distribution for the combined second, third, and
fourth-order gradients we may note from the covariance matrix (eqn. (65)) that
the third-order gradients and the second-order and fourth-order gradients are
independently distributed. In general, two derivatives of a Gaussian random
field $T_{i_{1}j_{1}}$, $T_{i_{2}j_{2}}$ at a given point $\bm{q}_{c}$ are
independently distributed in case $i_{1}+i_{2}$ or $j_{1}+j_{2}$ is odd are
independently distributed. It manifests itself in the checkerboard pattern in
the covariance matrix, such as for example seen in expression (65).

###  5.3 Filament number density

One of the most prominent aspects of the filament population concerns the
issue of the number of filaments and the evolution of the number of filaments.
Here, we equate the number density of filaments to that of the number of
centers of filament trunks (see section 4).

Figure 6: The number density of filament centers as a function of their
formation times $b_{c}$ and the smoothing scale $\sigma=0,0.25,0.5,\dots,2$
from top to bottom for a scale-invariant power spectrum with a Gaussian cutoff
at $R_{s}=1$.

Following Rice’s formula (Rice, 1944, 1945; Longuet-Higgins, 1957; Adler,
1981; Bardeen et al., 1986; Adler & Taylor, 2009), the number density of
filament trunk centres, at length scale $\sigma$, forming at growing mode
$b_{c}$ is given by the expectation value

$\displaystyle\mathcal{N}_{\sigma}(b_{c})=b_{c}^{-2}\bigg{\langle}\left|\det\begin{pmatrix}\partial_{1}(\bm{v}_{1}\cdot\nabla\lambda_{1})&\partial_{2}(\bm{v}_{1}\cdot\nabla\lambda_{1})\\\
\partial_{1}(\bm{v}_{2}\cdot\nabla\lambda_{2})&\partial_{2}(\bm{v}_{2}\cdot\nabla\lambda_{2})\\\
\end{pmatrix}\right|$
$\displaystyle\times\delta_{D}^{(1)}(\bm{v}_{1}\cdot\nabla\lambda_{1})\delta_{D}^{(1)}(\bm{v}_{2}\cdot\nabla\lambda_{2})\delta_{D}^{(1)}(\lambda_{1}-b_{c}^{-1})$
$\displaystyle\times\Theta_{H}(-\lambda_{2})\Theta_{H}\left(\bm{v}_{2}^{T}[\mathcal{H}\lambda_{2}]\bm{v}_{2}\right)\bigg{\rangle}$
(70)
$\displaystyle=b_{c}^{-2}\bigg{\langle}\left|\det\begin{pmatrix}T_{1111}+\frac{3T_{112}^{2}}{\lambda_{1}-T_{22}}&T_{1112}+\frac{3T_{112}T_{122}}{\lambda_{1}-T_{22}}\\\
T_{1222}-\frac{3T_{112}T_{122}}{\lambda_{1}-T_{22}}&T_{2222}-\frac{3T_{122}^{2}}{\lambda_{1}-T_{22}}\\\
\end{pmatrix}\right|$
$\displaystyle\times\pi|\lambda_{1}-T_{22}|\delta_{D}^{(1)}(T_{12})\Theta_{H}(-T_{22})\delta_{D}^{(1)}(T_{111})\delta_{D}^{(1)}(T_{222})$
$\displaystyle\times\Theta_{H}\left(T_{2222}-\frac{2T_{122}^{2}}{\lambda_{1}-T_{22}}\right)\delta_{D}^{(1)}(T_{11}-b_{c}^{-1})\bigg{\rangle}\,,$
(71)

In the second equality, we change coordinates from the eigenvalue fields to
the derivatives of the displacement field in the eigenframe. The resulting
number density curves are plotted in fig. 6 for smoothing scales
$\sigma=0,0.25,\dots,2$. The number density is the largest at the smallest
smoothing scale around the formation time $b_{c}=0.3$. As the smoothing scale
increases, the number density of filament trunks gradually decreases. Also, we
see that they form at later times, corresponding to larger growing modes.
Moreover, in the next section, we show that the large filaments emerge over a
wider range of formation times.

###  5.4 Formation time of a filament

Considering the emergence of the centre of the filament’s trunk as the
formation time of a filament, we may assess the formation times of the entire
filament population and explore its dependence on scale, smoothing scale,
power spectrum, and other characteristics of the filaments.

In our evaluation of the filament formation time properties, we therefore
focus on the emergence of the centers of the filaments. The caustic skeleton
framework provides a detailed quantitative definition for this. It is
identified with the moment shell-crossing occurs. Specifying time $t$ in terms
of the growing mode of structure formation, $b_{+}(t)$, for a location
$\bm{q}_{c}$ caustic conditions specify the occurrence of shell-crossing to
happen at growing mode factor $b_{c}=b_{+}(t_{c})$,

$\lambda_{1}(\bm{q}_{c})=1/b_{c}\,,$ (72)

with the additional condition that

$\bm{v}_{1}(\bm{q}_{c})\cdot\nabla\lambda_{1}(\bm{q}_{c})=0\,.$ (73)

Restricting our attention to specific filament center locations, we have the
additional requirement of a mass element at that location to expand maximally
along the ridge of the filaments. This means that

  * 1.

the second eigenvalue is negative

$\hskip 21.33955pt\lambda_{2}(\bm{q}_{c})<0\,,$ (74)

  * 2.

its first derivative in the $\bm{v}_{2}$ direction vanishes

$\hskip 21.33955pt\bm{v}_{2}(\bm{q}_{c})\cdot\nabla\lambda_{2}(\bm{q}_{c})=0$
(75)

  * 3.

its second-order derivative is positive,

$\hskip
21.33955pt\bm{v}_{2}(\bm{q}_{c})^{T}[\mathcal{H}\lambda_{2}(\bm{q}_{c})]\bm{v}_{2}(\bm{q}_{c})>0\,.$
(76)

Given the isotropy of the Gaussian random field, for the mathematical
treatment of the problem it is convenient to rotate to the eigenframe of the
Hessian $\mathcal{H}\Psi_{\sigma}(\bm{q}_{c})$. It allows us to infer the
corresponding identities without loss of generality. For the details of the
derivation of the corresponding identities, we refer to Feldbrugge & van de
Weygaert (2023). In terms of the second-order derivatives $T_{ij}$ of the
displacement potential the eigenvalues within this coordinate system are given
by the identities

$\displaystyle\lambda_{1}(\bm{q}_{c})$ $\displaystyle=T_{11}\,,$
$\displaystyle T_{12}$ $\displaystyle=0,$ (77)
$\displaystyle\lambda_{2}(\bm{q}_{c})$ $\displaystyle=T_{22}\,.$

The directional derivatives of the eigenvalue fields coincide with the third-
order derivatives of the displacement field

$\displaystyle\bm{v}_{1}\cdot\nabla\lambda_{1}$ $\displaystyle=T_{111}\,,$
(78) $\displaystyle\bm{v}_{2}\cdot\nabla\lambda_{2}$
$\displaystyle=T_{222}\,,$

while the second-order derivative of the second eigenvalue field $\lambda_{2}$
in the $\bm{v}_{2}$ direction yields the identity

$\displaystyle\bm{v}_{2}^{T}[\mathcal{H}\lambda_{2}]\bm{v}_{2}=T_{2222}-\frac{2T_{122}^{2}}{\nu-
T_{22}}\,.$ (79)

Using these inferred identities, we find that the probability distribution
function of the time at which the center of the filament emerges, in terms of
the corresponding growing mode $b_{c}$, is equivalent to the conditional
probability $p(\lambda_{1}\,|\,\dots)$ for the first eigenvalue $\lambda_{1}$,

$\displaystyle p(\text{center of filament emergest at growing mode }b_{c})$
(80)
$\displaystyle=\lambda_{1}^{2}p(\lambda_{1}\,|\,\lambda_{2}<0<\lambda_{1},\bm{v}_{1}\cdot\nabla\lambda_{1}=\bm{v}_{2}\cdot\nabla\lambda_{2}=0)\big{|}_{\lambda_{1}=1/b_{c}}\,,$

involving the conditions specified in the listing above. For simplicity we
dropped the positivity condition
$\bm{v}_{2}^{T}[\mathcal{H}\lambda_{2}]\bm{v}_{2}>0$, as a critical point of
$\lambda_{2}$ in the $\bm{v}_{2}$ direction is likely to be minimum when
$\lambda_{2}<0$. In this expression, the pre-factor $\lambda_{1}^{2}$ results
from the Jacobian of the transformation $b_{c}=1/\lambda_{1}$, preserving
probability i.e.,
$|p(b_{c})\mathrm{d}b_{c}|=|p(\lambda_{1})\mathrm{d}\lambda_{1}|$.

Figure 7: The PDF of the formation time $b_{c}$ of the centers of the
filaments as a function of the smoothing scale $\sigma$ for a scale-invariant
power spectrum with a Gaussian cutoff at $R_{s}=1$. The peak of the
distribution is plotted by the black curve.

Since the eigenvalue field and its first order derivatives along the
eigenvector directions are independent – following from their relations to the
derivatives of the primordial displacement potential in the eigenframe – we
find that

$\displaystyle p(\text{center of filament forms a growing mode }b_{c})$ (81)
$\displaystyle=\lambda_{1}^{2}p(\lambda_{1}|\lambda_{2}<0<\lambda_{1})\big{|}_{\lambda_{1}=1/b_{c}}$
(82)
$\displaystyle=\frac{\lambda_{1}^{2}\int_{-\infty}^{0}p(\lambda_{1},\lambda_{2})\mathrm{d}\lambda_{2}}{\int_{-\infty}^{0}\int_{0}^{\infty}p(\lambda_{1},\lambda_{2})\mathrm{d}\lambda_{1}\mathrm{d}\lambda_{2}}\bigg{|}_{\lambda_{1}=1/b_{c}}$
(83)
$\displaystyle=\frac{2\lambda_{1}^{2}\Theta_{H}(\lambda_{1})}{3\sigma_{2}^{2}}\bigg{[}2\lambda_{1}e^{-\frac{4\lambda_{1}^{2}}{3\sigma_{2}^{2}}}\text{erfc}\left(\frac{\lambda_{1}}{\sqrt{6}\sigma_{2}}\right)+\sqrt{\frac{6}{\pi}}\sigma_{2}e^{-\frac{3\lambda_{1}^{2}}{2\sigma_{2}^{2}}}\bigg{]}\bigg{|}_{\lambda_{1}=1/b_{c}}\,,$
(84)

using the Doroshkevich formula for the distribution of the eigenvalue fields
(Doroshkevich, 1970; Feldbrugge et al., 2023)

$\displaystyle p(\lambda_{1},\lambda_{2})=$
$\displaystyle\sqrt{\frac{2}{\pi}}\frac{2}{\sigma_{2}^{3}}(\lambda_{1}-\lambda_{2})e^{-\frac{1}{2\sigma_{2}^{2}}\left(3(\lambda_{1}+\lambda_{2})^{2}-8\lambda_{1}\lambda_{2}\right)}$
$\displaystyle\times\Theta_{H}(\lambda_{1}-\lambda_{2})\,.$ (85)

Figure 7 visually summarizes the resulting filament formation time. It shows
the probability value as a function of formation time $b_{c}$ and smoothing
scale $\sigma$. From the inferred formation time distribution, we may conclude
that in the case of no smoothing, i.e., $\sigma=0$, the filaments are likely
to form at a growing mode of $b_{c}\approx 0.36$. Most of the filament centers
have emerged at the current time $b_{c}=1$.

Another interesting observation is that filaments are more likely to form at
later times as the smoothing length $\sigma$ increases. Also, at higher values
of $\sigma$, we see a larger variation in formation time: large filaments have
a larger diversity in formation time and history. Moreover, at smoothing
length $\sigma=2$ a filament is unlikely to have formed yet in the cosmic web
at the current epoch.

##  6 Filaments in the cosmic web:
Constrained Realizations

Based on the dynamical definition of a cosmic filament, in principle, we may
use the corresponding conditions to generate constrained simulations of
filaments within any cosmological setting. In practice, we may do so by
investigating the primordial conditions and configurations out of which
filaments emerged given the primordial density and gravity field is a Gaussian
random field. Such a setup would enable a targeted "laboratory" study that
explores the setting, connections, characteristics, development, and evolution
of filaments in the cosmic web. The constrained filament simulations will
allow a systematic study of these aspects as a function of the spatial scale
$\sigma$ of the filaments and cosmic epoch (specified by growing mode
$b_{c}$). Particularly important is that such a setup would allow a systematic
exploration of the influence of background cosmology and the power spectrum of
primordial fluctuations on emerging filaments. It will be the most efficient
path towards gaining insight into the cosmological information content of the
cosmic web, and recognizing the cosmological signatures in its observed
structure and properties.

The formalism for constrained filament simulations has been described in
detail in Feldbrugge & van de Weygaert (2023). It is an elaboration of the
linear constraint Hoffman-Ribak formalism (Hoffman & Ribak, 1991) for density
and gravity fields towards non-linear constraints, following the definitions
and relations outlined in van de Weygaert & Bertschinger (1996). The non-
linear constraints are non-linear functions of Gaussian random field
variables. Tidal or deformation eigenvalues, and their gradients, are telling
examples of such non-linear functions of Gaussian quantities, in this case,
the components of the Hessian of the gravitational potential.

In the following, we outline the formalism for the constrained generation of
cosmic filaments in the cosmic web, based on their physical definition. The
latter includes the phase-space-based caustic conditions specifying their
spine in terms of the cusp $A_{3}$ singularity (Feldbrugge et al., 2018), in
combination with the tidal dilation of the mass content along their ridge
(section 4).

Figure 8: Six realizations of the cusp filaments with the formation time
$b_{c}=0.5$ at the scale $\sigma=0.5$. We plot the $N$-body particles and the
initial mesh on the corresponding density field $\log(\rho+1)$.

Figure 9: The median density field of the realizations of the centers of the
filament as a function of the formation time $b_{c}=0.3,0.5,0.8,1.0$ from left
to right and length scale $\sigma=0.1,0.5,1.0$ from top to bottom. We plot the
median density field on a logarithmic scale $\log(\rho_{median}+1)$.

Figure 10: The density profile perpendicular to the filament as a function of
formation time $b_{c}=0.3,0.5,0.8,1$ (blue yellow green and red) for the
length scale $\sigma=0.1,0.5,1.0$ (left to right).

Figure 11: The distribution of the number of streams as a function of the
formation time for the length scale $\sigma=0.1,0.5,1.0$ (left to right).

###  6.1 Constrained filament initial conditions:
specification

Initial conditions for filaments consist of two aspects: the cusp $A_{3}$
caustic constraints, in combination with supplementary constraints on the
tidally induced stretching of its mass content. Following the previous
section, we focus the filament constraint on the point representing the
stretching center of the filament.

In the eigenframe of the Hessian $\mathcal{H}\Psi_{\sigma}$ the cusp
conditions $\lambda_{1}=1/b_{c},$ and $\bm{v}_{1}\cdot\nabla\lambda_{1}=0$
(see sect. 5.4) assume the form

$\displaystyle T_{11}$ $\displaystyle=1/b_{c}\,,\,$ $\displaystyle T_{12}$
$\displaystyle=0\,,\,$ $\displaystyle T_{22}$ $\displaystyle<T_{11}\,,\,$ (86)
$\displaystyle T_{111}$ $\displaystyle=0\,.$

The normal $\bm{n}$ of the cusp line in Lagrangian space is given by

$\bm{n}=\left(T_{1111}+\frac{3T_{112}^{2}}{T_{11}-T_{22}},T_{1112}+\frac{3T_{112}T_{122}}{T_{11}-T_{22}}\right)\,.$
(87)

Along the cusp manifold, a mass element undergoes maximal stretch when it
fulfills the additional conditions $\lambda_{2}<0$,
$\bm{v}_{2}\cdot\nabla\lambda_{2}=0$,
$\bm{v}_{2}^{T}[\mathcal{H}\lambda_{2}]\bm{v}_{2}>0$ which in the eigenframe
assumes the form

$\displaystyle T_{22}$ $\displaystyle<0\,,\,$ $\displaystyle T_{222}$
$\displaystyle=0\,,\,$ (88) $\displaystyle 0$
$\displaystyle<T_{2222}-\frac{2T_{122}^{2}}{T_{11}-T_{22}}\,.$

###  6.2 Constrained filament initial conditions:
procedure

We generate initial conditions for a set of $N$-body simulations by sampling
Gaussian random fields subject to the constraints specified above. Outlined in
detail in Feldbrugge & van de Weygaert (2023) this is accomplished by the
combination of a rejection sampling scheme for the non-linear function of
second- third- and fourth-order derivatives of the smoothed displacement
potential with that of the Hoffman-Ribak algorithm for sampling Gaussian
random field Bertschinger (1987); Hoffman & Ribak (1991); van de Weygaert &
Bertschinger (1996). It results in the following procedure:

  1. 1.

Sample the statistic $Y=(T_{11},\dots,T_{2222})$. The probability of $Y$ is
that of a multidimensional normal distribution with vanishing mean and
covariance matrix (65) subject to the linear constraints

$(T_{11},T_{12},T_{111},T_{222})=(1/b_{+}(t_{c}),0,0,0)\,,$ (89)

multiplied with the Jacobian factor $|T_{11}-T_{22}|$.

  2. 2.

Reject samples that do not satisfy the conditions

$\displaystyle T_{22}<0\,,\quad
T_{2222}-\frac{2T_{122}^{2}}{T_{11}-T_{22}}>0\,.$ (90)

  3. 3.

For each remaining sample of $Y$ in the eigenframe, we evaluate the normal of
the cusp manifold in the eigenframe using expression (87),

  4. 4.

Subsequently, the derivatives of the deformation potential are rotated to
align the normal of the cusp manifold $\bm{n}$ with the vertical axis,
following the transformations derived in appendix E of Feldbrugge & van de
Weygaert (2023).

  5. 5.

For each sampled $Y$, we sample a corresponding Gaussian random field with
power spectrum $P_{\Psi}(k)$ using the Hoffman-Ribak algorithm (Bertschinger,
1987; Hoffman & Ribak, 1991; van de Weygaert & Bertschinger, 1996).

###  6.3 Constrained Filament Inventory

Having specified filament constraints for the primordial Gaussian density and
potential fields, we invoke these for an extensive and systematic study of the
properties and large-scale environment of the specified filaments in their
cosmological setting.

Here we make an inventory of filaments for a range of configurations within a
global cosmological setting of an Einstein-de Sitter Universe. Our interest
concerns in particular the the morphology of the filament configurations and
the multi-stream nature of the emerging filaments.

####  6.3.1 Filament Configurations

For each constrained filament configuration we set up $1000$ random
realizations in a two-dimensional box of size $25R_{s}\times 25R_{s}$, with
$R_{s}$ the spectral Gaussian cutoff scale (see eqn. (20)), specifying the
presence of the specific filaments at the centre of the box according to the
constrained field formalism outlined above. We also arrange the filament
orientation, such that all realizations contain a filament running through the
center of the box, oriented in the vertical direction. This is arranged using
a rather special procedure. For each sampled initial condition, two different
$N$-body simulations are run:

  1. 1.

The first simulation evolves the sampled initial conditions. The second
simulation evolves a smoothed version of the initial conditions, convolving
the entire field with a Gaussian kernel (27) whose filter scale $\sigma$ is
equal to that of the imposed filament constraint.

  2. 2.

From the smoothed $N$-body simulation, we evaluate the position
$\bm{x}_{t}(\bm{q}_{c})$ and orientation of the cusp filament in Eulerian
space using the equation $\nabla\bm{x}_{t}(\bm{q}_{c})\bm{T}$. Subsequently,
we center and orient the filament in the unsmoothed $N$-body simulation along
the vertical direction using the position and orientation of the filament in
the smoothed $N$-body simulation.

The filament configurations of which we generate each $1000$ realizations have
the following range of smoothing lengths $\sigma$: $\sigma=0.1,0.5,1.0$ In
addition, for each of these we explore five different amplitudes, specified in
terms of the implied formation (emergence) time $b_{c}=0.3,0.5,0.8,1.0$.

####  6.3.2 Filament Realizations

A representative sample of six filament realizations is shown in figure 8. The
density and particle map concern the matter distribution in and around the
specified filament at the current epoch.

The filament realization in the four central and righthand panels show fully
developed filaments, assembled around the stretching center at the center of
the box, with their trunk neatly oriented along the vertical axis. While the
mass distribution at the top lefthand panel stretches along a long vertical
ridge, it has not yet reached the advanced stage in which it developed into a
multi-stream region. In other words, the filament has not yet emerged as such.
In this situation, the Zel’dovich approximation apparently did not properly
predict the filament formation time.

A more complex situation is seen in the bottom lefthand panel. Strong non-
linear evolution of the surrounding mass distribution has entangled the
specified central filament, oriented vertically, with a horizontally running
filament. A massive and compact cluster has formed at their crossing point, at
the center of the simulation box: the specified filament is one of the
principal high-density filamentary extensions of the cluster. We should also
note that the cusp filament of the Zel’dovich approximation is unable to fully
capture the highly non-linear dynamical evolution of this region. In
principle, one might opt for a post-processing step in which such
configurations are removed from further consideration, by e.g. restricting the
filamentary regions in the smoothed $N$-body simulation to $3$\- or $5$-stream
regions. In the present analysis, we do not perform such a post-selection.

####  6.3.3 Filament Density Profiles

On the basis of the large, statistically representative, sample of constrained
filament realizations, we assess the generic mass distribution in the vicinity
of filaments. The density field around each filament realization is evaluated
on a lattice around the stretching center of the filament. They are computed
by means of the high-resolution self-adaptive phase-space Delaunaty
Tessleation Field Estimator formalism (Feldbrugge (2024) see appendix A for an
outline).

For each set of $1000$ filament simulation realizations we determine the
median density, as a function of the specified smoothing scales and formation
times. The median density profile is preferred over the mean density profile,
as the latter tends to be dominated by the caustics where the density spikes
to infinity. The median density is insensitive to these outliers. In addition,
the median also removes effects from the filamentary cluster extensions.
Figure 11 presents the median filament mass distribution in a matrix of 12
panels, with formation time running along the horizontal direction – for
$b_{c}=0.3,0.5,0.8,1.0$ – and with smoothing scale running along the vertical
direction – for $\sigma=0.1,0.5,1.0$.

The median filamentary mass distribution has an oval shape extending along the
vertical direction. Going from the lefthand panels to the right, i.e., as we
proceed from filaments forming early to ones forming later, we see that
filaments that form later are thinner and have a concentrated density spike at
the center. We may understand this as filaments that form at early times have
more time to develop and grow into a denser structure phase-mixed structure.
Proceeding from small-scale to large-scale filaments, i.e., from the top to
the bottom panels, we find that larger smoothing lengths $\sigma$ entail
thicker, more massive filaments. Moreover, on average they define more
extended structures along the vertical direction.

The density profiles along the width of the filaments are presented in figure
11. These underline the impression given by the median filament mass
distribution. The filaments that form at earlier times have more time to phase
mix, and produce thicker objects that have a density plateau, followed by a
rapid decline in density. By contrast, filaments that form at a later time end
up with a more peaked density profile.

####  6.3.4 Filament Multi-stream Character

The PS-DTFE formalism also allows us to determine the local streaming
structure, i.e., the local number of streams. We find that in the majority of
realizations, the central stretching region of the filament is situated in a
$3$\- or $5$\- and even $7$-stream region. For three different values of
smoothing length – $\sigma=0.1,0.5,1.0$ – four bar histograms depict the
distribution of a number of streams for four different filament formation
times.

Visual inspection indicates that the configurations with $3$, $5$ or $7$
streams are filaments. For late formation times and small length scales, we
find that a significant fraction of realizations entails filaments whose
stretching center resides in a $1$-stream region. These configurations tend to
be similar to the upper left panel of fig. 8, where the filament is about to
shell-cross, but has not yet reached that phase. In these instances, the
Zel’dovich approximation underpredicts the first shell-crossing time. By
contrast, visual inspection reveals that regions with more than $9$ streams
tend to correspond to filaments whose stretching center resides in regions
that would classify as the extension of a cluster (cf. the lower left panel of
fig. 8).

##  7 Alternative filament conditions

In addition to the physical definition for the nature of cosmic filaments
forwarded in the present study, there is a range of different definitions in
the literature for what constitutes a filament. Here we compare the caustic
definition of filaments with that of the two most common definitions. Both
entail singularities in the matter distribution, be it singularities in the
Eulerian primordial matter distribution.

The alternative filament definitions that we consider here are of an explicit
topological nature, involving the connections between singularities in the
mass distribution. They are based on the observation that filaments are often
the bridges between clusters in the matter and galaxy distribution.
Mathematically speaking, the contention is that clusters are maxima in the
density field or, following a stronger assumption, that they emanate from
maxima in the primordial density field. Within this context, cosmic voids
descend from the minima in the primordial density field. Dynamically even more
relevant is the identification of the minima in the primordial gravitational
potential with the dominant mass concentrations in the Universe, and the
maxima of the primordial gravitational potential with the gravitationally
dominant voids in the large-scale matter distribution. A major issue within
this view remains the determination of the relation between the scale of
primordial singularities and that of the resulting structures in the matter
distribution, in particular in view of the hierarchical buildup of structure
(see e.g. van de Weygaert & Bond, 2008; Aragon-Calvo, 2024).

Following the identification of clusters with the maxima in the density field,
or with the minima in the gravitational potential field, filaments are seen as
the bridges that join the clusters. Mathematically they are then defined to be
the integral lines running between the clusters, passing through the
corresponding saddle point(s) of the density or potential field. The
eigenvectors of the Hessian of the primordial density field, or the Hessian of
the gravitational potential, at the saddle point, determine the orientation of
the proto-filament. Arguably the most widely known version of this view
concerns the density field. The cosmic skeleton model (Pogosyan et al., 2009b;
Codis et al., 2018a) analyzes this mathematical model of the cosmic web in
considerable detail, while it forms the incentive for the Disperse formalism
for cosmic web classification and identification (Sousbie, 2011; Sousbie et
al., 2011). A similar identification of filaments with the saddles in the
gravitational potential field had been forwarded in the earlier studies of the
linear constrained gravitational field realizations (van Haarlem & van de
Weygaert, 1993; van de Weygaert & Bertschinger, 1996).

These alternative filament definitions share a clear geometric and topological
similarity, linking the central regions of the filaments with saddle points.
However, they operate on different levels. Whereas the saddle point of the
primordial gravitational potential is governed by the first and second
derivatives of the primordial gravitational potential (or displacement
potential), the saddle points of the density field follow from the third- and
fourth-order derivatives of the density field. The caustic skeleton definition
proposed in this paper is most closely linked with the former definition, as
it is defined in terms of the same eigenvalues and eigenvectors of the Hessian
of the gravitational field.

Nonetheless, there are some major conceptual differences between the caustic
skeleton description and the specification in terms of saddle points. Most
important is that the saddle point definitions do not allow to set the
formation time. As such, they do not allow the differentiation of filaments
with different formation times and hence lack the ability to describe the
hierarchical buildup of the cosmic web. This also implies a problem of a more
practical nature, that of the inability to guarantee that the proposed proto-
filaments have formed by the current epoch.

Below we treat the constrained formalism to set up constrained filament
realization for the two alternatives of a filament specified on the basis of
the saddle point location and orientation in the primordial Gaussian
gravitational potential field, and that in the primordial Gaussian density
field.

Figure 12: Realizations of saddle points in the smoothed primordial density
perturbation at the scale $\sigma=0.5$. We plot the $N$-body particles and the
initial mesh on the corresponding density field $\log(\rho+1)$.

Figure 13: The median density field of the realizations of the saddle points
of the primordial density perturbation as a function of the length scale
$\sigma=0.1,0.5,1.0$ from left to right. We plot the median density field on a
logarithmic scale $\log(\rho_{median}+1)$. Figure 14: The distribution of the
number of streams as a function of the formation time for the length scale
$\sigma=0.1,0.5,1.0$ (left to right)

###  7.1 Saddle points of the primordial density

One of the most widely used methods for describing and/or analyzing the cosmic
web center around the assumption that filaments are ridges in the cosmic
density field, connecting clusters via the saddle points in the density field.
A variety of geometric and topological methods are based on this premise for
the identification of structure at the current epoch, such as main versions of
the MMF/Nexus algorithm (Aragón-Calvo et al., 2007a; Cautun et al., 2013) and
the Disperse topological formalism (Sousbie, 2011; Sousbie et al., 2011). The
related cosmic skeleton model emphasizes the central role of saddles in the
primordial density field as sites around which filaments form (Pogosyan et
al., 2009b; Codis et al., 2018a).

While this seems a reasonable assumption for any one moment in time, given the
hierarchical buildup of structure and corresponding multiscale nature of the
cosmic web it may not have the universal validity it might yield at first
sight. Cluster nodes and filamentary arteries may have evolved from proto-
structures at a variety of scales, so it may be far from straightforward to
identify saddles in the primordial density field with filaments in the current
cosmic mass distribution. In other words, while it may be useful for the
analysis of the structural pattern at a given cosmic epoch, it remains to be
assessed in how far it has any validity for considering the dynamical
evolution of the cosmic web.

In this section, we investigate, using constrained realizations, how far
saddle points in the primordial density field are related to the emerging
filaments in the cosmic web.

####  7.1.1 Saddle Point Definition

For the primordial density perturbation field smoothed by a Gaussian kernel at
length scale $\sigma$, the density perturbation can be written in terms of the
displacement potential using the Poisson equation

$\delta_{\sigma}(\bm{q}_{c})=\nabla^{2}\Psi_{\sigma}(\bm{q}_{c})=T_{11}+T_{22}\,.$
(91)

Imposing singularity conditions on the density field thus involves constraints
on the third-order and fourth-order gradients of the potential field: the
point $\bm{q}_{c}$ is a critical point of the primordial density perturbation
$\delta_{\sigma}$ when its gradient vanishes, i.e.,

$\displaystyle\nabla\delta_{\sigma}(\bm{q}_{c})$
$\displaystyle=(T_{111}+T_{122},T_{112}+T_{222})=\bm{0}\,.$ (92)

It is a saddle point when the determinant of the Hessian

$\displaystyle\mathcal{H}\delta_{\sigma}(\bm{q}_{c})$
$\displaystyle=\begin{pmatrix}T_{1111}+T_{1122}&T_{1112}+T_{1222}\\\
T_{1112}+T_{1222}&T_{1122}+T_{2222}\end{pmatrix}$ (93)

assumes a negative value (the Hessian has one positive and one negative
eigenvalue). For the constrained realizations that we generate for our
analysis, we orient the associated filament by requiring that the eigenvector
of the Hessian $\mathcal{H}\delta_{\sigma}(\bm{q}_{c})$ associated with the
positive eigenvalue is aligned with the vertical direction.

####  7.1.2 Implementation

For the generation of the density field-constrained conditions for a proto-
filament, the implementation steps concern the generation of properly sampled
third-order and fourth-order gradients of the deformation potential. To this
end, we use the fact that the third-order and fourth-order gradients and
independently distributed Gaussian variables. Hence, we sample them
independently in two subsequent steps.

  1. 1.

The third-order derivatives set up the space
$\bm{Y}_{3}=(T_{111},T_{112},T_{122},T_{222})$, and obey the linear
constraints

$\displaystyle T_{111}+T_{122}$ $\displaystyle=0\,,$ $\displaystyle
T_{112}+T_{222}$ $\displaystyle=0\,.$ (94)

To exploit the linear constraints on the third-order derivates, we consider
the third-order statistic
$Y=(T_{111}+T_{122},T_{112}+T_{222},T_{111}-T_{122},T_{112}-T_{222})$, which
is normally distributed $\mathcal{N}(\bm{0},M)$ with covariance matrix

$\displaystyle M=\frac{\sigma_{3}^{2}}{4}\begin{pmatrix}2&0&1&0\\\ 0&2&0&-1\\\
1&0&1&0\\\ 0&-1&0&1\end{pmatrix}.$ (95)

  2. 2.

Given the constraints $T_{111}+T_{122}=T_{112}+T_{222}=0$, the stochastic
variables $T_{111}-T_{122}$ and $T_{112}-T_{222}$ are independent and have an
identical normal distribution $\mathcal{N}(0,\sigma_{2}^{2}/8)$. Sampling two
variables $A$ and $B$ from this normal distribution, a representative sample
of the third-order derivatives $T_{ijk}$ is obtained via the relation,

$\displaystyle T_{111}$ $\displaystyle=A/2\,,$ $\displaystyle T_{122}$
$\displaystyle=-A/2\,,$ $\displaystyle T_{112}$ $\displaystyle=B/2\,,$ (96)
$\displaystyle T_{222}$ $\displaystyle=-B/2\,.$

  3. 3.

We sample the fourth-order derivatives space
$\bm{Y}_{4}=(T_{1111},\dots,T_{2222})$, which is normally distribution
$\mathcal{N}(\bm{0},M_{4,4})$, with zero mean and covariance matrix $M_{4,4}$,

$\displaystyle M_{4,4}=\frac{\sigma_{4}^{2}}{8}\begin{pmatrix}35&0&5&0&3\\\
0&5&0&3&0\\\ 5&0&3&0&5\\\ 0&3&0&5&0\\\ 3&0&5&0&35\\\ \end{pmatrix}.$ (97)

  4. 4.

Having sampled the fourth-order derivatives, we evaluate the Hessian
$\mathcal{H}\delta_{\sigma}(\bm{q}_{c})$ and reject samples for which its
determinant assumes a positive value.

  5. 5.

For the positive eigenvalue of the Hessian $\mathcal{H}\delta_{\sigma}$, the
eigenvector is rotated to the vertical direction, so that the third- and
fourth-order derivatives are aligned with the vertical direction. This is
accomplished following the transforms derived in appendix E of Feldbrugge &
van de Weygaert (2023).

  6. 6.

Subsequently, we generate a realization of the constrained field with the
Hoffman-Ribak algorithm for Gaussian random fields fulfilling a set of linear
constraints (Bertschinger, 1987; Hoffman & Ribak, 1991; van de Weygaert &
Bertschinger, 1996).

  7. 7.

Finally, for each sampled initial condition, two dark-matter-only 2D $N$-body
simulations are run (Hidding, 2020). The 2D simulations gave $256\times 256$
particles in a $25R_{s}\times 25R_{s}$ box in an Einstein-de Sitter cosmology.
The first simulation evolves the sampled initial conditions. The second
simulation evolves a smoothed version of the initial conditions with a
Gaussian kernel (27) and the smoothing length $\sigma$ at which the
constraints are set. The $N$-body simulation without the smoothing is centered
and aligned using the $N$-body simulation with the smoothed initial
conditions.

Figure 15: Realizations of saddle points in the smoothed primordial
gravitational potential at the scale $\sigma=0.5$. We plot the $N$-body
particles and the initial mesh on the corresponding density field
$\log(\rho+1)$.

Figure 16: The median density field of the realizations of the saddle points
of the primordial potential field as a function of the length scale
$\sigma=0.1,0.5,1.0$ from left to right. We plot the median density field on a
logarithmic scale $\log(\rho_{median}+1)$. Figure 17: The distribution of the
number of streams as a function of the formation time for the length scale
$\sigma=0.1,0.5,1.0$ (left to right)

####  7.1.3 Realizations

A random set of realizations of filaments implied by the imposed primordial
density field saddle constraints is presented in figure 13.

The figure panels clearly reveal the existence of a correlation between the
non-linear structure of the cosmic web and the saddle points of the primordial
density field. Nonetheless, the correlation is not as strong as in the case of
the caustic skeleton filament conditions. In the next subsection, we will see
it is also considerably weaker than the filament constraints emanating from
saddle conditions in the primordial potential field.

The density field saddle conditions yield configurations that may not even
resemble that of filaments. In some cases, the saddle point ends up in a
cluster (righthand panels fig. 13). In other cases, it may reside in a region
that might develop into a filament in the future, but would at the current
epoch most likely be classified as a void (lefthandand central panels). A
major factor is that the density field saddle constraints do not specify the
formation time so the time at which the filament emerges is left a free
parameter. It may also be noticed that the emerging structures do not reveal
any of the intended clear alignment with the vertical axis.

To assess the statistical properties of the proto-filaments implied by the
primordial density field constraints, we study the corresponding properties of
an ensemble consisting of $1000$ $N$-body simulations. This is accomplished
for the length scale $\sigma=0.1,0.5$ and $1.0$. When assessing the number of
streams in the Eulerian region in which the saddle resides, we find a
quantitative confirmation of the impressions obtained from figure 13. For a
large fraction of the simulations, the saddle point resides in a single stream
region. Visual inspection indicates that these are largely void regions. When
the saddle point resides in a $3$\- or $5$-stream region, we often observe a
filamentary structure. Nonetheless, only rarely this structure is aligned with
the intended vertical axis. This is emphasized by the median density profiles
of the realizations, shown in figure 13. Reflecting the fact that the
orientation of the proto-filaments is not really constrained, and in fact has
a random isotropic distribution, is a median density field that has an
approximate rotational symmetry.

###  7.2 Saddle points primordial gravitational potential

From a dynamical perspective, one may assume that the minima of the primordial
gravitational potential correspond to regions that are due to gravitational
collapse and form clusters, while the maxima will evolve into voidlike
regions. From this dynamical perspective, proto-filaments are the ridges in
the primordial gravitational potential landscape, centered around the saddle
point of the (displacement) potential field $\Psi$ (see e.g. van Haarlem & van
de Weygaert, 1993; van de Weygaert & Bertschinger, 1996). The eigenvectors of
its Hessian, $\mathcal{H}\Psi$, determine the orientation of the proto-
filament.

Mathematically, this translates into the following conditions with respect to
the saddle points in the smoothed primordial gravitational potential. For
these, the following conditions hold:

  * 1.

the gradient of the displacement potential vanishes

$\displaystyle\nabla\Psi_{\sigma}(\bm{q}_{c})=(T_{1},T_{2})=\bm{0}\,,$ (98)

  * 2.

the determinant of the Hessian

$\displaystyle\mathcal{H}\Psi_{\sigma}(\bm{q}_{c})=\begin{pmatrix}T_{11}&T_{12}\\\
T_{12}&T_{22}\end{pmatrix}$ (99)

assumes a negative value (having one positive and one negative eigenvalue).

We orient the proto-filament such that the eigenvector corresponding to the
positive eigenvalue of the Hessian $\mathcal{H}\Psi_{\sigma}(\bm{q}_{c})$ is
aligned vertically.

####  7.2.1 Implementation

Following the realization that the first-order derivatives $T_{i}$ and second-
order derivates $T_{jk}$ of a Gaussian random field themselves also define a
Gaussian random field, and the first-order and second-order derivatives are
independently distributed, we set up the following procedure for setting up a
proto-filament:

  1. 1.

The first step consists of sampling the second-order derivatives from the
normal distribution $\mathcal{N}(\bm{0},M_{2,2})$. The covariance matrix
$M_{2,2}$ is given by equation (59).

  2. 2.

Evaluate the Hessian $\mathcal{H}\Psi_{\sigma}(\bm{q}_{c})$, rejecting sampled
Hessians for which the determinant assumes a positive value.

  3. 3.

Evaluate the eigenvector of the Hessian corresponding to the positive
eigenvalue. Align the eigenvector with the vertical direction, and rotate the
first and second-order derivatives accordingly (following the transforms
derived in appendix E of Feldbrugge & van de Weygaert (2023)).

  4. 4.

Generate a realization of the initial Gaussian random field using the Hoffman-
Ribak algorithm for Gaussian random field generation with linear constraints
(Bertschinger, 1987; Hoffman & Ribak, 1991; van de Weygaert & Bertschinger,
1996).

  5. 5.

For each sampled initial field realization, two dark-matter only 2D $N$-body
simulations (Hidding, 2020) are run with $256\times 256$ particles in a
$25R_{s}\times 25R_{s}$ box, within an Einstein-de Sitter cosmology setting.
The regular simulation evolves the sampled initial conditions. The second
simulation concerns the evolution of a smoothed version of the initial
conditions, with a Gaussian kernel (27) of smoothing length $\sigma$, the
scale of the imposed constraints. The smoothed $N$-body simulation is used to
center and align the regular simulation.

For this class of potential field saddle point configurations, we generated an
ensemble of $1000$ simulations for length scales $\sigma=0.1,0.5,1$.

####  7.2.2 Realizations

The filament emanating from the primordial gravitational potential criterion
appears to represent a more genuine filamentary structure than the one
resulting from the primordial density field. This is clearly borne out by the
3 filament realizations shown in figure 16. In all cases, we recognize clear
filamentary structures aligned along the vertical axis. Nonetheless, there
exists a far larger variation than in the case of the caustic skeleton
filaments: some of the filaments are fully formed, while other structures
still reside in an embryonic stage and will only form later.

The median density field of the filamentary realizations, shown in figure 16,
is that of an oval-shaped configuration aligned with the vertical axis. Not
surprisingly, the length – and width – of the median filament increases with
the length scale $\sigma$, before widening at its tips.

The visual impression yielded by figure 16 and the median density fields are
supported by the analysis of the multi-stream character of the configurations.
In a significant fraction of the realizations, the saddle points lie in a
single-stream region. This tends to be either a filament that has not yet
formed or one located in a void region. Overall, the $3$-, $5$\- and
$7$-stream regions appear to represent a wide variety of vertically oriented
filaments. Occasionally, they may even end up in a cluster in a cluster in
Eulerian space.

##  8 The three-dimensional cosmic web

The present paper studies the formation of filaments in a two-dimensional
model of large structure formation. We propose a non-linear condition on the
primordial initial conditions that identifies the proto-filaments forming at a
set time and length scale in the cosmic web. In an upcoming paper, we will
extend this formalism to the three-dimensional cosmic web consisting of voids,
walls, filaments and clusters.

While the details of the three-dimensional caustic skeleton conditions are
described in the upcoming paper, here we briefly summarize how the two-
dimensional analysis generalizes to the three-dimensional setting.

In the three-dimensional cosmic web, the cusp caustic is associated with the
walls in the matter distribution. It is described by the manifold,

$\displaystyle A_{3}=\\{\bm{q}\in L\,|\,\lambda_{1}(\bm{q})\geq
1/b_{+}(t_{0})\,,\ \bm{v}_{1}(\bm{q})\cdot\nabla\lambda_{1}(\bm{q})=0\\}\,.$
(100)

Like the filaments of the two-dimensional cosmic web, it is natural to define
the center of the wall as a point on the wall that is always being stretched,
and – in the Zel’dovich approximation – is identified with the point of
maximum local stretch.

Formally, in Lagrangian space, the two-dimensional tangent space
$T_{\bm{q}_{c}}A_{3}$ of the cusp manifold is spanned by the orthogonal basis
$\bm{T}_{1}$ and $\bm{T}_{2}$. In Eulerian space, the stretching or
contraction of the wall is governed by the outer product
$\|(\nabla\bm{x}_{t}\bm{T}_{1})\times(\nabla\bm{x}_{t}\bm{T}_{2})\|$.

In the three-dimensional cosmic web, the caustic skeleton entails the
existence of two classes of filaments. The first one of these is the
swallowtail filaments,

$\displaystyle A_{4}=\\{$ $\displaystyle\bm{q}\in
L\,|\,\lambda_{1}(\bm{q})\geq 1/b_{+}(t_{0})\,,\
\bm{v}_{1}(\bm{q})\cdot\nabla\lambda_{1}(\bm{q})=0\,\ $
$\displaystyle\bm{v}_{1}(\bm{q})\cdot\nabla(\bm{v}_{1}(\bm{q})\cdot\nabla\lambda_{1}(\bm{q}))=0\\}\,,$
(101)

the second class those of the umbilic filaments,

$\displaystyle D_{4}=\\{$ $\displaystyle\bm{q}\in
L\,|\,\lambda_{1}(\bm{q})=\lambda_{2}(\bm{q})\geq 1/b_{+}(t_{0})\\}\,.$ (102)

By determining the typical orientation of the filament $\bm{T}$, the center of
the filament is the point that is always stretching – according to the
Zel’dovich approximation – is the location of maximum stretch. Its location
follows from the evaluation of the stretching rate
$\|\nabla\bm{x}_{t}\bm{T}\|$.

##  9 Summary and Conclusions

We propose an unambiguous definition of filaments based on their formation
histories rather than their morphology, in the two-dimensional setting.

The cosmic web can be characterized by its multi-stream nature. Mass elements
that have undergone shell-crossing form multi-stream regions where non-linear
gravitational collapse takes hold. These multi-stream regions result from the
folding of the dark matter sheet in phase-space and can be identified in terms
of caustics, famously classified by catastrophe theory. The two-dimensional
filament was previously linked to the cusp caustic (Arnol’d et al., 1982;
Arnol’d, 1982; Hidding et al., 2014). For the Zel’dovich approximation, this
cusp caustic is expressed in terms of the first eigenvalue field of the
primordial tidal tensor and its derivative in the direction of the
corresponding eigenvector field. Although the cusp curve neatly follows the
filamentary structure of the cosmic web, we show that it by itself does not
suffice to identify proto-filaments, as the cusp curve stretches in the
filament regions and tangles up in the cluster regions where the Zel’dovich
approximation fails to capture structure formation. For this reason, we
augment the cusp condition to find the centre of the filament, where the cusp
curve is maximally stretched. These additional conditions are expressed in
terms of the second eigenvalue field and its derivative in the direction of
its corresponding eigenvector field. This yields a family of proto-filaments
specified in terms of their orientation, formation time, and length scale.

Using constrained Gaussian random field theory (Feldbrugge et al., 2023) and
dark matter only $N$-body simulations (Hidding, 2020), we demonstrate that
this filament condition successfully identifies proto-filaments. We compare
the structures with previously proposed proto-filament based on the saddle
points of the primordial density and gravitational potential. The saddle point
of the gravitational potential does a better job than the saddle points of the
density perturbation. As it turns out, the saddle point of the primordial
density field is correlated to the filaments but the eigenvectors of the
Hessian are not a reliable indication of the orientation of the proto-
filament. The orientation of the saddle point of the primordial gravitational
potential leads to a filamentary structure, however, it does not identify the
time at which it is likely to form. The filament definition in terms of the
eigenvalue fields performs better than the definition in terms of the
gravitational potential. Note that the here proposed filament condition
includes the formation time.

We believe that this unambiguous definition will aid the systematic study of
filaments in the cosmic web and their relation to the embedded galaxies and
baryonic matter. In the future, we plan to use this definition to study the
recently identified filament spin (Wang et al., 2021; Xia et al., 2021),
galaxy alignment (Porciani et al., 2002; Aragón-Calvo et al., 2007b; Hahn et
al., 2007; Aragón Calvo, 2007; Schäfer, 2009; Hahn et al., 2009; Tempel &
Libeskind, 2013; Codis et al., 2012, 2015, 2018b; Ganeshaiah Veena et al.,
2018, 2019; López et al., 2021; López, 2024), and to the question in how far
the mass, stellar content and star formation activity are influenced (Cautun
et al., 2014; Kraljic et al., 2018; Kuutma et al., 2017; Laigle et al., 2018;
Hellwing et al., 2021; Parente et al., 2023).

The present discussion focuses on two-dimensional models of large-scale
structure formation, as this enables clear visualization. The work will be
extended to three-dimensional structure formation in a sequel paper. Moreover,
we focused on the caustic skeleton of the Zel’dovich approximation. In an
upcoming paper, we will extend the analysis to higher-order Lagrangian
perturbation theory.

## Acknowledgements

We thank Johan Hidding for providing the two-dimensional $N$-body simulation
code used in this project. We thank Benjamin Hertzsch Maé Rodriguez for
conversations about the physics underlying cosmic structure formation. We are
also grateful to Nynke Niezink for many discussions and for proofreading this
manuscript and for an incisive manuscript comment to Lucia van de Weijgaert.

The work is supported by the STFC Consolidated Grant ‘Particle Physics at the
Higgs Centre’ and by a Higgs Fellowship at the University of Edinburgh.

## Data Availability

No new data were generated or analyzed in support of this research.

## References

  * Abel et al. (2012) Abel T., Hahn O., Kaehler R., 2012, MNRAS, 427, 61
  * Adler (1981) Adler R. J., 1981, The Geometry of Random Fields. Society for Industrial and Applied Mathematics
  * Adler & Taylor (2009) Adler R., Taylor J., 2009, Random Fields and Geometry. Springer Monographs in Mathematics, Springer New York
  * Alpaslan et al. (2014) Alpaslan M., et al., 2014, MNRAS, 440, L106
  * Aragón Calvo (2007) Aragón Calvo M. A., 2007, PhD thesis, University of Groningen, Netherlands
  * Aragon-Calvo (2016) Aragon-Calvo M. A., 2016, MNRAS, 455, 438
  * Aragon-Calvo (2024) Aragon-Calvo M. A., 2024, MNRAS, 529, 74
  * Aragon-Calvo & Szalay (2013) Aragon-Calvo M. A., Szalay A. S., 2013, MNRAS, 428, 3409
  * Aragón-Calvo et al. (2007a) Aragón-Calvo M. A., Jones B. J. T., van de Weygaert R., van der Hulst J. M., 2007a, A&A, 474, 315
  * Aragón-Calvo et al. (2007b) Aragón-Calvo M. A., van de Weygaert R., Jones B. J. T., van der Hulst J. M., 2007b, ApJ, 655, L5
  * Aragón-Calvo et al. (2010a) Aragón-Calvo M. A., van de Weygaert R., Jones B. J. T., 2010a, MNRAS, 408, 2163
  * Aragón-Calvo et al. (2010b) Aragón-Calvo M. A., Platen E., van de Weygaert R., Szalay A. S., 2010b, ApJ, 723, 364
  * Aragon Calvo et al. (2019) Aragon Calvo M. A., Neyrinck M. C., Silk J., 2019, The Open Journal of Astrophysics, 2, 7
  * Arnol’d (1982) Arnol’d V. I., 1982, Trudy Seminar imeni G Petrovskogo, 8, 21
  * Arnold (1984) Arnold V. I., 1984, Catastrophe theory. Springer Berlin, Heidelberg
  * Arnol’d et al. (1982) Arnol’d V. I., Shandarin S. F., Zel’dovich I. B., 1982, Geophysical and Astrophysical Fluid Dynamics, 20, 111
  * Bardeen et al. (1983) Bardeen J. M., Steinhardt P. J., Turner M. S., 1983, Phys. Rev. D, 28, 679
  * Bardeen et al. (1986) Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, The Astrophysical Journal, 304, 15
  * Bermejo et al. (2024) Bermejo R., Wilding G., van de Weygaert R., Jones B. J. T., Vegter G., Efstathiou K., 2024, MNRAS, 529, 4325
  * Bertschinger (1987) Bertschinger E., 1987, ApJ, 323, L103
  * Bi et al. (2024) Bi D., Shlosman I., Romano-Díaz E., 2024, MNRAS, 527, 11095
  * Bond & Myers (1996) Bond J. R., Myers S. T., 1996, ApJS, 103, 1
  * Bond et al. (1996) Bond J. R., Kofman L., Pogosyan D., 1996, Nature, 380, 603
  * Bonnaire et al. (2022) Bonnaire T., Aghanim N., Kuruvilla J., Decelle A., 2022, A&A, 661, A146
  * Cadiou et al. (2021) Cadiou C., Pontzen A., Peiris H. V., Lucie-Smith L., 2021, MNRAS, 508, 1189
  * Canducci et al. (2022) Canducci M., et al., 2022, Astronomy and Computing, 41, 100658
  * Cautun et al. (2013) Cautun M., van de Weygaert R., Jones B. J. T., 2013, MNRAS, 429, 1286
  * Cautun et al. (2014) Cautun M., van de Weygaert R., Jones B. J. T., Frenk C. S., 2014, MNRAS, 441, 2923
  * Codis et al. (2012) Codis S., Pichon C., Devriendt J., Slyz A., Pogosyan D., Dubois Y., Sousbie T., 2012, MNRAS, 427, 3320
  * Codis et al. (2015) Codis S., Pichon C., Pogosyan D., 2015, MNRAS, 452, 3369
  * Codis et al. (2018a) Codis S., Pogosyan D., Pichon C., 2018a, MNRAS, 479, 973
  * Codis et al. (2018b) Codis S., Jindal A., Chisari N. E., Vibert D., Dubois Y., Pichon C., Devriendt J., 2018b, MNRAS, 481, 4753
  * Colberg et al. (2005) Colberg J. M., Krughoff K. S., Connolly A. J., 2005, MNRAS, 359, 272
  * Colless et al. (2003) Colless M., et al., 2003, arXiv e-prints, pp astro–ph/0306581
  * Colombi (2015) Colombi S., 2015, MNRAS, 446, 2902
  * Courtois et al. (2013) Courtois H. M., Pomarède D., Tully R. B., Hoffman Y., Courtois D., 2013, AJ, 146, 69
  * de Lapparent et al. (1986) de Lapparent V., Geller M. J., Huchra J. P., 1986, ApJ, 302, L1
  * Dhawalikar & Paranjape (2024) Dhawalikar S., Paranjape A., 2024, arXiv e-prints, p. arXiv:2402.18669
  * Doroshkevich (1970) Doroshkevich A. G., 1970, Astrophysics, 6, 320
  * Eisenstein et al. (1997) Eisenstein D. J., Loeb A., Turner E. L., 1997, ApJ, 475, 421
  * Elek et al. (2021) Elek O., Burchett J. N., Prochaska J. X., Forbes A. G., 2021, IEEE Transactions on Visualization and Computer Graphics, 27, 806
  * Elek et al. (2022) Elek O., Burchett J. N., Prochaska J. X., Forbes A. G., 2022, Artificial Life, 28, 22
  * Feldbrugge (2024) Feldbrugge J., 2024, arXiv e-prints, p. arXiv:2402.16234
  * Feldbrugge & van de Weygaert (2023) Feldbrugge J., van de Weygaert R., 2023, J. Cosmology Astropart. Phys., 2023, 058
  * Feldbrugge et al. (2016) Feldbrugge J. L., Hidding J., van de Weygaert R., 2016, in van de Weygaert R., Shandarin S., Saar E., Einasto J., eds, Vol. 308, The Zel’dovich Universe: Genesis and Growth of the Cosmic Web. pp 107–114, doi:10.1017/S1743921316009704
  * Feldbrugge et al. (2018) Feldbrugge J., van de Weygaert R., Hidding J., Feldbrugge J., 2018, J. Cosmology Astropart. Phys., 2018, 027
  * Feldbrugge et al. (2019) Feldbrugge J., van Engelen M., van de Weygaert R., Pranav P., Vegter G., 2019, J. Cosmology Astropart. Phys., 2019, 052
  * Feldbrugge et al. (2023) Feldbrugge J., Yan Y., van de Weygaert R., 2023, MNRAS, 526, 5031
  * Forero-Romero et al. (2009) Forero-Romero J. E., Hoffman Y., Gottlöber S., Klypin A., Yepes G., 2009, MNRAS, 396, 1815
  * Galárraga-Espinosa et al. (2024) Galárraga-Espinosa D., et al., 2024, A&A, 684, A63
  * Ganeshaiah Veena (2020) Ganeshaiah Veena P., 2020, PhD thesis, University of Groningen, Netherlands
  * Ganeshaiah Veena et al. (2018) Ganeshaiah Veena P., Cautun M., van de Weygaert R., Tempel E., Jones B. J. T., Rieder S., Frenk C. S., 2018, MNRAS, 481, 414
  * Ganeshaiah Veena et al. (2019) Ganeshaiah Veena P., Cautun M., Tempel E., van de Weygaert R., Frenk C. S., 2019, MNRAS, 487, 1607
  * Gingold & Monaghan (1977) Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375
  * Giovanelli et al. (1986) Giovanelli R., Haynes M. P., Chincarini G. L., 1986, ApJ, 300, 77
  * Guth & Pi (1982) Guth A. H., Pi S. Y., 1982, Phys. Rev. Lett., 49, 1110
  * Hahn et al. (2007) Hahn O., Porciani C., Carollo C. M., Dekel A., 2007, MNRAS, 375, 489
  * Hahn et al. (2009) Hahn O., Porciani C., Dekel A., Carollo C. M., 2009, MNRAS, 398, 1742
  * Hahn et al. (2010) Hahn O., Teyssier R., Carollo C. M., 2010, MNRAS, 405, 274
  * Harlow et al. (1955) Harlow F., Evans M., Richtmyer R., 1955, A Machine Calculation Method for Hydrodynamic Problems. LAMS (Los Alamos Scientific Laboratory), Los Alamos Scientific Laboratory of the University of California
  * Harrison (1970) Harrison E. R., 1970, Phys. Rev. D, 1, 2726
  * Hellwing et al. (2021) Hellwing W. A., Cautun M., van de Weygaert R., Jones B. T., 2021, Phys. Rev. D, 103, 063517
  * Hidding (2020) Hidding J., 2020, jhidding/nbody2d: 2d PM n-body code, doi:10.5281/zenodo.4158731, https://doi.org/10.5281/zenodo.4158731
  * Hidding et al. (2014) Hidding J., Shandarin S. F., van de Weygaert R., 2014, MNRAS, 437, 3442
  * Hoffman & Ribak (1991) Hoffman Y., Ribak E., 1991, ApJ, 380, L5
  * Hoffman et al. (2012) Hoffman Y., Metuki O., Yepes G., Gottlöber S., Forero-Romero J. E., Libeskind N. I., Knebe A., 2012, MNRAS, 425, 2049
  * Huchra et al. (2012) Huchra J. P., et al., 2012, ApJS, 199, 26
  * Icke (1973) Icke V., 1973, A&A, 27, 1
  * Icke & van de Weygaert (1987) Icke V., van de Weygaert R., 1987, A&A, 184, 16
  * Jaber et al. (2024) Jaber M., Peper M., Hellwing W. A., Aragón-Calvo M. A., Valenzuela O., 2024, MNRAS, 527, 4087
  * Joeveer & Einasto (1978) Joeveer M., Einasto J., 1978, in Longair M. S., Einasto J., eds, 1 Vol. 79, Large Scale Structures in the Universe. p. 241
  * Kraljic et al. (2018) Kraljic K., et al., 2018, MNRAS, 474, 547
  * Kuchner et al. (2020) Kuchner U., et al., 2020, MNRAS, 494, 5473
  * Kuchner et al. (2021) Kuchner U., et al., 2021, MNRAS, 503, 2065
  * Kugul (2020) Kugul R., 2020, Master’s thesis, University of Groningen, Netherlands
  * Kuutma et al. (2017) Kuutma T., Tamm A., Tempel E., 2017, A&A, 600, L6
  * Laigle et al. (2018) Laigle C., et al., 2018, MNRAS, 474, 5437
  * Libeskind et al. (2018) Libeskind N. I., et al., 2018, MNRAS, 473, 1195
  * Longuet-Higgins (1957) Longuet-Higgins M. S., 1957, Philosophical Transactions of the Royal Society of London Series A, 250, 157
  * López (2024) López P., 2024, PASP, 136, 037001
  * López et al. (2021) López P., Cautun M., Paz D., Merchán M., van de Weygaert R., 2021, MNRAS, 502, 5528
  * Lu et al. (2024) Lu Y. S., Mandelker N., Oh S. P., Dekel A., van den Bosch F. C., Springel V., Nagai D., van de Voort F., 2024, MNRAS, 527, 11256
  * Mukhanov & Chibisov (1981) Mukhanov V. F., Chibisov G. V., 1981, Soviet Journal of Experimental and Theoretical Physics Letters, 33, 532
  * Neyrinck (2012) Neyrinck M. C., 2012, MNRAS, 427, 494