# Evolution of cosmic filaments in the MTNG simulation*

Daniela Galárraga-Espinosa<fnref n="1" />, Carentin Cadiou<fnref n="2" />, Céline Gouin<fnref n="3" />, Simon D. M. White<fnref n="1" />, Volker Springel<fnref n="1" />, Rüdiger Pakmor<fnref n="1" />, Boryana Hadzhiyska<fnref n="4" />,<fnref n="5" />, Sownak Bose<fnref n="6" />, Fulvio Ferlito<fnref n="1" />, Lars Hernquist<fnref n="7" />, Rahul Kannan<fnref n="8" />, Monica Barrera<fnref n="1" />, Ana Maria Delgado<fnref n="7" />, and César Hernández-Aguayo<fnref n="1" />

Received 15 September 2023 / Accepted 10 January 2024

## ABSTRACT

We present a study of the evolution of cosmic filaments across redshift with an emphasis on some important properties: filament lengths, growth rates, and radial profiles of galaxy densities. Following an observation-driven approach, we built cosmic filament catalogues at $z = 0, 1, 2, 3$, and 4 from the galaxy distributions of the large hydro-dynamical run of the MillenniumTNG project. We employed the extensively used DisPerSE cosmic web finder code, for which we provide a user-friendly guide, including the details of a physics-driven calibration procedure, with the hope of helping future users. We performed the first statistical measurements of the evolution of connectivity in a large-scale simulation, finding that the connectivity of cosmic nodes (defined as the number of filaments attached) globally decreases from early to late times. The study of cosmic filaments in proper coordinates reveals that filaments grow in length and radial extent, as expected from large-scale structures in an expanding Universe. But the most interesting results arise once the Hubble flow is factored out. We find remarkably stable comoving filament length functions and over-density profiles, showing only little evolution of the total population of filaments in the past ~12.25 Gyr. However, by tracking the spatial evolution of individual structures, we demonstrate that filaments of different lengths actually follow different evolutionary paths. While short filaments preferentially contract, long filaments expand along their longitudinal direction with growth rates that are the highest in the early, matter-dominated Universe. Filament diversity at a fixed redshift is also shown by the different (~$5\sigma$) density values between the shortest and longest filaments. Our results hint that cosmic filaments can be used as additional probes for dark energy, but further theoretical work is still needed.

**Key words.** methods: numerical – methods: statistical – galaxies: clusters: general – cosmology: theory – large-scale structure of Universe

## 1. Introduction

The first observations of galaxy distributions around the Perseus cluster (Jöeveer et al. 1978), around the Coma/A1367 supercluster (Gregory & Thompson 1978), and in the Center for Astrophysics (CfA) galaxy survey (de Lapparent et al. 1986) showed that on the largest scales, matter in the Universe is organised in clusters, filaments, walls, and voids. These structures form the cosmic web (White et al. 1987; Bond et al. 1996), a multi-scale network of dark matter (DM) and gas, which emerged under the action of gravity from the anisotropic collapse of initial fluctuations of the density field (Zel'dovich 1970). Following its discovery, the cosmic web as traced by galaxies has been observed in a variety of galaxy surveys, such as the Sloan Digital Sky Survey (SDSS, York et al. 2000), the two degree Field Galaxy Redshift Survey (2dFGRS, Colless et al. 2003), Cosmic

Evolution and COSMOS2015 (Scoville et al. 2007; Laigle et al. 2016), the 6dF Galaxy Survey (6dFGS, Jones et al. 2009), the Galaxy and Mass Assembly survey (GAMA, Driver et al. 2011), the VIMOS VLT deep survey (VVDS, Le Fèvre et al. 2005), VIPERS (Guzzo et al. 2014), WISExSuperCOSMOS (Bilicki et al. 2016), and the SAMI survey (Bryant et al. 2015), each with a different resolution and statistical power.

While the nodes of the web, traced by the most massive clusters of galaxies, have been intensively characterised (e.g. Navarro et al. 1997; Nagai et al. 2007; Arnaud et al. 2010; Baxter et al. 2017; Bartalucci et al. 2017; Pintos-Castro et al. 2019; Ghirardini et al. 2019, among many others), cosmic filaments are still only poorly understood. This is principally due to their low matter densities which (combined with their complex morphology) result in a fainter signal with respect to that of the cosmic nodes. For example, the milestone analysis of Cautun et al. (2014) has shown that the probability distribution function (PDF) of the DM over-density field in filaments is, on average, two orders of magnitude lower than that in nodes, although the exact number is definition dependent (e.g.

<footnote n="1">¹ Max-Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, 85741 Garching, Germany
e-mail: danigaes@mpa-garching.mpg.de</footnote>

<footnote n="2">² Lund Observatory, Division of Astrophysics, Department of Physics, Lund University, Box 43, 221 00 Lund, Sweden</footnote>

<footnote n="3">³ Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France</footnote>

<footnote n="4">⁴ Physics Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA</footnote>

<footnote n="5">⁵ Miller Institute for Basic Research in Science, University of California, Berkeley, CA 94720, USA</footnote>

<footnote n="6">⁶ Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="7">⁷ Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="8">⁸ Department of Physics and Astronomy, York University, 4700 Keele Street, Toronto, ON M3J 1P3, Canada</footnote>

<footnote>* The MillenniumTNG simulations will be publicly released at https://www.mtng-project.org in the near future. The cosmic filament catalogues constructed and analysed in this work will be made available to the community at the same webpage.</footnote>
Busch & White 2020). Observing cosmic filaments in galaxy surveys is thus extremely challenging.

In addition to these observational difficulties, the lack of an unequivocal definition of cosmic filaments complicates the study of these structures even further. As illustrated by the detailed review of Libeskind et al. (2018), a variety of methods have been developed to identify filaments in observational and simulated datasets, each adopting a different approach. For example, some algorithms define filaments based on the topology of the matter density field (Sousbie et al. 2011; Sousbie 2011), on the phase-space caustics of DM (Neyrinck 2012), or through the relative strength between eigenvalues of the Hessian matrix of the density field (e.g. Aragón-Calvo et al. 2007, 2010a; Cautun et al. 2014; Pfeifer et al. 2022). Filaments can also be identified based on statistical representations of stochastic processes (e.g. Stoica et al. 2007; Tempel et al. 2016), using graph theory (e.g. Pereyra et al. 2020), regularised minimum spanning trees (e.g. Bonnaire et al. 2020), or via machine learning classification tools (e.g. Buncher & Carrasco Kind 2020; Inoue et al. 2022; Awad et al. 2023). This diversity of filament finders underlines the high complexity involved in the study and physical characterisation of these structures.

Importantly, studies of numerical simulations have shown that cosmic filaments today contain more than 50% of the matter in the Universe (Cautun et al. 2014; Ganeshaiah Veena et al. 2019). The understanding of matter at the largest scales is therefore inevitably linked to that of these structures. The last decade has thus witnessed an outstanding number of studies aiming to characterise the properties of cosmic filaments themselves at $z \sim 0$ (e.g. Colberg et al. 2005; Gheller et al. 2015, 2016; Ho et al. 2018; Connor et al. 2018; Galárraga-Espinosa et al. 2020, 2022; Rost et al. 2020, 2021; Pereyra et al. 2020; Malavasi et al. 2020a; Wang et al. 2021; Zhu et al. 2021; Vernstrom et al. 2021; Zakharova et al. 2023), the properties of the haloes (e.g. Aragón-Calvo et al. 2007; Hahn et al. 2007; Musso et al. 2018; Ganeshaiah Veena et al. 2018, 2021; Jhee et al. 2022) and galaxies living in these structures (e.g. Pandey & Bharadwaj 2006; Tempel et al. 2014; Chen et al. 2017; Kuutma et al. 2017; Malavasi et al. 2017, 2022; Bonjean et al. 2018, 2020; Laigle et al. 2018; Kraljic et al. 2018, 2019, 2020a; Sarron et al. 2019; Liao & Gao 2019; Ganeshaiah Veena et al. 2019; Santiago-Bautista et al. 2020; Welker et al. 2020; Singh et al. 2020; Gouin et al. 2020; Song et al. 2021; Kotecha et al. 2022; Kuchner et al. 2022), and those of gas around cosmic filaments (e.g. Klar & Mücket 2012; Nevalainen et al. 2015; Gheller & Vazza 2019; de Graaff et al. 2019; Tanimura et al. 2019, 2020a,b, 2022; Tuominen et al. 2021; Martizzi et al. 2019; Galárraga-Espinosa et al. 2021; Biffi et al. 2022; Holt et al. 2022; Gouin et al. 2022; Erciyes et al. 2023; Lu et al. 2024), both in numerical simulations and observations. From these studies, we have learned about the broad diversity of filaments in the late Universe: from hot and dense bridges of matter between clusters, to longer, cooler (but still $T \gtrsim 10^5$ K on average) cosmic filaments living in less dense environments, for example close to voids. We have understood that galaxies are pre-processed in cosmic filaments before falling into clusters, and that halo properties are affected by these large-scale environments, in particular their spins. We have reached the consensus that cosmic filaments are the hosts of most of the “missing baryons” (Cen & Ostriker 2006; Shull et al. 2012) in the late Universe, with the warm-hot intergalactic medium (WHIM) absolutely dominating the gas content of these structures in simulations, although this phase is still partially elusive in today’s observations.

In this context, this work aims at setting the stage for pushing those studies to higher redshifts, with the goal of building a global picture of the evolution of cosmic filaments and of matter within them. We use the outputs of the large hydro-dynamical MillenniumTNG simulation (e.g. Hernández-Aguayo et al. 2023; Pakmor et al. 2023) together with the filament finder code DisPerSE (Sousbie 2011; Sousbie et al. 2011), both introduced in Sect. 2, to construct cosmic filament catalogues at $z = 0, 1, 2, 3$ and 4. The extraction procedure of these catalogues is presented in Sect. 3, in which we also detail our robustness checks performed for the filament catalogues. The latter organically led to a first exploration of the evolution with redshift of the number of connections of massive haloes to the filaments of the web, shown in Sect. 4. Thanks to the computation of filament evolutionary tracks in Sect. 5, we analyse and interpret crucial properties such as the evolution of filament length functions and growth rates. In Sect. 6, we present the evolution of galaxy density profiles around cosmic filaments. Finally, we summarise and further discuss our results in Sect. 7.

The cosmic web being a multi-scale structure, it is important to point out that this work focuses on the filaments at the largest scales of the Universe (i.e. the cosmic filaments). These connect galaxy clusters and host galaxies and small groups, as opposed to the smaller-scale (~kiloparsec) filaments feeding individual galaxies with fuel for star-formation (the “cold stremes” e.g. Ramsøy et al. 2021; Galárraga-Espinosa et al. 2023; Lu et al. 2024). Throughout this paper, we adopt the values of the cosmological parameters given by Planck Collaboration XIII (2016), that are, $\Omega_{\Lambda,0} = 0.6911$, $\Omega_{m,0} = 0.3089$, $\Omega_{b,0} = 0.0486$, $\sigma_8 = 0.8159$, $n_s = 0.9667$, and $h = 0.6774$. All the error bars correspond to the errors on the mean values, derived from bootstrap resampling.

## 2. Data and methods

### 2.1. The MTNG simulations

We analyse the outputs of the full physics run of the MillenniumTNG (MTNG) project<fnref n="1" />. This hydrodynamical simulation combines the successful galaxy formation model of the IllustrisTNG suite (Pillepich et al. 2018; Nelson et al. 2019) with the large volume ($500\,\text{cMpc}\,h^{-1}$)³ or ($738.1\,\text{cMpc}\,h^{-1}$)³ of the iconic Millennium simulation (Springel et al. 2005). Compared to the widely used TNG300 simulation, the MTNG hydrodynamical run is ~15 times larger in volume and has a mass resolution that is only reduced by a factor of 2.8. Indeed, with $4320^3$ DM particles and the same number of gas cells, the DM resolution of this simulation is $m_{\text{DM}} = 1.12 \times 10^8\,\text{M}_\odot\,h^{-1}$, and the target baryonic mass is $m_{\text{target}} = 2.0 \times 10^7\,\text{M}_\odot\,h^{-1}$. The MTNG simulation is thus ideal to analyse the interplay between the large-scales of the cosmic web and the non-linear scales of galaxy formation and evolution. In particular, it is one of the best simulations up to date to perform statistical studies of cosmic filaments, another one being the Flamingo project (Schaye et al. 2023). Further information on the MTNG project including details on the model and available outputs are provided in the introductory papers of Hernández-Aguayo et al. (2023) and Pakmor et al. (2023), the latter with particular emphasis on the hydrodynamical run used in this work. The other introductory papers of the MTNG project cover a wide range of astrophysical and cosmological topics, such as the description of the lightcones (Barrera et al. 2023), clustering analyses (Bose et al. 2023), inference of cosmological

<footnote n="1">¹ https://www.mtng-project.org/</footnote>
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Sketch representing the large-scale density field and its topology using mountains as an analogy. The black lines following the ridges of the field represent the cosmic filaments. The red and green circles show, respectively, the positions of the CPmax, or peaks, and the (filament-type) saddles of the field. These are the endpoints of filaments within the DisPerSE framework. The blue circles locate the wall-type saddles of the field. Minimum density critical points are not represented.</figcaption></figure>

parameters (Contreras et al. 2023), galaxy intrinsic alignments (Delgado et al. 2023), weak-lensing convergence maps (Ferlito et al. 2023), halo occupation models (Hadzhiyska et al. 2023a,b), and high redshift ($z > 8$) galaxies (Kannan et al. 2023).

The detection and analysis of cosmic filaments in this work makes extensive use of the catalogues of galaxies and groups found respectively at the centres of the SUBFIND-HBT subhaloes and Friend-of-Friends (FoF) haloes, identified on-the-fly following Springel et al. (2021). We perform our study at redshifts $z = 0, 1, 2, 3,$ and $4$. These output times were chosen because their corresponding snapshots offer the full particle data, thus enabling future studies of, for example, gas properties around the filaments introduced in this work. We note that we made the conservative choice of stopping the analysis at $z = 4$, as galaxy properties at higher redshifts are expected to be affected by the non-uniform UV background due to the epoch of reionisation. Finally, this work analyses the distributions of objects in real space and so does not account for peculiar velocity effects.

## 2.2. The DisPerSE filament finder: A user's guide

In the following, we concisely describe the steps required for the extraction of filaments from galaxy catalogues using the DisPerSE (Sousbie 2011; Sousbie et al. 2011) cosmic web finder. With the goal of providing a complete procedure that could be helpful for future DisPerSE users, we focus explicitly on the practical details related to the application of the code. For the theoretical background at the core of DisPerSE, we refer to the presentation papers (Sousbie 2011; Sousbie et al. 2011) and website<fnref n="2" />. We recall that this work adopts an observation-driven approach, in which the large-scale density field is traced by the positions of galaxies (more details in Sect. 3.1).

The first step for extracting DisPerSE filaments from a distribution of galaxies is to estimate the underlying (continuous) density field. This is done by the Delaunay tessellation field estimator (hereafter DTFE, Schaap & van de Weygaert 2000; van de Weygaert & Schaap 2009) implemented in the

delaunay_3d function of DisPerSE. This function provides an estimate of the density at the position of each tracer of the input dataset, which are located by definition at vertices of the Delaunay tessellation. The user can then decide to “smooth” the resulting DTFE densities (hereafter $\rho_{DTFE}$) by applying the netconv function to the tessellation. The latter replaces the $\rho_{DTFE}$ value of a given vertex by the average density of the surrounding vertices in the tessellation, thus effectively reducing the contamination by shot noise in the density values (e.g. Malavasi et al. 2020a,b). The number of iterations of this averaging procedure can be decided by the user (using the N1 parameter, see below).

The identification of the filaments is then performed by the mse function (the main core of the DisPerSE code). This step identifies the critical points in the (possibly smoothed) Delaunay density field, which are locations where the density gradient vanishes. Filaments are then defined as the sets of segments connecting maximum-density critical points (hereafter CPmax, or peaks) to saddles, following the ridges of the density field. The most relevant topological features identified by DisPerSE are illustrated in the drawing of Fig. 1. In that sketch, the definition of the filaments (black lines) can be clearly understood. Importantly, the significance level of the detected filaments with respect to the noise can be controlled by adjusting the persistence ratio of the CPmax-saddle critical point pairs (i.e. of the filament’s ending points). This parameter, hereafter N2, is defined as the ratio of the density of the two critical points in the pair (Sousbie et al. 2011), and needs to be rigorously calibrated as we show in the following section.

Lastly, the skelconv function allows a conversion of the output of mse from binary to human readable format. The user can also choose to soften the angles between contiguous filament segments, exclusively for aesthetic reasons, without changing the position of the filaments’ sampling points. Reading the final output of DisPerSE needs additional post-processing codes. Those used in this paper are the same as in Garaldi et al. (2023) and are publicly available<fnref n="3" />.

The different DisPerSE keywords needed for the extraction of filaments are also explicitly provided in the following. These are -btype periodic, -smoothData field_value N1, -forceLoops -robustness -manifolds -upSkl -nsig N2, and -to NDskl_ascii -smooth 1, respectively, for the delaunay_3d, netconv, mse, and skelconv functions. Here N1 and N2 represent the only two parameters of the code, corresponding to the smoothing value of the Delaunay densities and the persistence ratio of the critical point pairs. A way of calibrating these parameters is introduced in Sect. 3.2. Finally, we mention that in this work we followed additional procedures in order to clean the DisPerSE outputs. We discarded filaments with zero lengths, which appear when the two ending points are located at the same position, and CPmax when associated to only one filament (together with the corresponding filament). The latter was motivated by the prior that nodes connected to a single filament are not expected at the scales of the cosmic web (Codis et al. 2018).

## 3. Extraction of cosmic filament catalogues at different redshifts

This section presents in detail the building of the cosmic web filament catalogues using DisPerSE on galaxy distributions. We first address the selection of the tracers of the cosmic skeletons at different redshifts in Sect. 3.1. Then, in Sect. 3.2, we explain the

<footnote n="2">² http://www2.iap.fr/users/sousbie/web/html/indexd41d.html</footnote>

<footnote n="3">³ https://github.com/EGaraldi/disperse_output_reader/</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Stellar mass functions of the $N_{\text{gal}}$ most massive subhaloes at redshifts $z = 0, 1, 2, 3$, and 4. These objects are used to trace the cosmic skeleton. The vertical dotted lines show the mean values, i.e. $\log_{10}(M_*/M_{\odot}) = 9.85, 9.70, 9.40, 9.01$, and 8.61, respectively from $z=0$ to 4.</figcaption></figure>

calibration procedure of DisPerSE parameters at $z=0$. Finally, the extraction of the higher-$z$ filament catalogues is presented in Sect. 3.3.

### 3.1. Choice of tracers

At redshift $z=0$, we select from the subhalo catalogue of the MTNG simulation all the objects with stellar masses larger than $10^9 M_{\odot}$. This lower limit is chosen following the observational limits of Brinchmann et al. (2004) and Taylor et al. (2011), and ensures a selection of subhaloes with masses large enough to trace the positions of the large-scale cosmic filaments (as in e.g. Galárraga-Espinosa et al. 2020). The number of objects at $z=0$ satisfying this criterion is $N_{\text{gal}} = 3\,083\,441$, resulting in a number density of $n_{\text{gal}} = 7.7 \times 10^{-3} \, \text{Mpc}^{-3}$ tracers in the total simulation volume. This is compatible with the expected number density of tracers in state-of-the-art galaxy surveys, such as DESI BGS at low redshift (Hahn et al. 2023). The stellar mass function of these subhaloes is presented by the dark red curve in Fig. 2.

The choice of using galaxies as tracers of the density field is entirely driven by observations, where it is relatively accessible to detect the cosmic skeleton from the distribution of galaxies (as done in surveys such as VIPERS, COSMOS2015, SDSS, and SAMI, respectively in Malavasi et al. 2017, 2020a; Laigle et al. 2018; Carrón Duque et al. 2022; Welker et al. 2020). Compared to skeletons extracted from the underlying DM distribution, this choice has a negligible impact on the positions of cosmic filaments (see Appendix C of Laigle et al. 2018), but can play a role in the number of recovered filaments, as expected from biased tracers (Zakharova et al. 2023).

At higher redshifts, the large-scale cosmic web is also traced by the most massive galaxies at a given cosmic time. We build the dataset of tracers for detecting the filaments at $z=1, 2, 3$, and 4, by selecting the $N_{\text{gal}}$ most massive subhaloes (in stellar mass) in the respective SUBFIND catalogues. The same number of subhaloes as for $z=0$ is adopted precisely with the goal of working at constant (comoving) density of tracers across cosmic time. This ensures that the detected high-$z$ filaments belong to the same spatial scales as the $z=0$ cosmic filaments.

The cosmic web being a multi-scale structure where smaller-scale objects can be found inside larger ones (e.g. galaxies connected by kiloparsecscale filaments, or tendrils inside galaxy clusters), the method adopted in this work enables the extraction and comparison between different filament catalogues without the contamination of structures at smaller-scales. One could also think about using the first progenitors of the $\log_{10}(M_*/M_{\odot}) \ge 9$ galaxies at $z=0$ as tracers of the higher redshift cosmic skeletons. The first reason against using them is that it would be impossible to perform the same analysis on observational data. The second is that, due to the complex pathways of mass assembly, the progenitors of massive galaxies at $z=0$ are not always the most massive objects at higher redshift. They are thus poor tracers of structure at the largest-scales, as shown in Appendix A.

The filamentary pattern is clearly present in the distribution of the most massive galaxies, visualised in the 15 cMpc thick slices of the left panels of Fig. 3. With increasing redshift from top to bottom, we can appreciate how the filamentary structure of the cosmic web is already apparent at $z=4$, and more sharply delineated at $z=0$. Finally, the stellar mass functions of the selected subhaloes at all redshifts are shown in Fig. 2.

### 3.2. Calibration of DisPerSE at $z=0$

As introduced in Sect. 2.2, the DisPerSE filament finder possesses two main parameters that need to be carefully calibrated. Starting from the galaxy distribution at $z=0$ presented above, we perform eight different extractions of the cosmic skeleton, each with a different combination of parameters, that is, with or without smoothing of the $\rho_{\text{DTFE}}$ distribution ($N1=0$ or 1), and persistence ratios $N2$ varying from 1 to $4\sigma$. The best configuration of DisPerSE parameters is the one which recovers the most salient features of the large-scale distribution while minimising the effects of noise (i.e. respectively maximising the completeness and purity of the skeleton). Nevertheless, testing the robustness of filaments is far from being a trivial task, since a “true” physical reference of the positions of filament spines can not be easily determined in simulations. To bypass this problem, one could think of applying different filament finders to the same field (e.g. those of Cautun et al. 2013; Tempel et al. 2016; Bonnaire et al. 2020; Pfeifer et al. 2022) and comparing the corresponding outputs. However, this exercise would be a characterisation of the similarities and differences between the different methods (as done by Libeskind et al. 2018), but would not provide conclusive results regarding the accuracy of the skeletons in detecting the actual filaments.

In this work, we assess the robustness of the skeletons by comparing the positions of the densest ending points of the filaments (the CPmax, which by construction correspond to the topological nodes) to those of the massive Friend-of-Friends (FoF) haloes of the simulation. The latter are classified into the three following mass bins: $M_{200c} \in 5 \times [10^{11}, 10^{12}]$, $M_{200c} \in 5 \times [10^{12}, 10^{13}]$, and $M_{200c} \ge 5 \times 10^{13} M_{\odot} h^{-1}$, with respectively 643 329, 70 392, and 6135 objects in the MTNG snapshot at $z=0$. Hereafter, we refer to them as galactic haloes, groups, and clusters. The full results of the matching between CPmax and FoF haloes can be found in Appendix B, together with further details on each skeleton. In the following, we discuss only the most relevant information of this matching.

The first panel of Fig. 4 shows, for each of the eight skeletons, the fraction of CPmax found within $R_{200c}$ of any of the massive FoF haloes defined above (i.e. the matched CPmax, in dashed lines with squares), and conversely, the fraction of clusters hosting a CPmax inside their $R_{200c}$ spheres (in dotted lines
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Slices of thickness 15 cMpc showing the galaxy distribution (left) and the associated filaments (right) at different redshifts of the MTNG simulation.</figcaption></figure>
<figure><img src="image_4.png" /></figure>

**Fig. 4.** Calibration of DisPerSE parameters at $z = 0$. Left: results of the matching between CPmax and massive haloes for eight different DisPerSE outputs, each with a different parameterisation (see main text). Blue and red colours show, respectively, the results for smoothed (S) and not smoothed (noS) density fields. Right: filament length distributions of the two best skeleton candidates.

with crosses). These fractions are maximised by the smoothed $2\sigma$ skeleton (S, $2\sigma$, at the intersection point in blue), and the $4\sigma$ skeleton with no smoothing (noS, $4\sigma$, in red). The other filament extractions are either too impure with too few matched CPmax caused by the lack of smoothing and low persistence ratio, or incomplete with excessively low fractions of clusters identified as CPmax because of overly high persistence values. The right panel of Fig. 4 presents the filament length distributions of the two best candidates. We notice that the one corresponding to the “noS, $4\sigma$” skeleton presents a contribution of extremely short structures with lengths as small as 0.03 Mpc, which are likely to be originating from the noise. Therefore, we conclude that, given the $n_{\text{gal}}$ density of tracers, the best possible extraction of large-scale cosmic filaments in MTNG is given by the smoothed-$2\sigma$ DisPerSE skeleton.

### 3.3. Higher redshift cosmic skeletons

We choose to detect the high-$z$ skeletons using the same DisPerSE parameterisation as for the $z=0$ case, that is, smoothing $\rho_{\text{DTFE}}$ and applying a $2\sigma$ persistence ratio between critical point pairs. This is done in order to retrieve structures that belong to the same scales and have similar density contrasts as the $z=0$ filaments. A visualisation of the resulting filaments at different redshifts can be appreciated in the right panels of Fig. 3.

Similarly to the $z=0$ case, we check the robustness of the high-$z$ skeletons by analysing the matching between the topological nodes (the CPmax) and the FoF haloes of the corresponding MTNG snapshots. We find that at all redshifts the most massive haloes of the simulation are always identified as CPmax, and are thus, by construction, at the topological nodes of our skeletons. This is shown by the nice overlap towards the high mass values between the mass functions of the nodes, $\Phi_{\text{node}}$ (in thick lines), and those of the FoF haloes (thin lines) presented in Fig. 5. To be clear, we here define the mass of a given node as that of its matched FoF halo, following the procedure described in the previous section.

The number of filaments and CPmax detected at each redshift are reported in Table 1. These decrease with redshift, as expected from the lower density contrasts in the early Universe compared to late times, thus resulting in fewer critical point pairs satisfying the fixed $2\sigma$ persistence ratio. For similar reasons, the fraction of filaments connected to massive FoF haloes also mildly decreases with redshift (from 91.5% at $z=0$ to 77.3% at $z=4$). Precisely, the more homogeneous density fields at high redshifts lead to more diluted density peaks that are less likely to be identified as nodes of the cosmic skeleton. We have checked that the disconnected filaments (i.e. the minority of filaments whose nodes are not matched with massive FoF haloes) do not have particular features in their length distributions, density

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

**Fig. 5.** Mass function of the nodes connected to the filaments at different redshifts (thick curves). The vertical marks show the mean masses, which are $\log_{10}(\bar{M}_{200c}/M_{\odot}) = 13.66, 13.30, 12.90, 12.47$, and $12.06$, respectively, from $z=0$ to 4. For comparison, the thin curves give the mass functions of all the FoF haloes of the MTNG simulation with a stellar mass content of at least $10^8 M_{\odot}$.

profiles, or position in the simulation box: they are just structures that connect haloes of low virial masses. Because of the cosmic scales studied in this work we choose to discard these few structures, so that the analyses of the following sections are performed on the connected filaments having more than one segment. This yields the final number of filaments reported in Table 1.

## 4. Connectivity of massive haloes across redshift

The cosmic filaments of the catalogues built in the previous section are all connected to massive FoF haloes (whose mass distributions were presented in Fig. 5). In this section, we perform a first analysis of the number of connections of these haloes, at the nodes of the cosmic web, with the detected cosmic filaments. Hereafter, we refer to this number as the connectivity, denoted $\kappa$.

We estimate the connectivity of a halo by measuring the number of filaments that intersect the sphere of radius $1.5 \times R_{200c}$ centred at the halo position, similar to Darragh Ford et al. (2019), Gouin et al. (2021). Our results of the mean connectivity as a function of halo mass are presented in Fig. 6. These are number-weighted averages, but we have checked that mass-weighted results are essentially the same.

At all redshifts, we find the expected trend that $\kappa$ increases with mass, in agreement with the studies of, for example,
<table><caption>**Table 1.** Details of the skeletons at different redshifts.</caption><thead><tr><th></th><th>z = 0</th><th>z = 1</th><th>z = 1.99</th><th>z = 3</th><th>z = 4</th></tr></thead><tbody><tr><td>N<sub>CPmax</sub></td><td>41 064</td><td>37 666</td><td>32 289</td><td>28 760</td><td>27 513</td></tr><tr><td>N<sub>fils</sub> (initial)</td><td>169 497</td><td>170 352</td><td>153 808</td><td>144 539</td><td>143 115</td></tr><tr><td>N<sub>fils</sub> (final)</td><td>155 172 (91.5%)</td><td>154 465 (90.7%)</td><td>134 138 (87.2%)</td><td>117 776 (81.5%)</td><td>110 539 (77.3%)</td></tr><tr><td>N<sub>segs</sub> (final)</td><td>790 844</td><td>837 304</td><td>791 972</td><td>742 837</td><td>719 291</td></tr><tr><td>L<sub>min</sub> [pMpc]</td><td>0.20</td><td>0.09</td><td>0.10</td><td>0.09</td><td>0.08</td></tr><tr><td>L<sub>max</sub> [pMpc]</td><td>115.06</td><td>59.02</td><td>36.41</td><td>34.97</td><td>21.00</td></tr><tr><td>L<sub>mean</sub> [pMpc]</td><td>12.38</td><td>6.98</td><td>5.32</td><td>4.35</td><td>3.64</td></tr><tr><td>L<sub>median</sub> [pMpc]</td><td>9.83</td><td>5.61</td><td>4.25</td><td>3.52</td><td>2.94</td></tr><tr><td>L<sub>min</sub> [cMpc]</td><td>0.20</td><td>0.18</td><td>0.30</td><td>0.37</td><td>0.37</td></tr><tr><td>L<sub>max</sub> [cMpc]</td><td>115.06</td><td>118.04</td><td>109.22</td><td>139.88</td><td>104.98</td></tr><tr><td>L<sub>mean</sub> [cMpc]</td><td>12.38</td><td>13.96</td><td>15.95</td><td>17.41</td><td>18.19</td></tr><tr><td>L<sub>median</sub> [cMpc]</td><td>9.83</td><td>11.21</td><td>12.76</td><td>14.07</td><td>14.70</td></tr></tbody></table>

Notes. The length statistics are for the final sample of filaments.

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Evolution of the connectivity $\kappa$ as a function of halo mass. The horizontal lines show the mean connectivity values at a given redshift. From $z = 0$ to 4, these are respectively 4.13, 4.52, 4.77, 5.02, and 5.19. The errorbars are derived from bootstrap resampling.</figcaption></figure>

Aragón-Calvo et al. (2010b), Codis et al. (2018), Sarron et al. (2019), Kraljic et al. (2020b), Lee et al. (2021), Gouin et al. (2021) in both simulations and observations, albeit with different normalisations due to different methods (as presented in the discussion of Malavasi et al. 2020b, 2023). Interestingly, we note the presence of an “elbow” at low halo masses ($M_{200c} < 10^{11} M_{\odot}$) for $z = 3$ and 4, in agreement with the findings of Lee et al. (2021) for the HR5 simulation.

Going beyond the well known trends of halo connectivity with mass, the large statistics of MTNG and the wide redshift range explored in this work allow us to clearly show that halo connectivity consistently decreases from early to late times. Indeed, for all masses, the $z = 0$ connectivity values (purple curve) are always lower than the $z = 1$ ones (in red), which in turn are lower compared to the $z = 2$ ones (in yellow), and so on. These results match the theoretical predictions of Codis et al. (2018, see their Figs. 19 and 23), and are explained by both the nonlinear growth of cosmic structures (causing a decrease of $\kappa$ in the increasingly non-Gaussian fields of the late Universe), and the action of dark energy (DE). Powering the increased expansion of voids at late times, DE can also act on the merging and disconnection of cosmic filaments (see e.g. Cadiou et al. 2020),

thus effectively reducing the connectivity of haloes with respect to the early Universe.

This exploration of halo connectivity as a function of mass and redshift was a natural outcome of the filaments’ robustness checks performed in the previous section. Additional studies should be performed in order to gain more insight on the evolution of connectivity across cosmic time, by, for instance, focusing on specific ranges of halo masses to better interpret the observed “elbow” at early times, or by analysing how these trends change with different DE models. This is however beyond the scope of this paper, whose main goal is to understand the evolution of cosmic web filaments, in particular of their lengths and densities, as presented in the next sections.

# 5. Filament length evolution

In this section we perform a complete characterisation of filament lengths and their evolution with cosmic time. We define the length of a filament, $L_f$, as the sum of the lengths of the individual segments forming the filament, from saddle to peak. We first introduce and analyse the filament length functions $\Psi_f$ in Sect. 5.1. Filament evolutionary tracks are built in Sect. 5.2, and we follow the evolution of filament progenitors in Sect. 5.3 by estimating filament growth rates. The full length statistics of the filaments at the different redshifts are reported in Table 1, in both proper and comoving coordinates.

## 5.1. Filament length functions

We define the filament length function, $\Psi_f$, as the number density of filaments in a bin of length rescaled by the bin width, in a similar spirit as halo mass functions. Figure 7 presents the resulting $\Psi_f$ at different redshifts, both in proper and comoving megaparsecs (hereafter pMpc and cMpc) in the left and right panels, respectively.

In proper coordinates, filament lengths increase with cosmic time. This general trend is evident both from the ranges covered by the different $\Psi_f$ and from the (number-weighted) mean values, marked by the vertical lines and reported in Table 1. For example, at $z = 4$ we measure a mean filament length of 3.64 pMpc while this value is 12.38 cMpc at $z = 0$. This global growth in proper coordinates can be easily understood as a result of the expanding Universe. Tied to the Hubble flow, the cosmic web expands and stretches, resulting in the general increase of
<figure><img src="image_7.png" /><figcaption>Fig. 7. Filament length function Ψf for different redshifts, in proper and comoving coordinates, in the left and right panels, respectively. The thin vertical dotted lines denote the (number-weighted) mean length, reported in Table 1.</figcaption></figure>

cosmic filament lengths. We mention that, in proper coordinates, cosmic filaments also expand along the radial direction, increasing their radial extension with cosmic time (this is shown by the radial density profiles of Appendix C).

In comoving coordinates, once the Hubble flow is factored out, we can probe the specific effects of gravity and DE on the evolution of filament lengths. As revealed by the right panel of Fig. 7, the filament length functions are quite similar and span overall the same ranges, indicating little evolution of the comoving lengths with cosmic time. We note, however, a slight shift of Ψf towards high Lf values at high redshift. This can be simply explained by the more homogeneous density fields at high redshifts, leading to fewer peaks (as seen in Sect. 3.3), and therefore to longer filaments.

These results characterise the lengths of all the cosmic filaments at a given cosmic time, but say little about the evolutionary pathways followed by individual structures. This question is addressed in the next sections, by building the filament histories and tracking the evolution of filament progenitors with time.

## 5.2. Filament evolutionary tracks

We gain further insights on the evolution of filament lengths by associating filaments across different snapshots. These associations aim at identifying the higher-z progenitors of the z = 0 filaments, so one can track the “same” filament back in time and explore its variations in length. In the following, we refer to these associations across cosmic time as the filament history, or evolutionary track. The code used to build them has been adapted from Cadiou et al. (in prep.) and makes use of the scipy and numpy libraries (Virtanen et al. 2020; Harris et al. 2020).

This code works as follows: (i) in comoving coordinates, we compute the pairwise distances of all the peaks (i.e. the CPmax) of the skeleton at redshift $z_1$ to the peaks of the skeleton at a higher redshift $z_2 > z_1$. We limit the computation to a distance of twice the mean nearest-neighbour separation<fnref n="4" />. (ii) We iterate over all pairs (CPmax₁, CPmax₂), ordered by increasing pairwise distance. If CPmax₁ or CPmax₂ have already been visited, we do nothing. Otherwise, we assign CPmax₂ as the progenitor of CPmax₁ and mark both as having been visited. After applying this loop to the skeletons of consecutive redshifts, we obtain lists of peaks tracked as a function of time. (iii) The same two steps are applied to the saddle points S of the skeletons so that we track the saddles (i.e. the second ending points of filaments) as a function of time. (iv) For each filament at z = 0 we locate

its CPmax and saddle point, and we climb up their respective trunks (going to higher redshifts) as long as there exists a filament connecting them. The identified filaments are saved as the progenitors of the filament at z = 0.

This association produces 5061 filaments at z = 0 with a progenitor up to z = 4 (i.e. with complete tracks from redshift zero to four). This is a very low number with respect to the total (3.3%), but can be interpreted by some of the processes that filaments undergo during their evolution, such as mergers, drifting, and twisting due to gravitational interactions with other structures (Cadiou et al. 2020). Studying these dynamical processes would require a robust model of filament formation which, of our knowledge, is currently lacking in the field and goes beyond the scope of this paper. The sample of 5061 filaments with complete histories is thus clearly biased towards structures that had little spatial evolution in the last ~12.25 Gyr (the lookback time at z = 4). We have checked that the length distributions of the associated filaments show similar trends in comoving coordinates as those of the total population at their corresponding redshift. The number of filaments at z = 0 with tracks stopping at z = 1, 2 and 3 is, respectively, 68 268, 31 659, and 14 052.

For further insight, we display in Fig. 8 some visualisations of resulting filament evolutionary tracks. For all the 3D boxes of this figure, we have plotted the z = 0 filament (in purple) and its progenitors at higher redshifts in the same comoving box. The slight drift in comoving position between filament progenitors shows the mitigated effect of gravitational tides in these filaments. In the top row, we see the example of three filaments that have followed an evolutionary path of collapse (or contraction) along the longitudinal direction. Indeed, their lengths (numbers in parenthesis) at z = 0 are smaller than those of their progenitors at higher redshifts. The opposite behaviour of filament growth (or expansion) is illustrated in the bottom row of Fig. 8, where we observe the increase of progenitors' lengths from z = 4 to 0.

## 5.3. Growth rates of cosmic filaments

We now perform a more quantitative investigation of the expansion and contraction along the longitudinal direction of filaments. We compute the filament growth rate Γ defined as:

$$ \Gamma(z_1, z_2) = \frac{L_f(z_1) - L_f(z_2)}{|t(z_1) - t(z_2)|}, \qquad (1) $$

where the redshifts are such that $z_2 > z_1$, the filament length $L_f$ is in units of cMpc, and t is the age of the Universe in units of Gyr at a given redshift. With this definition, collapse or contraction along the longitudinal direction (i.e. $L_f(z_1) < L_f(z_2)$) is characterised by $\Gamma < 0$. Filament growth or expansion with time [$L_f(z_1) > L_f(z_2)$], is measured when $\Gamma > 0$. Since the inputs of

<footnote n="4">⁴ This limit was set by exploring the results with different factors (from 0.5 to 4) of the mean nearest-neighbour separation. Stable results (i.e. small variations in the final number of trees) were found with mean nearest-neighbour separations of two or higher.</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 8. Example of $z = 0$ filaments having complete evolutionary tracks up to $z = 4$. The coloured lines represent the spines of the filaments at different redshifts. Top row: filaments that contracted with time. Bottom row: filaments that expanded with time. The numbers in parenthesis are the length of the filament at a given cosmic time, in comoving megaparsec.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 9. Growth rate distribution of cosmic filaments computed with Eq. (1). The arguments of $\Gamma(z_1, z_2)$ correspond to redshifts such that $z_2 > z_1$.</figcaption></figure>

this quantity are lengths at consecutive redshifts, we use all the available filament progenitors presented in the previous section, even those whose history track is not complete up to $z = 4$. In order to avoid any bias due to extremely wiggly filaments (which represent only a small fraction of the filaments), we take particular care in performing this analysis on relatively straight filaments.

These are the structures whose ratio between length and distance between endpoints (following a straight line) is smaller than 1.2. We have checked that lower values only impact the statistics, without changing the results.

Figure 9 exhibits the distributions of growth rates in intervals of consecutive redshift (i.e. $\Gamma(0, 1)$, $\Gamma(1, 2)$, etc.). Overall, these are rather symmetrically centred around zero, showing that filament contraction is as likely as expansion at all times. Interestingly, the $\Gamma$ ranges significantly narrow down at late times. For example, between redshifts three and four (blue histogram) the maximum magnitude of $\Gamma$ is 50.3 cMpc Gyr$^{-1}$, while it is only 6.6 cMpc Gyr$^{-1}$ for $|\Gamma|$ between redshift zero and one (in red). This suggests a steadier evolution of (comoving) filament lengths in the late Universe compared to early times, when Universe was matter dominated.

### 5.4. Normalised growth rates vs. filament length

Let us now explore the possible dependence of growth rate on filament length. We define the normalised growth rate, $\Gamma_n$, as the growth rate $\Gamma$ divided by the average length between the filament and its progenitor, thus giving the following full expression:

$$ \Gamma_n(z_1, z_2) \equiv \frac{2}{|t(z_1) - t(z_2)|} \times \frac{L_f(z_1) - L_f(z_2)}{L_f(z_1) + L_f(z_2)}. \quad (2) $$

The relation between filament expansion and length needs to be explored using $\Gamma_n$ instead of $\Gamma$, since the latter is intrinsically
<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Normalised growth rate $\Gamma_n$ (Eq. (2)) vs. filament length for four consecutive redshift intervals. For each panel, the x-axis corresponds to the filament lengths of the lowest redshift of the interval. The yellow lines indicate the 68% and 95% contours of the signal (pink colourbar), while the grey contours correspond to a random signal computed by associating filaments and progenitors in a random way, thus breaking the structure of the evolutionary tracks.</figcaption></figure>

correlated with $L_f$ (as longer filaments have higher $\Gamma$ values than short ones, by definition).

Figure 10 presents the distribution of filaments in the $\Gamma_n - L_f$ plane. Each of the four panels of this figure shows a different redshift interval, with the x-axis corresponding to the filament lengths of the lowest redshift. We compare the signal (pink pixels and their corresponding 68% and 95% yellow contours) to the results from random pairings of filaments, presented by the grey contours in the figure. While the latter are quite circular and symmetric along the $\Gamma_n = 0$ axis (as expected from random associations), the former exhibit quite some deviations from this symmetry, demonstrating a positive relation between the normalised growth rate and filament length. Indeed, the tilted and elongated yellow contours of filament progenitors shows that long filaments preferentially expand whereas the short ones contract (with redshift dependent speeds, as shown in the previous section). In particular, we find that, at all times, the longest filaments ($L_f \gtrsim 20$ cMpc) are always expanding in comoving coordinates, and the shortest filaments of lengths $\le 2.5$ cMpc are contracting, as they occupy exclusively the $\Gamma_n < 0$ region.

These results demonstrate that filament expansion and contraction are related with filament length. Further characterising this dependency (by e.g. fitting models to quantify the slope of $\Gamma_n$ with $L_f$) and performing the same analysis in simulations with different DE models could help us understand how to use cosmic filaments as probes of DE. While this approach is beyond the scope of the present work, it would be extremely timely given

the recent and forthcoming galaxy surveys such as DESI, *Euclid*, JPAS, and PFS, which will enable an unprecedented study and characterisation of the large-scale structure of the Universe.

# 6. Evolution of radial density profiles

We finally analyse the evolution of the matter distribution along the radial direction transverse to filaments. This direction is parameterised by the radial coordinate $r$, which is a measure of the perpendicular distance from a mass element to the filament spine. Taking the positions of galaxies as proxies of the distribution of matter, we compute the radial profiles of galaxy density around cosmic filaments and present them in the following.

## 6.1. All filaments

Adopting the same method as in Galárraga-Espinosa et al. (2020, see their Sect. 3.3), we build the average density profiles of galaxies around the filaments presented in Fig. 11. In a nutshell, the method consists of counting the number of galaxies in concentric hollow cylinders of increasing radii, using as axis the spine of the filament, and then dividing this number by the corresponding volume. This procedure is applied to every filament segment defined by DisPerSE (see sketch in Fig. 12), thus producing the same number of profiles as filament segments. For reference, the number of filament segments is reported in Table 1. We compute galaxy over-densities, hereafter $1 + \delta_{gal}$,
<figure><img src="image_11.png" /></figure>

**Fig. 11.** Radial density profile of galaxy over-densities, $1 + \delta_{\text{gal}}$, around cosmic filaments. The black dashed line corresponds to the average profile of all the filament segments, which is compared to the grey profile of Galárraga-Espinosa et al. (2020). The purple and mustard profiles correspond to the average densities in the main filament stem, and in the node outskirts, as described in the main text.

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

**Fig. 13.** Radial density profile of galaxy over-densities $1 + \delta_{\text{gal}}$ around the main stem of cosmic filaments of different redshifts. The radial coordinate $r$ is in units of cMpc.

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

**Fig. 12.** Sketch exemplifying a filament and its partitioning into the different regions considered for the radial density profiles. Most of the profiles of Sect. 6 focus on the main stem of filaments (purple segments).

by dividing the density profiles by the background density of galaxies ($n_{\text{gal}}$, see Sect. 3.1). The average profile of all the filament segments at $z=0$ is displayed by the dashed black curve in Fig. 11. This profile is in excellent agreement with the results of Galárraga-Espinosa et al. (2020) using the TNG300-1 simulation (grey dotted curve).

Since we focus solely on the signal from cosmic filaments, we remove from the analysis the filament segments associated with the nodes of the cosmic web. In all what follows we hence discard the segments inside spheres of radius $R_{200c}$ centred at the position of the FoF haloes, at the ending points of the filaments. These are the segments in the grey zone in the sketch of Fig. 12. The contributions from the node outskirts (defined by the $1-3 \times R_{200}$ regions) are also discarded in order to retrieve the signal of the main stem of filaments, i.e. only of the portions of filaments that are at distances larger than $3 \times R_{200}$ from the centres of the connected nodes. For the sake of clarity, Fig. 12 presents a schematic illustration of the partition of a filament into these different zones.

The average profile of segments along the main stem of filaments is plotted in purple in Fig. 11. This curve shows that densities in the main stem are more concentrated towards filament cores than in the node outskirts (in mustard). Indeed, the latter exhibits larger radial extensions, as expected for these particular environments where filaments and clusters intersect. We note

that these two profiles overlap in filament cores, at densities that are significantly lower than those obtained when the contributions from the nodes are included. In the latter case, the resulting curve (in black) is steeper and more concentrated towards small radii than both the purple and mustard curves, with a shape close to that of NFW haloes. This simply reflects how node contributions are completely dominating over filament densities when taken into account in the computation.

Now targeting the study of the main stem of filaments, we present the resulting average radial profiles at redshifts $z=0, 1, 2, 3$, and $4$ in Fig. 13. We choose to plot these profiles as a function of the comoving distance to the spine ($r$ in units of cMpc) in order to probe trends beyond the expected growth of filaments in the radial direction, caused by the expansion of the Universe. This radial growth can be seen in proper coordinates in Appendix C (which also includes the versions in proper coordinates of all the figures of this section). In Fig. 13, it is striking to see the close similarities between the different profiles. Besides the mild increase of densities at the cores ($r < 1$ cMpc) at late times – in qualitative agreement with Zhu et al. (2021) and probably due to the gravitational collapse along the radial direction – we see a large overlap at larger distances from the spines, at the outskirts of filaments ($r > 1$ cMpc) and further away. This is remarkable considering the large redshift range that is explored in this work, and hints at very little evolution of the comoving radial densities with cosmic time. This is further explored next.

## 6.2. Profiles by bins of comoving length

Let us now split the filaments in bins of comoving length and analyse their corresponding radial density profiles. The selection of filaments is again done according to their comoving lengths, because we aim at probing the evolution of a given population of filaments (e.g. the short, or the long population) across redshift. We compute the average profiles in the four bins of length $L_f < 4, [4, 9], [9, 20],$ and $\ge 20$ cMpc, and the results are shown in the panels of Fig. 14. Interestingly, for all bins, we observe the same trends with redshift as for the total population: the profiles are overall stable with time and show very little variation besides a slight increase of the densities at the cores.

These trends persist when restricting our analysis to the 5061 filaments with complete evolutionary tracks introduced in
<figure><img src="image_14.png" /><figcaption>Fig. 14. Same as Fig. 13 but with filaments separated in bins of comoving length, ordered from left to right by increasing lengths.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 15. Differences D between profiles of consecutive redshifts (Eq. (3)), in black for all the filaments, and in other colours for filaments in a particular bin of comoving length, as labelled.</figcaption></figure>

Sect. 5.2. For the filaments at $z = 0$ in a given bin of length, we compute the average profiles of their progenitors at higher redshifts, so that we track the radial evolution of the “same” filaments across cosmic time. The results of the progenitors are displayed by the coloured squares in Fig. 14. We observe quite a good agreement with the profiles of the general population and with the trends presented above.

For a more quantitative comparison we measure the difference D between the average profiles (p) obtained at $z_1$ and at $z_2$, with respect to their corresponding errors (e), as follows:

$$D(z_1, z_2) = \frac{p(z_1) - p(z_2)}{\sqrt{e(z_1)^2 + e(z_2)^2}} \quad (3)$$

The differences between profiles of consecutive redshifts are reported in Fig. 15, demonstrating that the variations in filament over-densities with time are indeed tiny (<2σ), even at filament cores.

## 6.3. Comment on filament populations

We finally compare the different panels of Fig. 14 with each other. For all the redshifts it is remarkable to see that the shortest filaments have always higher density values than the longest ones. This trend is even more visible in Fig. 16, where we have over-plotted the profiles of the extreme populations at all z, in orange and blue for the shortest ($L_f < 4$ cMpc) and longest filaments ($L_f \ge 20$ cMpc), respectively. The bottom panel of this figure quantifies the difference D (adapted from Eq. (3)) between these two extremes at the same redshift. This evidences a persistent

<figure><img src="image_16.png" /><figcaption>Fig. 16. Evolution of density profiles of extreme bins of comoving length (top panel). The bottom panel presents the significance of the differences (adapted from Eq. (3)) between the longest and shortest average filaments at a given redshift.</figcaption></figure>

~5σ difference between the shortest and the longest populations of all redshifts. In addition, we can observe some differences in the overall shape of the profiles: those of the long population present flatter cores than those of the short one, which rather show steeper slopes at the innermost parts.

The results presented here complement the findings of Galárraga-Espinosa et al. (2020, 2021, 2022) by demonstrating that the diversity of cosmic filaments observed at $z = 0$ is also present at higher redshifts, with the same significance as at $z = 0$. According to those studies, late time short filaments are preferentially located in the denser regions of the cosmic web (e.g. at the vicinity of nodes, see Fig. 1 of Galárraga-Espinosa et al. 2022) whereas long filaments inhabit less dense environments (e.g. delimiting voids). These different locations in the cosmic web are most likely at the origin of the different comoving radial densities between short and long filaments, at all redshifts.

Finally, the density profiles presented in this work may include contributions from objects located at significant radial
distances from a filament but which are most likely not physi-
cally associated with the filament. These might be, for example,
galaxies belonging to an intersecting wall or void, or to other
nodes. It might be interesting to separate the profiles of the fil-
aments themselves from those characterising their environment,
in a similar way as has been done by Ramsøy et al. (2021) for
smaller-scale gaseous filaments. We expect to find significant
differences between different populations of filaments, since the
shortest structures are preferentially located in denser environ-
ments than the longest ones.

7. Summary, discussion, and conclusions

We have used the hydro-dynamical runs of the MTNG simula-
tion at redshifts z = 0, 1, 2, 3, and 4 to build catalogues of cosmic
filaments. We have analysed the spatial evolution of these struc-
tures across cosmic time, as well as the evolution of filament
lengths, growth rates, and radial density profiles.

Following an observational approach, we have used the most
massive galaxies (in stellar mass) at each redshift as tracers of
the density field (Sect. 3.1), and extracted the cosmic skeleton
using the DisPerSE filament finder. The full practical details
related to running DisPerSE on discrete tracers were presented
in Sect. 2.2 with the aim of guiding future users through the steps
of this code. By working at fixed number density of tracers we
have implicitly imposed a scale for the detected filaments. This
allowed us to perform our study on the large-scale cosmic struc-
tures without the contamination of smaller-scale filaments (such
as, e.g. filaments linking satellites with their central galaxies).

An important outcome of this work is the introduction of a
method based on physical priors that allows us to test and quant-
tify the robustness of the cosmic web skeletons that we identify.
This was presented in full detail in Sect. 3.2 with the hope that
future work in this field will address these questions with similar
care. Keeping these considerations in mind can help the commu-
nity move towards studies that are less dependent on the specific
techniques used to detect filaments in the cosmic web.

We have quantified the growth of proper filament lengths and
radial extensions in Figs. 7 (left) and C.1. This growth in proper
coordinates is expected for structures at these large scales as a
result of the cosmic expansion. The most interesting conclusions
have however arisen when working in comoving coordinates,
once the Hubble flow is factored out and similar types of fila-
ment can be compared at different redshifts. These conclusions
are summarised and discussed in the following:

(i) We have found surprisingly stable filament length func-
tions (Fig. 7, right) and over-density profiles (Fig. 13) across
redshift, demonstrating very little evolution of the global pop-
ulation of cosmic filaments over the past ~12.25 Gyr. These
results complement those of Gheller et al. (2016), who have
shown only mild variations of filament scaling relations (e.g.
mass vs. temperature) from redshift 0 to 1. The very small red-
shift dependence of filament over-densities found in this work
shows that filaments are not becoming thicker or thinner with
time in comoving coordinates. This could be directly related
to the continuous accretion of matter from walls or voids, thus
building up mass with time, and compensating for the flow of
matter along the filaments towards nodes. Indeed, filaments can
act as highways for matter transport as suggested by the theory of
large-scale structure formation (Zel'dovich 1970; Peebles 1980;
Zeldovich et al. 1982; Bond et al. 1996). Yet, to our knowledge,
direct and quantitative evidence for this picture is still missing.

(ii) Going beyond the study of global trends, we have associ-
ated filaments across redshifts and built their evolutionary tracks

from *z* = 0 to 4 in Sect. 5.2. This filament association, performed
for the first time in this work, allowed us to follow the evolution
of a relatively small but significant number of structures across
cosmic time. By tracking individual filaments we have shown
that some of them actually significantly change their length with
time (they can contract or elongate, Fig. 10), while the proper-
ties of the global population are preserved (see the point above).
We have quantified the expansion or contraction along the fila-
ment's longitudinal axis by computing growth factors Γ and Γn
(Figs. 9 and 10) finding that, at all times, filament expansion
along the longitudinal direction is related with length. Indeed,
we have shown that filaments longer than a certain length are
preferentially expanding whereas the shortest structures always
collapse. We anticipate that further characterising these trends
by modelling the filament growth in, for example, different DE
scenarios might open doors towards new cosmological probes.

(iii) Conclusions regarding the diversity of cosmic filaments at a given time have also been reached. At a fixed redshift, filament length distributions span more than two orders of magnitude in comoving megaparsecs (Fig. 7, right panel), and the over-density profiles of populations of extreme lengths, i.e. the short and long filaments, always differ by ~5σ (Fig. 16). In the light of previous work (Galárraga-Espinosa et al. 2020, 2021, 2022), we argue that the existence of this diversity of filaments at all redshifts is likely to arise from environmental effects. Different populations inhabit different locations within the cosmic web, where either gravity or DE dominate the local evolution. Within this picture, high redshift cosmic filaments located in high density regions would be more likely to collapse and become part of the short population at later times. Conversely, if the same high-z filaments are found in less dense environments (e.g. next to cosmic voids), they would rather expand, due to the stretching of the cosmic web caused by the DE forces.

Future perspectives for cosmic filament studies beyond those
already mentioned throughout the paper lie in the analysis of
the evolution of gas properties. In particular, investigating the
emergence of the WHIM phase, which dominates the gas frac-
tion (>80%) of cosmic filaments at z = 0, can be insightful for
the study of the late-time missing baryons. The evolution of the
cold gas phase and its relation with filamentary structures should
also be further explored to improve our understanding of galaxy
evolution.

In conclusion, the results presented in this work are an
important step forward in the characterisation of cosmic fil-
aments. With further theoretical work these structures can
potentially become tools for cosmological studies, supplement-
ing their well-established importance in the galaxy evolution
field. Better understanding cosmic filaments is crucial to fully
exploit the potential of this decade's experiments (such as the
already mentioned DESI, *Euclid*, JPAS, *Vera C. Rubin*, or
PFS galaxy surveys) but also the potential of on-going and
future experiments focusing on the detection of gas in different
phases, such as CLAMATO, eROSITA, *Athena*, and the Simons
Observatory.

Acknowledgements. We thank the anonymous referee for their interesting and insightful comments. D.G.E. would like to thank Charlotte Welker, Marcelo Musso, Guinevere Kauffmann, Toni Tuominen, Enrico Garaldi, Fabrizio Arrigoni-Battaia, Joanne Cohn, and the ByoPiC team (https://byopic.eu/) for insightful comments and discussions during the preparation of this work. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958, and a fraction of it was carried out during the KITP Cosmic Web workshop. S.B. is supported by the UK Research and Innovation (UKRI) Future Leaders Fellowship [grant number MR/V023381/1]. C.C. acknowledges support from the Knut and Alice Wallenberg Foundation and the Swedish Research Council (grant 2019-04659). C.G. acknowledges
support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme grant agreement ERC-2015-AdG 695561 (ByoPiC project). Author contributions: DGE led the project: project conceptualisation, methodology, software, validation, formal analysis, investigation, data curation, writing (original draft, review and editing), and visualisation. C.C. helped with project conceptualisation, methodology, software, writing (review and editing). C.G. helped with project conceptualisation and writing (review and editing). S.W. helped with methodology, writing (review and editing), and resources. V.S. and R.P. helped with data curation, writing (review and editing), and resources. B.H., S.B., F.F., R.K., L.H. helped with writing (review and editing), and resources. M.B., A.M.D., C.H.A. helped with resource provision.

## References

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

Aragón-Calvo, M. A., Platen, E., van de Weygaert, R., & Szalay, A. S. 2010a, ApJ, 723, 364

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

Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92

Awad, P., Peletier, R., Canducci, M., et al. 2023, MNRAS, 520, 4517

Barrera, M., Springel, V., White, S. D. M., et al. 2023, MNRAS, 525, 6312

Bartalucci, I., Arnaud, M., Pratt, G. W., et al. 2017, A&A, 608, A88

Baxter, E., Chang, C., Jain, B., et al. 2017, ApJ, 841, 18

Bifvi, V., Dolag, K., Reiprich, T. H., et al. 2022, A&A, 661, A17

Bilicki, M., Peacock, J. A., Jarrett, T. H., et al. 2016, ApJS, 225, 5

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

Bonjean, V., Aghanim, N., Salomé, P., Douspis, M., & Beelen, A. 2018, A&A, 609, A49

Bonjean, V., Aghanim, N., Douspis, M., Malavasi, N., & Tanimura, H. 2020, A&A, 638, A75

Bonnaire, T., Aghanim, N., Decelle, A., & Douspis, M. 2020, A&A, 637, A18

Bose, S., Hadzhiyska, B., Barrera, M., et al. 2023, MNRAS, 524, 2579

Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151

Bryant, J. J., Owers, M. S., Robotham, A. S. G., et al. 2015, MNRAS, 447, 2857

Buncher, B., & Carrasco Kind, M. 2020, MNRAS, 497, 5041

Busch, P., & White, S. D. M. 2020, MNRAS, 493, 5693

Cadiou, C., Pichon, C., Codis, S., et al. 2020, MNRAS, 496, 4787

Carrón Duque, J., Migliaccio, M., Marinucci, D., & Vittorio, N. 2022, A&A, 659, A166

Cautun, M., van de Weygaert, R., & Jones, B. J. T. 2013, MNRAS, 429, 1286

Cautun, M., van de Weygaert, R., Jones, B. J. T., & Frenk, C. S. 2014, MNRAS, 441, 2923

Cen, R., & Ostriker, J. P. 2006, ApJ, 650, 560

Chen, Y.-C., Ho, S., Mandelbaum, R., et al. 2017, MNRAS, 466, 1880

Codis, S., Pogosyan, D., & Pichon, C. 2018, MNRAS, 479, 973

Colberg, J. M., Krughoff, K. S., & Connolly, A. J. 2005, MNRAS, 359, 272

Colless, M., Peterson, B. A., Jackson, C., et al. 2003, arXiv e-prints [arXiv:astro-ph/0306581]

Connor, T., Kelson, D. D., Mulchaey, J., et al. 2018, ApJ, 867, 25

Contreras, S., Angulo, R. E., Springel, V., et al. 2023, MNRAS, 524, 2489

Darragh Ford, E., Laigle, C., Gozaliasl, G., et al. 2019, MNRAS, 489, 5695 de Graaff, A., Cai, Y.-C., Heymans, C., & Peacock, J. A. 2019, A&A, 624, A48

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

Delgado, A. M., Hadzhiyska, B., Bose, S. D. M., et al. 2023, MNRAS, 523, 5899

Driver, S. P., Hill, D. T., Kelvin, L. S., et al. 2011, MNRAS, 413, 971

Erciyes, C., Basu, K., Kim, S., & Rey, S.-C. 2023, A&A, 675, A63

Ferlito, F., Springel, V., Davies, C. T., et al. 2022, MNRAS, 524, 5591

Galárraga-Espinosa, D., Aghanim, N., Langer, M., Gouin, C., & Malavasi, N. 2020, A&A, 641, A173

Galárraga-Espinosa, D., Aghanim, N., Langer, M., & Tanimura, H. 2021, A&A, 649, A117

Galárraga-Espinosa, D., Langer, M., & Aghanim, N. 2022, A&A, 661, A115

Galárraga-Espinosa, D., Garaldi, E., & Kauffmann, G. 2023, A&A, 671, A160

Ganeshaiah Veena, P., Cautun, M., van de Weygaert, R., et al. 2018, MNRAS, 481, 414

Ganeshaiah Veena, P., Cautun, M., Tempel, E., van de Weygaert, R., & Frenk, C. S. 2019, MNRAS, 487, 1607

Ganeshaiah Veena, P., Cautun, M., van de Weygaert, R., Tempel, E., & Frenk, C. S. 2021, MNRAS, 503, 2280

Garaldi, E., Kannan, R., Smith, A., et al. 2023, MNRAS, submitted [arXiv:2309.06475]

Gheller, C., & Vazza, F. 2019, MNRAS, 486, 981

Gheller, C., Vazza, F., Favre, J., & Brüggen, M. 2015, MNRAS, 453, 1164

Gheller, C., Vazza, F., Brüggen, M., et al. 2016, MNRAS, 462, 448

Ghirardini, V., Ettori, S., Eckert, D., & Molendi, S. 2019, A&A, 627, A19

Gouin, C., Aghanim, N., Bonjean, V., & Douspis, M. 2020, A&A, 635, A195

Gouin, C., Bonnaire, T., & Aghanim, N. 2021, A&A, 651, A56

Gouin, C., Gallo, S., & Aghanim, N. 2022, A&A, 664, A198

Gregory, S. A., & Thompson, L. A. 1978, ApJ, 222, 784

Guzzo, L., Scodeggio, M., Garilli, B., et al. 2014, A&A, 566, A108

Hadzhiyska, B., Eisenstein, D., Hernquist, L., et al. 2023a, MNRAS, 524, 2507

Hadzhiyska, B., Hernquist, L., Eisenstein, D., et al. 2023b, MNRAS, 524, 2524

Hahn, O., Carollo, C. M., Porciani, C., & Dekel, A. 2007, MNRAS, 381, 41

Hahn, C., Wilson, M. J., Ruiz-Macias, O., et al. 2023, AJ, 165, 253

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Hernández-Aguayo, C., Springel, V., Pakmor, R., et al. 2023, MNRAS, 524, 2556

Ho, A., Gronke, M., Falck, B., & Mota, D. F. 2018, A&A, 619, A122

Holt, P., Tuominen, T., Nevalainen, J., et al. 2022, MNRAS, 513, 3387

Inoue, S., Si, X., Okamoto, T., & Nishigaki, M. 2022, MNRAS, 515, 4065

Jöveer, M., Einasto, J., & Tago, E. 1978, MNRAS, 185, 357

Jhee, H., Song, H., Smith, R., et al. 2022, ApJ, 940, 2

Jones, D. H., Read, M. A., Saunders, W., et al. 2009, MNRAS, 399, 683

Kannan, R., Springel, V., Hernquist, L., et al. 2023, MNRAS, 524, 2594

Klar, J. S., & Mücket, J. P. 2012, MNRAS, 423, 304

Kotecha, S., Welker, C., Zhou, Z., et al. 2022, MNRAS, 512, 926

Kraljic, K., Arnouts, S., Pichon, C., et al. 2018, MNRAS, 474, 547

Kraljic, K., Pichon, C., Dubois, Y., et al. 2019, MNRAS, 483, 3227

Kraljic, K., Davé, R., & Pichon, C., 2020a, MNRAS, 493, 362

Kraljic, K., Pichon, C., Codis, S., et al. 2020b, MNRAS, 491, 4294

Kuchner, U., Haggar, R., Aragón-Salamanca, A., et al. 2022, MNRAS, 510, 581

Kuutma, T., Tamm, A., & Tempel, E. 2017, A&A, 600, L6

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Laigle, C., Pichon, C., Arnouts, S., et al. 2018, MNRAS, 474, 5437

Lee, J., Shin, J., Snaith, O. N., et al. 2021, ApJ, 908, 11

Le Fèvre, O., Vettolani, G., Garilli, B., et al. 2005, A&A, 439, 845

Liao, S., & Gao, L. 2019, MNRAS, 485, 464

Libeskind, N. I., van de Weygaert, R., Cautun, M., et al. 2018, MNRAS, 473, 1195

Lu, Y. S., Mandelker, N., Oh, S. P., et al. 2024, MNRAS, 527, 11256

Malavasi, N., Arnouts, S., Vibert, D., et al. 2017, MNRAS, 476, 3817

Malavasi, N., Aghanim, N., Douspis, M., Tanimura, H., & Bonjean, V. 2020a, A&A, 642, A19

Malavasi, N., Aghanim, N., Tanimura, H., Bonjean, V., & Douspis, M. 2020b, A&A, 634, A30

Malavasi, N., Langer, M., Aghanim, N., Galárraga-Espinosa, D., & Gouin, C. 2022, A&A, 658, A113

Malavasi, N., Sorce, J. G., Dolag, K., & Aghanim, N. 2023, A&A, 675, A76

Martizzi, D., Vogelsberger, M., Artale, M. C., et al. 2019, MNRAS, 486, 3766

Musso, M., Cadiou, C., Pichon, C., et al. 2018, MNRAS, 476, 4877

Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Nelson, D., Springel, V., Pillepich, A., et al. 2019, Comput. Astrophys. Cosmol., 6, 2

Nevalainen, J., Tempel, E., Liivamägi, L. J., et al. 2015, A&A, 583, A142

Neyrinck, M. C. 2012, MNRAS, 427, 494

Pakmor, R., Springel, V., Coles, J. P., et al. 2023, MNRAS, 524, 2539

Pandey, B., & Bharadwaj, S. 2006, MNRAS, 372, 827

Peebles, P. J. E. 1980, The Large-scale Structure of the Universe (Princeton: Princeton University Press)

Pereyra, L. A., Sgró, M. A., Merchán, M. E., Stasyszyn, F. A., & Paz, D. J. 2020, MNRAS, 499, 4876

Pfeifer, S., Libeskind, N. I., Hoffman, Y., et al. 2022, MNRAS, 514, 470

Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077

Pintos-Castro, I., Yee, H. K. C., Muzzin, A., Old, L., & Wilson, G. 2019, ApJ, 876, 40

Planck Collaboration XIII. 2016, A&A, 594, A13

Ramsøy, M., Slyz, A., Devriendt, J., Laigle, C., & Dubois, Y. 2021, MNRAS, 502, 351

Rost, A., Stasyszyn, F., Pereyra, L., & Martínez, H. J. 2020, MNRAS, 493, 1936

Rost, A., Kuchner, U., Welker, C., et al. 2021, MNRAS, 502, 714

Santiago-Bautista, I., Caretta, C. A., Bravo-Alfaro, H., Pointecouteau, E., & Andernach, H. 2020, A&A, 637, A31

Sarron, F., Adami, C., Durret, F., & Laigle, C. 2019, A&A, 632, A49

Schaap, W. E., & van de Weygaert, R. 2000, A&A, 363, L29

Schaye, J., Kugel, R., Schaller, M., et al. 2023, MNRAS, 526, 4978

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Shull, J. M., Smith, B. D., & Danforth, C. W. 2012, ApJ, 759, 23

Singh, A., Mahajan, S., & Bagla, J. S. 2020, MNRAS, 497, 2265

Song, H., Laigle, C., Hwang, H. S., et al. 2021, MNRAS, 501, 4635

Sousbie, T. 2011, MNRAS, 414, 350

Sousbie, T., Pichon, C., & Kawahara, H. 2011, MNRAS, 414, 384
Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629

Springel, V., Pakmor, R., Zier, O., & Reinecke, M. 2021, MNRAS, 506, 2871

Stoica, R. S., Martínez, V. J., & Saar, E. 2007, J. R. Stat. Soc.: Ser. C (Appl. Stat.), 56, 459

Tanimura, H., Hinshaw, G., McCarthy, I. G., et al. 2019, MNRAS, 483, 223

Tanimura, H., Aghanim, N., Bonjean, V., Malavasi, N., & Douspis, M. 2020a, A&A, 637, A41

Tanimura, H., Aghanim, N., Kolodzig, A., Douspis, M., & Malavasi, N. 2020b, A&A, 643, L2

Tanimura, H., Aghanim, N., Douspis, M., & Malavasi, N. 2022, A&A, 667, A161

Taylor, E. N., Hopkins, A. M., Baldry, I. K., et al. 2011, MNRAS, 418, 1587

Tempel, E., Stoica, R. S., Martínez, V. J., et al. 2014, MNRAS, 438, 3465

Tempel, E., Stoica, R. S., Kipper, R., & Saar, E. 2016, Astron. Comput., 16, 17

Tuominen, T., Nevalainen, J., Tempel, E., et al. 2021, A&A, 646, A156

van de Weygaert, R., & Schaap, W. 2009, in The Cosmic Web: Geometric Analysis, eds. V. J. Martínez, E. Saar, E. Martínez-González, & M. J. Pons-Bordería (Berlin: Springer), 665, 291

Vernstrom, T., Heald, G., Vazza, F., et al. 2021, MNRAS, 505, 4178

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Meth., 17, 261

Wang, P., Libeskind, N. I., Tempel, E., Kang, X., & Guo, Q. 2021, Nat. Astron., 5, 839

Welker, C., Bland-Hawthorn, J., Van de Sande, J., et al. 2020, MNRAS, 491, 2864

White, S. D. M., Frenk, C. S., Davis, M., & Efstathiou, G. 1987, ApJ, 313, 505

York, D. G., Adelman, J., Anderson, J. E., Jr., et al. 2000, AJ, 120, 1579

Zakharova, D., Vulcani, B., De Lucia, G., et al. 2023, MNRAS, 525, 4079

Zel'dovich, Y. B. 1970, A&A, 5, 84

Zeldovich, I. B., Einasto, J., & Shandarin, S. F. 1982, Nature, 300, 407

Zhu, W., Zhang, F., & Feng, L.-L. 2021, ApJ, 920, 2
## Appendix A: Spatial distribution of first progenitors of the subhaloes at $z=0$

<figure><img src="image_17.png" /><figcaption>**Fig. A.1.** Spatial distribution at different redshifts of the first progenitors of the subhaloes of stellar masses $log_{10}(M_*/M_\odot) \ge 9$ at $z=0$. The different panels show the same slice of thickness 15 cMpc as the one in Fig. 3.</figcaption></figure>

Figure A.1 shows, for the different redshifts explored in this work, the spatial distribution of the first progenitors of the $z=0$ galaxies of stellar masses $log_{10}(M_*/M_\odot) \ge 9$. It is clearly apparent that these galaxies are not optimal tracers of the cosmic skeleton. Indeed, as we move to higher redshift, the filamentary pattern (finely delineated at $z=0$) becomes noisier, hard to distinguish at $z=3$, and absent by $z=4$. This figure can be compared with Fig. 3 of the main text, which shows, in the same slice, that the more massive galaxies at any redshift are better tracers of the cosmic skeleton. More quantitatively, we find that the first progenitors of the selected $z=0$ galaxies have stellar masses as low as $log_{10}(M_*/M_\odot) = 6.47$ at $z=1$. By $z=4$, only

76% of them have a non zero stellar component, the rest being only low mass DM clumps.

## Appendix B: Matching the CPmax with massive haloes at $z=0$

The full results of the matching of CPmax with FoF haloes are presented in Table B.1, together with all the details regarding the number of filaments, CPmax, and length statistics of each of the eight skeletons extracted at $z=0$ with different DisPerSE parameterisation.

<table><caption>**Table B.1.** Details of DisPerSE skeletons extracted with different parameterisation at $z=0$.</caption><thead><tr><th>ρ<sub>DTFE</sub> smoothing</th><th>S</th><th>S</th><th>S</th><th>S</th><th>noS</th><th>noS</th><th>noS</th><th>noS</th></tr><tr><th>Persistence</th><th>1σ</th><th>2σ</th><th>3σ</th><th>4σ</th><th>1σ</th><th>2σ</th><th>3σ</th><th>4σ</th></tr></thead><tbody><tr><td>Number filaments</td><td>318 337</td><td>169 497</td><td>114 269</td><td>75 065</td><td>697 234</td><td>395 033</td><td>255 200</td><td>172 280</td></tr><tr><td>Number CPmax</td><td>61 035</td><td>41 064</td><td>29 600</td><td>18 157</td><td>173 370</td><td>114 399</td><td>76 777</td><td>51 040</td></tr><tr><td>Min length [Mpc]</td><td>0.0001</td><td>0.20</td><td>0.23</td><td>0.23</td><td>0.0001</td><td>0.001</td><td>0.003</td><td>0.03</td></tr><tr><td>Max length [Mpc]</td><td>82.16</td><td>115.61</td><td>127.41</td><td>146.77</td><td>72.0</td><td>71.5</td><td>105.5</td><td>147.4</td></tr><tr><td>Mean length [Mpc]</td><td>10.5</td><td>12.3</td><td>14.7</td><td>20.6</td><td>6.4</td><td>7.3</td><td>8.8</td><td>11.0</td></tr><tr><td>Median length [Mpc]</td><td>8.9</td><td>9.8</td><td>11.3</td><td>15.9</td><td>5.1</td><td>5.5</td><td>6.3</td><td>7.9</td></tr><tr><td>Galactic haloes hosting CPmax</td><td>2.0%</td><td>0.9%</td><td>0.51%</td><td>0.25%</td><td>10.1%</td><td>5.9%</td><td>3.1%</td><td>1.5%</td></tr><tr><td>Groups hosting CPmax</td><td>48.3%</td><td>36.4%</td><td>27.3%</td><td>16.9%</td><td>74.2%</td><td>65.3%</td><td>55.4%</td><td>43.4%</td></tr><tr><td>Clusters hosting CPmax</td><td>97.0%</td><td>91.3%</td><td>80.9%</td><td>59.1%</td><td>99.7%</td><td>98.7%</td><td>96.7%</td><td>91.3%</td></tr><tr><td>CPmax in galactic haloes</td><td>21.6%</td><td>14.6%</td><td>11.2%</td><td>8.7%</td><td>37.5%</td><td>33.2%</td><td>26.3%</td><td>18.7%</td></tr><tr><td>CPmax in Groups</td><td>55.7%</td><td>62.3%</td><td>64.8%</td><td>65.6%</td><td>30.4%</td><td>40.3%</td><td>50.9%</td><td>59.9%</td></tr><tr><td>CPmax in Clusters</td><td>9.9%</td><td>13.7%</td><td>16.7%</td><td>20.0%</td><td>4.7%</td><td>5.7%</td><td>7.8%</td><td>11.0%</td></tr><tr><td>Total fraction of matched CPmax</td><td>87.1%</td><td>90.6%</td><td>92.8%</td><td>94.4%</td><td>72.7%</td><td>79.1%</td><td>84.9%</td><td>89.6%</td></tr><tr><td>Comments</td><td>impure</td><td>ok</td><td>incomplete</td><td>incomplete</td><td>impure</td><td>impure</td><td>impure</td><td>pb lengths</td></tr></tbody></table>

**Notes.** Results after cleaning the skeletons.
**Appendix C: Density profiles in proper radial coordinates**

In proper coordinates, we observe the expected trends that the density values and radial extensions of filaments increase with time. This is shown by Figs. C.1 and C.2 of this appendix.

<figure><img src="image_18.png" /><figcaption>Fig. C.1. Same as Fig. 13 but with *r* in units of proper megaparsec.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. C.2. Same as Fig. 14 but with *r* in units of proper megaparsec.</figcaption></figure>