We utilize a large amount of simulations that was simulated using the COmoving Lagrangian Acceleration (COLA) code (Tassev et al., 2013). COLA calculates the movement of dark matter particles in a comoving frame, following the trajectories forecasted by Lagrangian Perturbation Theory (LPT), enabling it to accurately handle small-scale structures while preserving accuracy on large scales. Despite being hundreds of times faster than N-body simulations, it maintains high precision from very large to highly non-linear scales. To explore various aspects of filament length statistics and evaluate the impact of numerical noise on analysis results and test robustness, we utilize COLA to generate 700 pure dark matter cosmological simulations with different parameters. The parameters of these simulation are listed in Table.1 and they are categorized as follows,

1. 1)

   In the fiducial dataset, we conduct 50 standard cosmological simulations with the cosmological parameters of Planck15, with $\Omega_{m}=0.308$, $\Omega_{b}=0.05$, $\sigma_{8}=0.83$, $n_{s}=0.968$, and $H_{0}=67.8,\text{km s}^{-1}\text{Mpc}^{-1}$. Each simulation box has a size of $512^{3}(h^{-1}\text{Mpc})^{3}$ and contains $256^{3}$ particles.

2. 2)

   P1 and P2 datasets have the same cosmological parameters as the fiducial one, but with different number of particles: $128^{3}$ and $512^{3}$ respectively. These two datasets are used to investigate the impact of varying particle numbers on filament searches and statistics. Each dataset contains five realizations with different random seeds for producing initial conditions.

3. 3)

   The C1 dataset comprises 80 cosmological simulations, with varying $\Lambda$CDM parameters: $\Omega_{m}$, $\sigma_{8}$, and $n_{s}$. Our goal is to analyze changes in filament length statistics across various cosmological parameters and determine their potential for constraining these parameters. To achieve this, we conduct simulations where the value of a single parameter is varied with respect to the fiducial simulations. Each simulation utilizes the same box size and particle number as the fiducial dataset. To ensure robustness, we also generate five realization for each simulation using distinct random seeds. The cosmological parameters for the C1 dataset ranged from [0.2075, 0.4075] for $\Omega_{m}$, [0.60, 1.00] for $\sigma_{8}$, and [0.60, 1.20] for $n_{s}$.

4. 4)

   The SC1 dataset closely resembles C1, except for a box size of $300^{3}(h^{-1}{\rm Mpc})^{3}$ and $512^{3}$ particles, encompassing 80 cosmological simulations. This dataset serves to assess the robustness of our findings across varying simulation box sizes and precision levels.

5. 5)

   The I1 dataset comprises 240 cosmological simulations which have modified initial power spectrum. The goal is to study how filament length statistics are affected by initial power spectrum. We vary the initial fluctuations at eight different scales, with each scale containing six different fluctuation amplitudes, with a total of 48 simulations. For each modified power spectrum, we generate five realizations with different random seeds, using the same box size and particle number as the fiducial simulation.

6. 6)

   The SI1 dataset closely resembles I1, except for a box size of $300^{3}(h^{-1}{\rm Mpc})^{3}$ and $512^{3}$ particles, consisting of 240 cosmological simulations. This dataset serves to assess the robustness of our findings across varying simulation box sizes and precision levels.

By running COLA, we gain access to a diverse range of cosmological simulations, enabling us to extract valuable physical insights and comprehend the intricate, nonlinear nature of cosmic filaments. This process aids in the development of more robust statistical methods for observational constraints on cosmological models.

After acquiring a cosmological simulation, we utilize the CIC (Cloud-in-Cell) scheme (Bertschinger, 1998) to allocate dark matter particles onto a mesh grid, thereby generating density fields for specific snapshots similar to Galárraga-Espinosa et al. (2022). In this work we only focus on the properties of filaments at five specific redshift, from z=0, 0.5,1,2,3. We then apply the filament finder, DisPerSE, to the density field to identify filaments and study their properties and dependence on cosmological parameters.

It is worth noting that, since we can only observe galaxies, filament finders should ideally be applied to the distribution of galaxies, as done by (Galárraga-Espinosa et al., 2020, 2021, 2022, 2023, 2024). One major concern in these studies arises from the need to construct the underlying dark matter density using galaxy distribution, which is inevitably affected by the uncertainty of galaxy bias. In our current studies, we aim to investigate the cosmological dependence of filament distribution by using the dark matter density field directly from simulations. This approach avoids any bias from the sparse tracers of the cosmic density field, such as galaxies or haloes. However, to extract constraints on cosmology, one will ultimately have to use observed galaxies to identify cosmic filaments, which will be a focus of future studies.

We employ the filament finder DisPerSE to extract filamentary structure from the DM density field. DisPerSE (Sousbie, 2011) is a publicly available code that identifies the critical points of the density field using the Morse method, where $S=4,3,2,1$ corresponds to bifurcation, node, saddle, and void, respectively. The saddle connects two nodes to form a filament (Kraljic et al., 2019), and bifurcation is where multiple filaments intersect, which is considered unresolved clusters and the location where superclusters exist. Building upon these considerations and inspired by the cluster-filament-cluster configuration proposed by Bond et al. (1996), we define a complete filament as a high-density line connecting two nodes or bifurcations, or linking a cluster and a bifurcation which is the same as the definition proposed by Wang et al. (2024).

DisPerSE segments the identified filaments into multiple parts, which are subsequently reconnected to form complete filaments based on the definition above. Fig. 1 illustrates the critical points and segments of filaments obtained using DisPerSE, similar to the figure by Galárraga-Espinosa et al. (2020), although with a distinct definition of an individual filament. We use black and brown colors to distinguish different filaments and numbers to tag each filament in the figure, where there are a total of 11 complete filaments.

Ultimately, we define the length of a filament, denoted as $L$, as the sum of the lengths of the individual segments comprising the complete filament. This summation encompasses segments between nodes to nodes, nodes to bifurcations, or bifurcations to bifurcations. By computing the lengths of all complete filaments, we derive the filament length function. Here, we define the number density of filaments in each length bin as $\Phi_{L}$.

Here, $N(L)$ represents the number of filaments in each length bin, and $V_{boxsize}$ is the volume of the simulation box. Note that in this paper all quantitative relationships are derived with this number density $\Phi_{L}$, but one should be aware that cosmological fluctuation is significant when $N(L)$ is smaller. In addition, we define the relative difference of $\Phi_{L}$ as,

to show the comparison of filament length distribution between different parameters. For example, $T(L)$ denotes the deviation of the filament length distribution respect to the reference value.

The filament identification by DisPerSE is influenced by several factors: the number of grid points, denoted as $N_{grid}$, the smoothing scale of the field, denoted as $N_{smooth}$, and topological persistence, denoted as $N_{\sigma}$. Given our emphasis on the cosmological dependence of filament length distribution, it is vital to ensure that our field structure is accurate and robust, enabling us to draw essential physical conclusions while minimizing the introduction of numerical noise.

Coarse fields and smaller $N_{grid}$ leads to the elimination of small-scale information and minor structural details. we have employed a boxsize of 512 $h^{-1}\text{Mpc}$ with 512-grid point setup, ensuring the accuracy on small scales and providing ample resolution for the cosmological scale under investigation.

On the other hand, the smoothing scale $N_{smooth}$ involves eliminating the shot noise introduced by the lattice process from point cloud to field. Following previous studies on the optimization of smoothing length (Hahn et al., 2007; Cautun et al., 2014), we use a Gaussian smoothing with $N_{smooth}=2h^{-1}{\rm Mpc}$ to get the density field. This smoothing effectively eliminates shot noise in the density field below a mass scale of $10^{13}M_{\odot}h^{-1}$ and filaments with lengths below $2h^{-1}{\rm Mpc}$.

$N_{\sigma}$ characterizes topological persistence, aiding in the identification of robust components of the filamentary network in comparison to a discrete random Poisson distribution. To determine the optimal value for the topological persistence parameter $N_{\sigma}$, we employ the $CP_{max}-Halo$ calibration method to refine our pipeline settings. This method has been widely utilized in previous studies (e.g., (Galárraga-Espinosa et al., 2024, 2023; Boldrini & Laigle, 2024)). The primary concept of this method involves comparing the positions of nodes ($CP_{\text{max}}$) to those of the most massive halos identified by the Friend-of-Friends (FoF) algorithm (Davis et al., 1985), which is widely used to find dark matter halo by linking all particles into one halo when the distance between any two particles is less than f times the average distance between the particles in the simulation box, where f is usually taken as 0.2. The results are depicted in Fig. 13. Following the approach of Boldrini & Laigle (2024), we present the ratio of the fraction of $CP_{\text{max}}$ located within the virial radius ($r_{\text{vir}}$) of any massive halo to the fraction of massive halos hosting a $CP_{\text{max}}$ within their $r_{\text{vir}}$ spheres. This comparison aims to identify the topological persistence value that brings the ratio closest to unity, thus minimizing both the number of halos lacking nodes association and the number of nodes lacking halo association at $z=0$. Consequently, we select $N_{\sigma}=3$ which corresponds to the ratio we defined as unity. Notably, we designate the most massive halos as those with masses exceeding $10^{14}M_{\odot}h^{-1}$ to ensure compatibility with the smoothing scale encompassing $10^{13}M_{\odot}h^{-1}$. Detail analysis of the choosing $N_{\sigma}$ can be found in Appendix.A.

This complete pipeline, starting from the initial power spectrum and given different cosmological parameters, allows us to obtain the statistics of filament lengths. A short summarize of the pipeline involves the following steps,

1. 1)

   Input the initial power spectrum and cosmological parameters, then employ COLA to generate snapshots at various redshifts.

2. 2)

   Utilize the Cloud-in-Cell (CIC) method on the snapshot acquired at a specific redshift to derive a smoothed density field.

3. 3)

   Employ the DisPerSE method to extract the diverse components of filaments.

4. 4)

   Construct a complete list of filaments according to the cluster-filament-cluster definition.

5. 5)

   Obtain various statistical measures for each filament.

   Therefore, we can utilize this pipeline to examine numerous filament statistics, and the pipeline indicates that the statistics are solely dictated by the initial power spectrum and cosmological parameters.

In N-body simulation, one have to sample the power spectrum of the underlined cosmological density field, thus the final mass distribution is inevitable to be affected by the parameters used in N-body simulation, such as the random seed, particle number and box size of the simulation. Here we analyze their impacts one by one. In our fiducial simulation, we have 5 realization of a give cosmological parameters, namely with $\Omega_{m}=0.308$, $\Omega_{b}=0.05$, $\sigma_{8}=0.83$, $n_{s}=0.968$, and $H_{0}=67.8{\rm kms}^{-1}{\rm Mpc}^{-1}$ In the left panel of Fig. 2, we show the filament length distributions of these realization. As can be seen from the left panel, the cosmological fluctuation of the filament length distribution is very small, with only a small scatter in filaments longer than $85h^{-1}{\rm Mpc}$.

To assess the influence of box size on filament length statistics, we examined the normalized filament length distribution using simulations with box sizes of $300(h^{-1}{\rm Mpc})^{3}$ from the SC1 dataset and $512(h^{-1}{\rm Mpc})^{3}$ within fiducial dataset, as depicted in the middle panel of Fig. 2, with 5 random seeds for each box size. Our analysis reveals that the variance resulting from box size is negligible and does not affect our results significantly. In principle, a larger box is preferable for analysis as it provides more filament samples and longer filaments. However, due to limited computing resources, we opt for a $512(h^{-1}{\rm Mpc})^{3}$ box, which is sufficient for our current study.

In the right panel of Fig. 2, we compare the effects of varying particle numbers $N_{p}$ in numerical simulations on filament length statistics. Within a $512(h^{-1}{\rm Mpc})^{3}$ box, the distribution of filament lengths exceeding $8h^{-1}{\rm Mpc}$ is similar for particle numbers of $256^{3}$ and $512^{3}$, while simulations with $128^{3}$ particles differ significantly. This is because a low number density in the simulation results in a different distribution of the CIC density field. Therefore, to ensure both the accuracy of the density field and the efficiency of the simulation, we use $N_{p}=256^{3}$.

The above results indicate that cosmological fluctuations including different random seed, different boxsize and different particle numbers do not significantly influence the filament length distribution. This outcome reinforces the viability of utilizing filament length distribution as a method for studying cosmological constraints.

As mentioned in Section 1, the large-scale mass distribution is significantly affected by the cosmological parameters, as does the cosmic filament. However, few studies have focused on the effects of cosmology on filament length distribution. Using the filament pipeline discussed in Sec.2.2, we have investigated the relationship between filament length distribution and three different cosmological parameters in Fig. 3 where we show the results by varying $\sigma_{8}$, $\Omega_{m}$ and initial power spectrum index $n_{s}$, respectively. Overall, it is found that filament length distribution exhibits varying degrees of sensitivity to these three cosmological parameters, and the strength of this signal is significantly beyond the cosmic variance.

In the left panel of Fig. 3, we observe that as $\sigma_{8}$ increases, the number of filaments, particularly shorter ones, increases, while the number of longer filaments decreases. The decrease of $\sigma_{8}$ causes the universe to become more homogeneous, while an extremely low value of $\sigma_{8}$ would leave the filament structure nearly empty. As $\sigma_{8}$ increases and the clustering of each scale strengthens, more voids and peaks emerge, leading to an increase in the number of short filaments. However, as $\sigma_{8}$ reaches higher values, strong clustering results in increased merging, causing long filaments to break and additional peaks to form within them. This phenomenon may lead to the fragmentation of longer filaments into smaller ones, thereby contributing to an elevated number of short filaments. To demonstrate this, we generate two-dimensional slice plots of the density field and employ the color rendering method from $pmwd$ (Li et al., 2022) to effectively depict the cosmic structures of $\sigma_{8}=0.59$, $\sigma_{8}=0.82$, and $\sigma_{8}=1.66$ at $z=0$ in a $256^{3}(h^{-1}{\rm Mpc})^{3}$ box. This is illustrated in Fig. 4, where the absence of filamentary structure at low $\sigma_{8}$ is notable on the left panel. As $\sigma_{8}$ increases, more peaks form and the filament structure becomes apparent. However, with further increases in $\sigma_{8}$ to much higher values, more structures merge, leading to the break of long filaments and the emergence of additional peaks that break down the filaments. Consequently, the number of long filaments decreases while the number of short filaments increases.

Seen from the middle panel of Fig. 3, it is evident that as $\Omega_{m}$ increases, short filaments increase while long filaments diminish. This phenomenon is explicable: the greater proportion of matter fosters nearby merging, each culminating in a small peak, thereby impeding the formation of extended structures.

While an increase in $n_{s}$ exhibits similar behavior to $\Omega_{m}$, it results in a stronger change in filament length distribution, as shown on the right panel of Fig. 3. The features are very similar to those of $\Omega_{m}$, with extensive structures present at small $n_{s}$, and as $n_{s}$ increases, smaller structures displace the extensive ones, accompanied by the emergence of numerous peaks. The similar behavior of $n_{s}$ and $\Omega_{m}$, and the significant change in filament length distribution for $n_{s}$, originate from the initial power spectrum, where a larger $n_{s}$ or $\Omega_{m}$ corresponds to a smaller amplitude in small k modes.

The sensitivity of filament length distribution to different cosmological parameters persists when we change the number of grid $N_{grid}$ or the smoothing scale $N_{smooth}$ as well as the topological persistence $N_{\sigma}$ in Appendix.B.

Notably, the abundance of filaments exhibits a crossover at approximately $60h^{-1}{\rm Mpc}$ with varying cosmological parameters ($\sigma_{8}$, $\Omega_{m}$, $n_{s}$). This specific value is slightly dependent on different parameter selections, such as $N_{grid}$ and $N_{smooth}$. Therefore, this crossover point signifies a characteristic scale for filaments that warrants further exploration.

The filament length distribution is shown in Fig. 5, and it can be well fitted with the following equation Eq. 3,

Here the parameters $b$ and $n_{0}$ are dependent on the cosmological parameters, and we show the best fitted $b$ and $n_{0}$ to different cosmological parameters in Fig. 6. For each cosmology, we use five realization with distinct initial matter distributions to calculate the cosmic variance for $n_{0}$ and $b$ values and plot 1 $\sigma$ dispersion in the figure. As evident from the left and right panels, the relation is tighter for $\sigma_{8}$ and $n_{s}$, while the relation for $\Omega_{m}$ on the middle panel shows slight stronger sensitivity to cosmic variance.

It is well known from the CDM theory that the number density of massive dark matter haloes is sensitive to cosmological parameters, thus being a powerful probe of cosmology although it is difficult to obtain halo mass from observations. As an interesting theoretical test, here we compare the sensitivity of filament length as a cosmological probe in comparison to the halo mass function. To determine the halo mass function, we employ the Friends-of-Friends (FOF) method. Here we use the number density of massive haloes with two mass bins with mass in the range of $\log_{10}(M/M_{\odot})\in[13.5,15]$ and $[13.8,15]$, respectively. Similarly, for the filament length function, we utilize the fitting parameters $b$ and $n_{0}$. We define a statistical function as follows:

Where $\mathbf{p}$ represents the parameter of the statistical quantity, denoting the number density in given mass bin for the halo mass function and the fitting parameters $b$ and $n_{0}$ for the filament length function. For our analysis, we adopt the Planck 2015 cosmology as the fiducial geometric background. We utilize 40 simulations from the $\Lambda$CDM dataset to compute the covariance matrix and an additional 5 simulations from the $\Lambda$CDM dataset to calculate $\mathbf{p}_{\text{fiducial}}$ for both the halo mass function and the filament length function. Simultaneously, we employ simulations from the $C1$ dataset to calculate $\mathbf{p}_{\text{target}}$.

In Fig. 7 we compare the ability of constraining cosmological parameters using halo mass function and filament length distribution. The solid lines are for halo mass function with halo mass within [13.5, 15] and the dashed lines are for halo in a higher mass bin of [13.8, 15]. It is found that in general the halo mass function shows stronger constraints on the cosmological parameters of $\sigma_{8}$ (left panel) and $\Omega_{m}$ in the right panel. However, the filament length function demonstrates stronger sensitivity to $n_{s}$ than the halos within the massive bin. (dashed line in the right panel).

We note that this comparison could just serve as a qualitative consideration. In practice, it is difficult to derive halo mass although it could be done with either x-ray or weak lensing measurements. On the other hand, as mentioned before, one has to use observed galaxies to extract filament. We leave more quantitative comparison to future work.

In this section we investigate the evolutionary trend in the distribution of filament lengths and its sensitivity to cosmological parameters. To conduct this analysis, we use 5 redshift snapshots in our numerical simulation: $z=3,2,1,0.5,0$. Utilizing our filament pipeline, we obtain the filament length distributions for the Planck15 cosmology at each snapshot, with $z=0$ serving as a reference. Fig.8 shows the redshift evolution of the filament length distribution, and it reveals that as redshift decreases from z=3, the number density of filament of all length firstly rises until approximately $z\sim 1$, and from $z\sim 1$ to $z\sim 0$, the number of long filaments declines while short filaments keep increase. To provide a more intuitive understanding, we use $N(L)$ instead of $\Phi_{L}$, as $N(L)$ less than 1 lack physical meaning. Our findings are in contrast with the previous study by Cautun et al. (2014), who reported a predominance of long filaments and fewer short ones at low redshifts. Our results are more consistent with Galárraga-Espinosa et al. (2024), despite their use of galaxies as tracers and a different definition of individual filament. The disagreement with Cautun et al. (2014) is mainly because they focused only on bifurcation and curvature, neglecting instances where peaks manifest and disrupt the filament, while we use a more broader definition. Our results reveals a three-step process in filament formation: 1) A first rise in filament formation with time, with rapid formation of longer filament by connections of shorter filaments. 2) After $z\sim 2$, long filaments reach peak efficiency in their formation, leading to a dynamical balance between the disruption and formation of long filament. 3) At lower redshift, increased matter clustering causes some long filaments to fragment, making shorter filaments more prevalent.

Additionally, we also find that Eq. 5 can effectively capture the redshift dependence of the filament length function with the following forms,

The relationship between $\beta$ and redshift is illustrated in the left panel of Fig.9. Specifically, the power-law index within the filament length function demonstrates a power-law correlation with redshift. By employing five redshift values to fit the power-law index and redshift relationship (namely $z=3.0,2.0,1.0,0.5,0.0$) in black, we successfully predict the behavior for another set of nine redshift values (namely $z=2.2,1.8,1.5,1.2,0.8,0.6,0.3,0.2,0.1$) in blue. This suggests that the fitting function effectively describes the power-law index $\beta$ and its relation to redshift, as shown in Fig.9. Furthermore, we ascertain that the relationship between $n_{1}$ and redshift can be adequately approximated by a cubic function, which is shown in the right panel of Fig. 9. Similarly, it is seen that the parameters obtained from five redshift can well predict the parameters at other redshifts.

Moreover, we observe that an increase in the parameter $\beta$ corresponds to a steep filament length function, indicating a relatively pronounced difference between the numbers of long and short filaments. This suggests a transition from extensive filaments to shorter ones. Furthermore, an increase in the parameter $n_{1}$ results not only in shorter filaments but also in more extensive ones, thereby increasing the overall presence of filaments. By combining these two relationships, we depict a simple scenario for the formation and evolution of filament in the following.

Under the assumption that long filaments are formed through the merging of shorter ones, we propose a simple picture to describe the evolution of filament. Basically, the formation of short filament is mainly propelled by the growth of density field. Concurrently, short filaments contribute to the maturation of long filaments, whereas localized density enhancements prompt the fragmentation of long filaments into shorter segments. Fig. 10 illustrates our scenario on the formation of both short and long filaments. The (F) process signifies formation of filament arising from the density perturbations, while (B) represents the breakdown of long filament into shorter ones. Additionally, (G) denotes the growth of long filament by merger/connection of short filaments.

Before $z\sim 2$, the $\beta$ value remains below 1, indicating a relatively minor difference between the numbers of long and short filaments. During this period, the dominant processes are filament formation (F) and growth (G), while filament breakdown processes (B) are negligible. As the redshift approaches $z=2$, the $\beta$ value increases towards 1, signifying that short filaments become more predominant than long filaments. The increasing slope of $n_{1}$ suggests a significant enhancement in the filament formation process (F). At approximately $z=1.5$, the slope of $n_{1}$ reaches its peak, indicating the greatest disparity between the growth process (G) and the breakdown process (B), signifying the maximum efficiency of formation of long filaments. Subsequently, the growth process (G) diminishes, while the breakdown process (B) becomes dominant, resulting in the fragmentation of long filaments and the prevalence of shorter ones.

In the above analysis, we have investigated the impacts of different cosmology parameters on the filament length distribution using simulations in which the initial power spectrum is produced from the CLASS code (Blas et al., 2011). As filament is a distinct structure on large scale, one naive thought is that the number of filaments with a given length could be determined by the power at a given scale. It is then interesting to test how a change of the initial power spectrum at a specific scale will impact the filament length distribution, and inspection of its impact will give more insight of how filament of given length is formed. Here we carry out an experiment by artificially changing the amplitude of the initial power spectrum at different scales with $k\in[0.02,0.1,0.2,0.3,0.4]$. In Fig.11 we show an example where we change the amplitude of the initial power spectrum with k in the range $k\in[0.02,0.1]$. Note that the black line with scale=1 is the standard initial power spectrum, for other lines the power spectrum is increased by a factor given by the parameter scale.

In Fig. 12, we show the resulting final non-linear power spectrum at z=0 and the filament length distribution for two cases, with the left panels for changing the power on large scale with $k\in[0.02,0.1]$ and the right panels for changing power on small scale with $k\in[0.3,0.4]$. The upper panels show the power spectrum at z=0. It is seen from the left panel that increasing the initial power on large scale, $k\in[0.02,0.1]$, results in strong increase of non-linear power spectrum on small scales. More interesting is that the boost of power spectrum is nearly constant at $k>0.1{\rm Mpc}^{-1}h$, as can be seen from the lower panel where the deviation of power spectrum respect to the original power spectrum (scale=1) is shown. Similar results on the power spectrum are also seen for the right panel with $k\in[0.3,0.4]$. The lower panels of Fig. 12 show the filament length distribution. It is found that increasing the power on large scales leads to formation of more longer filaments. This is because more short filaments are now connected to form long ones, so the total number of filaments, mainly dominated by short filaments, is decreasing. However, the lower right panel shows that increasing the power on small scales has the opposite effect, it leads to formation of more short filaments and less number of long filaments.

The above results indicate that the effect of power spectrum on the length distribution of filament is complicated. Overall, there is a dynamical equilibrium between the formation of short filament, formation of long filament by connecting short filaments and break of long filament into short ones. Unfortunately, there is lack of theory or model to describe how the length of filament is determined by the power spectrum. Our experimental simulation could serve as a starting point for future development of such a theory or model.