Abstract

Standard analyses of large-scale cosmic structure rely primarily on two-point correlation functions $\xi \left(r\right)$ and power spectra $P\left(k\right)$ . Distinct networks can share identical two-point statistics while differing fundamentally in graph-theoretic properties. Here, we analyze the topology of the cosmic web using join distributions $P\left(k\right)$ , which count independent connections (joins of order $k$ ) within the percolating filament graph and provide a direct, non-metric probe of effective dimensional organization $D_{\text{eff}}$ . We apply join analysis to observational galaxy surveys spanning a broad range of redshifts $0.1\le z\le 5.0$ . Across five statistically independent sky regions (each 500³ Mpc³, ~10⁶ galaxies), the observed $P\left(k\right)$ are statistically uniform (${\chi }^2/\nu <1.1$ ) and quantitatively fit by a Galton-Watson branching process with offspring probability generating function $f\left(s\right)=\left(1-p\right)+p{s}^m$ , mean offspring $m=1.054$ , yielding $D_{\text{eff}}={\log }_{2}m\approx 1.054$ . In contrast, the LCDM model predicts and $N$ -body simulations exhibit a broad $P\left(k\right)\propto k^{-\alpha }$ ($1.5\le \alpha \le 2.0$ ) with substantial high-order tails ($k=5-8$ ) from hierarchical halo mergers. These high-order joins are strongly suppressed (<1% observed vs. ~30% predicted) or absent in the data, yielding a $>25\sigma$ discrepancy (${\chi }^2>1700$ ) that cannot be resolved by parameter tuning (${\Omega }_m,{\sigma }_8$ ) or baryonic feedback. Three further independent empirical probes reinforce this conclusion: i) Box-counting analyses across redshift slices spanning ~12 billion years yield effective dimension $D_{\text{eff}}=-{\lim }_{ϵ\to 0}\log N\left(ϵ\right)/\log ϵ$ and density $\rho$ converging gradually toward 1.0 from above, $D_{\text{eff}}\left(z\right)\approx 1+0.02{\left(1+z\right)}^{0.1}$ . LCDM simulations do not exhibit comparable convergence. ii) Observational data show matter and voids percolate at the same threshold ${\rho }_c\approx 0.20\overline{\rho }$ (within 2% - 3%), whereas LCDM generically exhibits separated thresholds ${\rho }_c^{\text{matter}}\approx 0.22\overline{\rho }$, ${\rho }_c^{\text{void}}\approx 0.15\overline{\rho }$. iii) Finally, matter and voids display identical scale-free power-law behavior $N\left(>R\right)\propto R^{-\gamma }$ with $\gamma \approx 2.89$ ($D=3-\gamma \approx 0.11$ for void filling factor) in the scaling regime 1 < R < 20 Mpc∙h⁻¹, indicating phase symmetry consistent with critical phenomena.

Keywords

Cosmic Web Topology, Join Distributions, Percolation Analysis, Fractal Dimension, Dimensional Attractors, ΛCDM Discrepancy, Branching Process, Box-Counting Dimension, Matter-Void Symmetry, Scale-Free Scaling

Share and Cite:

1. Introduction

Observations over the past several decades have revealed that matter in the Universe is organized into a complex network of filaments, nodes, sheets, and voids, commonly referred to as the cosmic web (see simulations in Figure 1 and Figure 2; for comprehensive reviews, see [1]-[3]). This filamentary architecture spans scales from megaparsecs to the observable horizon, with galaxies tracing the luminous backbone of underlying dark matter distributions [4]. Within the LCDM paradigm, this structure emerges from gravitational instability acting on nearly Gaussian primordial density fluctuations $\delta \left(k\right)~\mathcal{N}\left(0,P\left(k\right)\right)$ with a near-scale-invariant power spectrum $P\left(k\right)\propto k^{n_s-4}$ ($n_s\approx 0.96$ ; [5]), with cold dark matter driving hierarchical clustering via the Zel’dovich approximation

Figure 1. A simulation of the large-scale galaxy distribution illustrating the filamentary cosmic web over 300 million years as modeled by IllustrisTNG. Galaxies (gold) are surrounded by blown-off shocked gas (white). Image credit: TNG Collaboration.

Figure 2. Enlarged view of the cosmic web as it appears in simulations based on the standard model of cosmology, emphasizing filament connectivity and void structure. Image credit: V. Springle, Max-Planck Institut für Astrophysik, Garching bei München.

$x\left(q,t\right)=q-D\left(t\right)\nabla \Psi \left(q\right)$ (where $D\left(t\right)\propto a\left(t\right)$ is the linear growth factor, $\Psi$ the displacement potential, and $a\left(t\right)$ the scale factor; [6] [7]). This framework has achieved remarkable success in reproducing observations from the cosmic microwave background anisotropies ($C_{\ell }$ power spectrum; [5]) to large-scale galaxy clustering $\xi \left(r\right)\approx {\left(r/r_0\right)}^{-\gamma }$ ($\gamma \approx 1.8$ , $r_0\approx 5$ Mpc∙h⁻¹; [8]) and weak lensing shear correlations [9].

The statistical characterization of large-scale structure has traditionally relied on correlation functions $\xi \left(r\right)=\langle \delta \left(x\right)\delta \left(x+r\right)\rangle$ and power spectra $P\left(k\right)={\displaystyle \int \xi \left(r\right){\text{e}}^{-ik\cdot r}{\text{d}}^{3}r}$ , which quantify clustering strength and scale dependence and have been instrumental in constraining cosmological parameters (${\Omega }_m\approx 0.31$ , ${\sigma }_8\approx 0.81$ ; [5]). However, these measures are fundamentally metric: they depend on distances, amplitudes, and scaling relations but do not directly encode topology, that is, how structures are connected (via the filament graph $G=\left(V,E\right)$ ), the number of independent paths (Betti number $b_1$ ), or how connectivity is organized across scales (persistence diagrams; [10]-[13]). Recent studies highlight that metric statistics alone cannot distinguish between tree-like ($b_1\approx 0$ ) and looped ($b_1\gg 1$ ) configurations with similar $\xi \left(r\right)$ [14] [15].

Topology matters. In graph-theoretic terms, two networks can share identical two-point statistics (e.g., the same degree distribution $P\left(\text{deg}\left(v\right)\right)$ or edge-length correlations $\langle l_e\rangle$ ) yet differ profoundly in connectivity, redundancy (excess degree $\Delta \mathrm{deg}=\mathrm{deg}\left(v\right)-2$ ), and global organization (genus $g=1-\chi /2$ for surfaces; [16]). A tree-like structure (acyclic, $b_1=0$ , $\nu =0$ ) and a highly interconnected mesh (high cyclomatic number $\nu =e-v+c>0$ ) may exhibit similar clustering amplitudes while representing radically different physical systems, as demonstrated in synthetic cosmic web models [17] [18]. In the cosmic web context, visual inspection of galaxy surveys reveals a richly connected network containing junctions (high-degree vertices), loops ($b_1>0$ ), and alternative paths, particularly within the percolating structure that spans survey volumes via the giant component at $p>p_c\approx 0.2$ [19]. Quantifying this connectivity requires observables beyond traditional clustering statistics, such as persistent homology tracking $b_i\left(R\right)$ across scales $R$ [20] [21].

Several topological descriptors have been explored in cosmology, including genus statistics $g\left(\nu \right)\propto H_1\left(\delta >\nu \right)$ [22] [23], Minkowski functionals $V_i={\displaystyle \int {\mu }_i\left(x\right){\text{d}}^{3}x}$ ($i=0\text{-}3$ for volume, surface area, mean curvature, and Euler characteristic [24] [25]), and percolation analyses with critical exponents $\beta \approx 0.41$ , $\nu \approx 0.88$ from finite-size scaling in 3D lattices [26] [27]. These provide valuable insights into morphology and connectivity but often integrate over information or are sensitive primarily to global transitions (e.g., ${\rho }_c$ ; [28]). In this work, we focus on join distributions $P\left(k\right)$ , which count independent connections (joins of order $k$ ) within connected structures and directly measure topological redundancy via local excess connectivity $k=\text{deg}\left(v\right)-2$ ; see [15] [28] [29]. Joins provide a natural probe of effective dimensional organization: tree-like networks exhibit minimal joins ($k=0$ , $D_{\text{eff}}=1$ ), while increasingly interconnected structures exhibit higher $k$ as shown in [1] [24].

In the LCDM model, hierarchical structure formation is formalized via the extended Press-Schechter (excursion-set) theory [30] [31], where the first-crossing distribution $f\left(\sigma ,m\right)=\frac{{\delta }_c}{\sqrt{2\pi }\sigma }\text{exp}\left[\frac{{\delta }_c^2}{2{\sigma }^2\left(m\right)}\right]\left|\frac{\text{d}\sigma }{\text{d}m}\right|$ yields merger trees with multi-progenitor branching $n_p\left(m_1,\cdots ,m_n|M,z\right)\propto {\displaystyle \prod f\left({\sigma }_i,m_i\right)}$ , leading to a broad join-order distribution $P\left(k\right)\propto k^{-\alpha }$ ($\alpha \approx 1.5\text{-}2.0$ due to cumulative ancestral paths in filamentary inheritance [18] [32]). This yields a percolating network with high cyclomatic complexity $\nu ~L^3/{\xi }^3\gg 1$ ($L$ box size, $\xi$ correlation length ≈100 Mpc∙h⁻¹), incompatible with near-tree-like topologies ($b_1\approx 0$ ; [33]). Recent simulations confirm this prediction but highlight discrepancies in void-matter symmetry [34].

In this paper, we demonstrate that LCDM makes a clear, testable topological prediction for $P\left(k\right)$ that fails against observational data ($>25\sigma$ ; see Section 3.4). We further show that four independent empirical probes converge on near-one-dimensional organization ($D_{\text{eff}}\approx 1$ ) within a three-dimensional embedding, echoing fractal analyses suggesting non-homogeneous evolution [35] [36]. This identifies a topological constraint on structure formation not captured by standard clustering statistics and not reproduced by LCDM, aligning with emerging tensions in cosmic topology [37].

2. Data and Observational Material

2.1. Galaxy Surveys and Sky Coverage

Our analysis is based on large-scale galaxy survey data spanning a wide range of redshifts and sky coverage, drawing from the Sloan Digital Sky Survey (SDSS [38]) and Baryon Oscillation Spectroscopic Survey (BOSS [39]) for low-to-intermediate $z$ , supplemented by higher-$z$ data from eBOSS [40]. To assess statistical robustness and cosmic variance, the sky is divided into five independent regions, each a cube of side $L=500$ Mpc containing $N_g~{10}^6$ galaxies with mean density $\overline{n}=N_g/L^3\approx {10}^{-2}{\text{Mpc}}^{-3}$ (comoving). These regions are spatially disjoint ($\left|r_i-r_j\right|>L$ ) and statistically independent ($\langle {\delta }_i{\delta }_j\rangle =0$ ), allowing tests of whether observed topological features are localized anomalies or global properties of the cosmic web, mitigating cosmic variance estimated at ${\sigma }_v\approx 5\text{\%}$ per region according to [41].

The surveys provide sufficient volume $V=L^3\approx 1.25\times {10}^8{\text{Mpc}}^3$ and sampling density $\bar{n}$ to robustly identify percolating structures over a range of density thresholds ${\rho }_{\text{th}}$ , with completeness >95% after selection function corrections [42]. Standard preprocessing includes survey mask correction (inverse weighting $w\left(\widehat{n}\right)$ to account for angular completeness); Data Release notes, DR12 as in [43] with selection-function weighting $S\left(z\right)\propto n\left(z\right)b\left(z\right)$ (where $b\left(z\right)$ is bias evolution [44]; Beutler et al., 2017), and redshift-space distortion mitigation via distant-observer approximation $s=r+v\cdot \widehat{z}/H\left(z\right)$ with Kaiser factor $f\approx {\Omega }_m^{0.55}$ [44] [45]. Details are documented in Notes A-E of Chapter 5 of Origins [46], reproduced in the Supplementary Methods.

2.2. Density Field Construction

Galaxy distributions are converted into continuous density fields $\rho \left(x\right)=\overline{\rho }\left[1+\delta \left(x\right)\right]$ using standard smoothing and thresholding techniques, following the Cloud-in-Cell (CIC) interpolation on a 512³ grid [47]. Smoothing employs a Gaussian kernel $W\left(\left|x\right|\right)={\left(2\pi {\sigma }^2\right)}^{-3/2}\text{exp}\left(-{\left|x\right|}^2/2{\sigma }^2\right)$ with $\sigma =2$ Mpc h⁻¹ to suppress shot noise while preserving filamentary features on $\gtrsim 5$ Mpc scales [48], yielding ${\delta }_s\left(x\right)={\displaystyle \int \delta \left(y\right)W\left(\left|x-y\right|\right){\text{d}}^{3}y}$ . While choices of $\sigma$ and threshold ${\rho }_{\text{th}}=\bar{\rho }\left(1+{\delta }_{\text{th}}\right)$ influence local morphology (e.g., filament thickness; [49]), percolation ($p\left(\rho >{\rho }_c\right)=0.5$ ) and join statistics are robust across $\text{Δ}\sigma /\sigma <25\text{\%}$ , $\text{Δ}{\delta }_{\text{th}}<0.1$ , as verified in robustness tests (Note A; compare with [50]). All join analyses are performed only after crossing the percolation threshold ${\rho }_c$ , ensuring global topology (spanning cluster size $S~V$ ) is meaningful and free from finite-size effects [51].

2.3. LCDM Simulations

For comparison, we analyze LCDM $N$ -body simulations matched in volume $V=L^3$ , resolution $N_p~{10}^9$ particles (${\overline{n}}_p={200}^3{\text{Mpc}}^{-3}$, particle mass $m_p\approx {10}^9{M}_{\odot }{h}^{-1}$ ), and selection effects to the observational data, using suites like MillenniumTNG [52] and IllustrisTNG [53] with Planck cosmology (${\Omega }_m=0.3089$ , ${\sigma }_8=0.8159$ ). These incorporate hierarchical dark-matter halo formation and mergers via the merger rate $\frac{{\text{d}}^{2}N}{\text{d}{m}_1\text{d}{m}_2\text{d}z}=\sqrt{\frac{2}{\pi }}\frac{{\delta }_c\left(z\right)\left(m_1+m_2\right)}{{\sigma }^2\left(m_1\right)\sqrt{{\sigma }^2\left(m_2\right)-{\sigma }^2\left(m_1\right)}}\text{exp}\left[\frac{{\delta }_c^2\left(z\right){\left(m_1+m_2\right)}^2}{2{\sigma }^2\left(m_2\right)}\right]\left|\frac{\text{d}{\sigma }^2}{\text{d}{m}_1}\right|$ as in [31] extended to semi-analytic models in [54]. The same analysis pipeline that is, density field construction ${\rho }_p\left(x\right)$ , percolation identification (Hoshen-Kopelman algorithm with 26-connectivity; [55]), and join counting is applied to simulations and data to ensure consistency [4] [56].

3. Join Distribution and Topological Measures

Methods:

This section describes the construction of the density field, the identification of percolating structures, the definition and counting of joins, the box-counting analysis used to estimate effective dimension, and the estimation of uncertainties. All methods are applied identically to observational data and to ΛCDM simulations unless otherwise stated.

1) Density Field and Smoothing: Embed galaxies in comoving volume; smooth with Gaussian $\sigma =2$ Mpc∙h⁻¹. Robust to ±25%variations (Note A; [49]).

2) Percolation: Threshold at ${\rho }_{\text{th}}$ ; identify spanning cluster via Hoshen-Kopelman ($p_c\approx 0.5$ sites); finite-size scaling validation [51].

3) Joins: DisPerSE skeletons [15]; $k=\text{deg}\left(v\right)-2$ at nodes $\text{deg}\left(v\right)\ge 3$ (Note E). Matched to LCDM via MillenniumTNG/IllustrisTNG; ${\chi }^2={\displaystyle \sum {\left(O_i-E_i\right)}^2/E_i}$ [57].

4) Box-Counting: $\log N\left(ϵ\right)=-D\log ϵ+C$ ; 10 - 12 logarithmic sizes ($0.5<\epsilon <20$ Mpc∙h⁻¹), $R^2>0.98$ ; bootstrap errors (×100 resamples; Note C; [58]).

5) Errors: Dispersion across skies ($\sigma =\sqrt{{\displaystyle \sum {\left(\text{Δ}{x}_i\right)}^2/4}}$ ); cosmic variance via simulations (Note E; [41]).

3.1. Definition of Joins

The filamentary skeleton extracted with the Discrete Persistent Structures Extractor (DisPerSE) is represented as a graph $G=\left(V,E\right)$ , where the vertex set V contains critical points, filament endpoints, and junctions, and the edge set E contains the filament segments connecting them. DisPerSE constructs such skeletons from the Morse-Smale complex of the density field and uses persistence to suppress features that are not robust against noise [15]. Related topological analyses of the cosmic web use persistent homology and Betti numbers to characterize connectivity on a global rather than purely local level [20] [21].

In the present work, a join is defined operationally as excess local branching at a filament junction. If $d_v=\text{deg}\left(v\right)$ is the graph degree of a vertex $v$ , then a junction is a vertex with $d_v\ge 3$ , and its join order is defined by $k_v=d_v-2$ ,$d_v\ge 3$ . This definition measures the local excess branching relative to an unbranched filament, whose internal vertices have degree two. Thus, a degree-three junction has $k_v=1$ , a degree-four junction has $k_v=2$ , and so forth. The join order is therefore a local connectivity statistic; it should not be identified with the number of independent cycles or with the first Betti number of the graph.

The join-order distribution over the set of junctions $V_J=\left\{v\in V:d_v\ge 3\right\}$ may be written as $P_{Join\left(k\right)}=\left(\frac{1}{N_J}\right){\displaystyle {\sum }_{v\in V_J}1\text{l}\left(k_v=k\right)}$ , $N_J=\left|V_J\right|$ . This distribution quantifies how strongly the skeleton branches locally and is insensitive to the metric lengths of the individual filament segments once the skeleton has been constructed. It can therefore be compared between observational and simulated skeletons extracted with the same persistence, smoothing, and segment-selection criteria [15].

Local branching must be distinguished from global topological redundancy. For a finite graph with E edges, $V$ vertices, and $C$ connected components, the number of independent graph cycles is the cyclomatic number, equivalently the first Betti number ${\beta }_1$ of the graph: ${\beta }_1=E-V+c$ .

For a single connected component $C=1$ , and hence ${\beta }_1=E-V+1$ . This quantity is global: it counts independent cycles after the full connectivity of the graph is known. Persistent-homology analyses extend this idea by following the birth and death of one-dimensional homology classes as the filtration scale changes [20] [21]. In the language of discrete Morse theory, the skeleton and its topological features may likewise be related to the critical-cell structure of the underlying complex [59].

The distinction is essential because branching does not imply a loop. A tree may contain many vertices of degree three or higher and therefore many non-zero join orders, while still satisfying ${\beta }_1=0$ . For example, a Y-shaped tree has one degree-three junction and hence $k_v=1$ but $E=3$ $V=4$ , and $C=1$ , giving ${\beta }_1=0$ . Conversely, the presence of ${\beta }_1>0$ establishes closed-loop redundancy, which cannot be inferred from the local degree distribution alone.

Accordingly, we use join statistics and global topology as complementary, rather than interchangeable, descriptors. The join distribution $P_{Join\left(k\right)}$ measures local branching complexity, whereas ${\beta }_1$ (or its persistent-homology generalization) measures global loop topology. A network may therefore be strongly branched but acyclic, weakly branched but looped, or possess both substantial branching and substantial cycle rank. This separation permits a cleaner comparison between the observed cosmic web and ΛCDM simulations [60]: differences in join order test the local junction structure, while differences in ${\beta }_1$ or persistent Betti numbers test global topological redundancy [20] [21].

3.2. Joins and Percolations

Positioning Join analysis is meaningful only for percolating structures. Below the percolation threshold $p_c$ (bond/site occupation fraction ≈ 0.2488 for 3D bond percolation; [27]), the system consists of disconnected clusters (finite $S\left(p<p_c\right)\ll V$ ), and join counts reflect local properties rather than global topology ${\nu }_{\text{global}}\approx 0$ . We restrict analyses to the percolating matter structure ($p>p_c$ , giant component $S~V$ ) and do not analyze the non-percolating regime, using finite-size scaling $S\left(p\right)/V~\mid p-p_c{\mid }^{\beta }{L}^{\gamma /\nu }$ near criticality ($\beta \approx 0.41$ , $\gamma \approx 1.80$ , $\nu \approx 0.88$ in 3D continuum; [51] [61] [62]). Thresholding is at ${\delta }_{\text{th}}=0.2$ , consistent with void finder algorithms (e.g., ZOBOV; [63] [64]).

3.3. LCDM Predictions from Hierarchical Merging

In LCDM, structure formation proceeds through hierarchical mergers of dark-matter halos, formalized by the conditional mass function $n\left(mM,z\right)=\frac{m}{{\sigma }^2\left(m\right)}\frac{\text{d}\sigma }{\text{d}m}\exp \left[\frac{{\left({\delta }_c\left(z\right)-{\delta }_c\left(z_i\right)\right)}^2}{2\left({\sigma }^2\left(m\right)-{\sigma }^2\left(M\right)\right)}\right]$ (extended Press-Schechter; [30]). As small halos ($m\ll M$ ) merge into larger systems ($M$ ), filaments inherit multiple ancestral connections, yielding merger trees with progenitor distribution $P\left(n_p\right)~\text{Poisson}\left(b\left(z\right)\right)$ , $b\left(z\right)~\text{ln}\left(1+z\right)$ as in [65]. This generically produces structures with many independent paths, resulting in a broad join-order spectrum $P\left(k\right)\propto k^{-\alpha }$ ($\alpha \approx 1.5\text{-}2.0$ , from power-law tails in branching as in [66]), including substantial high-order joins ($k=5-8$ , ${\displaystyle {\int }_5^{\infty }P\left(k\right)\text{d}k}~0.3$ ; [4] [18] [67] [68]). Simulations exhibit $\langle k\rangle \approx 3\text{-}4$ in percolating components, with $b_1\left(R\right)~R^{2.5}$ at megaparsec scales as explained in [32] [50].

This prediction is structural, following from hierarchical merging ($n_p>1$ ), and insensitive to baryonic physics (AGN feedback; [69]) or modest parameter variations ($\text{Δ}{\sigma }_8<0.1$ ; [70]).

3.4. Observed Join Distributions

Figure 3 shows join-order distributions $P\left(k\right)$ for the percolating matter structure

Figure 3. Percolating join-order distributions $P\left(k\right)$ measured across five independent sky regions (brown), compared with LCDM predictions (blue). Observed distributions are peaked at low $k$ and fit by a $D_{\text{eff}}=1.054$ branching process (tan). LCDM predicts high-$k$ tails from hierarchical mergers, strongly suppressed in data.

across five sky regions, compared to LCDM. Key features: i) uniformity across regions ($\frac{{\chi }^2}{\nu }<1.1$ , $p>0.05$ via KS test), indicating a stable global property robust to cosmic variance; ii) dominance by low-order joins ($k\le 4$ ,${\displaystyle {\int }_0^{4}P\left(k\right)\text{d}k}>0.99$ ), with $k\ge 5$ rare/absent (<0.01, consistent with supercritical branching); iii) quantitative fit to a critical Galton-Watson branching process $P\left(k\right)=\left(1-m\right){\delta }_{k,0}+m{\left(1-m\right)}^{k-1}$ ($m=1.054<2$ supercritical; mean fitted via maximum likelihood, ${\chi }^2/\text{dof}\approx 0.9$ ), $D_{\text{eff}}={\log }_{2}m\approx 1.054$ ($R^2>0.99$ , Kolmogorov-Smirnov $p>0.95$ across skies; compare with fractal fits in LSS).

The effective dimension $D_{\text{eff}}$ quantifies the self-similarity of the branching structure, distinguishing near-one-dimensional (tree-like) from higher-dimensional (mesh-like) topologies. It is derived from the similarity dimension of the Galton-Watson process [71] [72] as follows. Assume self-similarity: at each generation $g$ , the structure replicates $m$ times (mean offspring), with linear scale reduced by $r=1/2$ (dyadic scaling for resolution doubling in hierarchical grids). The mass $M_g$ (e.g., branch count or node number) scales as $M_g=M_0{m}^g$ . Resolution ${ϵ}_g={ϵ}_0{\left(1/2\right)}^g$ , so $g={\text{log}}_2\left(1/{ϵ}_g\right)$ . Thus, $M\left(ϵ\right)\propto m^{{\text{log}}_2\left(1/ϵ\right)}={\left(2^{{\log }_{2}m}\right)}^{{\text{log}}_2\left(1/ϵ\right)}={ϵ}^{-{\log }_{2}m}$ , yielding the fractal scaling $M\left(ϵ\right)\propto {ϵ}^{-D_{\text{eff}}}$ with $D_{\text{eff}}={\log }_{2}m=\ln m/\ln 2$ . For a straight line ($m=2$ segments at half-scale), $D_{\text{eff}}=1$ ; for plane-filling ($m=4$ ), $D_{\text{eff}}=2$ . The fit $m\approx 2.07$ ($D_{\text{eff}}\approx 1.054>1$ ) indicates slight bushiness beyond 1D, consistent with low-$k$ dominance and suppressed high-$k$ tails, differing from LCDM’s $\alpha \approx 1.8$ [18].

LCDM exhibits broad $P\left(k\right)$ with $\langle k\rangle >5$ , reflecting merger-induced redundancy ($b_1\gg 1$ as in [32]) (Table 1).

Consistency with Two-Point Correlation Functions and Power Spectra:

Topological results recover standard $\xi \left(r\right)\approx {\left(r/5\right)}^{-\gamma }$ ($\gamma =1.8$ ) and $P\left(k\right)\propto k^{n_s-4}$ ($n_s\approx 0.96$ ;[46] [74]). Two-point statistics quantify metric clustering but not global connectivity (distinct $G$ with same $P\left(\text{deg}\right)$ but different $\nu$ ; [12] [13] [75]). Discrepancies arise at the topological level ($b_i\left(R\right)$ ; as in [21]) (Table 2).

Table 1. Raw percolating join counts ($k=3\text{-}8$ ) for observed skies, 1.054D model, and LCDM simulations (averaged over 10 realizations).

Table 2. ${\chi }^2$ tests of discrepancies ($>25\sigma$ each). Each sky rejects LCDM individually, ruling out fluctuations [73].

4. Independent Empirical Constraints

4.1. Percolation of Matter and Voids

Percolation analysis probes global connectivity via the spanning cluster fraction $S\left(p\right)/V$ and void fraction $f_v=1-S\left(p\right)$ [76]. In data, matter ($\rho >{\rho }_{\text{th}}$ ) and voids ($\rho <-{\rho }_{\text{th}}$ ) percolate at the same threshold ${\rho }_c\approx 0.20\overline{\rho }$ (within 2% - 3%, $\text{Δ}{\rho }_c/\overline{\rho }<0.006$) across all skies, implying symmetric embedding in 3D with shared critical exponents ($\beta \approx \nu \approx 0.88$ as in [28] [29] [63]). This symmetry suggests a balanced phase transition, contrasting with Gaussian random fields where ${\rho }_c^{\text{void}}<{\rho }_c^{\text{matter}}$ due to skewness [77]. LCDM exhibits separated ${\rho }_c^{\text{matter}}\approx 0.22\overline{\rho }$ (log $p_c\approx 0.718$ ), ${\rho }_c^{\text{void}}\approx 0.15\overline{\rho }$ (log $p_v\approx -0.816$ ) due to non-Gaussian tails in $\delta \left(x\right)$ from gravitational amplification and hierarchical clustering [60] [61] [78] [79]. Same-threshold joint percolation is therefore not reproduced in LCDM and constitutes an additional empirical discriminator, as noted in void studies [80].

Joint percolation supports near-1D organization ($b_0\approx b_1$ , low $\chi \approx 1$ ; [63]) when combined with $P\left(k\right)$ and box-counting, implying a filament-dominated web with minimal loops [33].

4.2. Scale-Free Symmetry between Matter and Voids

Scaling analyses reveal power-law behavior $N\left(>R\right)={\displaystyle {\int }_R^{\infty }4\pi r^2\overline{n}\left(r\right)\text{d}r}\propto R^{-\gamma }$ with $\gamma \approx 2.89\pm 0.02$ for both phases over $1<\log R/\text{Mpc}<1.3$ (multiple decades; Figure 4(a) & Figure 4(b); fitted via least-squares on cumulative counts, $R^2>0.995$ ). This coincidence is highly non-trivial ($\text{Δ}\gamma <0.01$ , $p<{10}^{-3}$ via $t$ -test) and indicates scale-free organization and symmetry between phases, with void filling factor $f_v\approx R^{D-3}$ ($D=3-\gamma \approx 0.11$ as in [81] [82]). The exponent $\gamma \approx 2.89$ aligns with fractal void hierarchies in galaxy catalogs, as in [83] [84].

Figure 4. Scale-free power-law slopes $\gamma \approx 2.89$ for matter (a) and voids (b) in the non-percolating regime, indicating phase symmetry consistent with critical near-1D structure ($D\approx 1$ as in [85]).

Scaling probes the full density field; joins probe percolating connectivity. Power-law scaling without a characteristic scale ($\text{d}\gamma /\text{d}\log R=0$ ) hallmarks self-organized criticality; shared $\gamma =2.89$ is consistent with 1D embedding in 3D ($D_{\text{eff}}=3-\gamma +ϵ$ , $ϵ\ll 1$ [71] [82]), as in generalized void models [63] [86].

4.3. Box-Counting and Temporal Convergence

Box-counting estimates $D_{\text{eff}}=-{\lim }_{ϵ\to 0}\frac{\log N\left(ϵ\right)}{\log ϵ}$ , where $N\left(ϵ\right)$ is the number of occupied boxes of size $ϵ$ in the galaxy distribution, providing a robust measure of fractal scaling in LSS [58]. Table 3 summarizes results across $z$ -slices spanning ~12 Gyr (using DESI-like binning as in [87]), showing $D_{\text{eff}}$ and $\rho$ converge gradually to 1.0 from above, parametrized as $D_{\text{eff}}\left(z\right)\approx 1+0.02{\left(1+z\right)}^{0.1}$ (slow, persistent; least-squares fit $R^2=0.98$ ; [88]). Figures 5-7 illustrate attractor-like behavior $\dot{D}\propto -{\left(D-1\right)}^2$ (dynamical systems analogy as in [89]), absent

Table 3. Box-counting $D_{\text{eff}}$ and $\rho$ across $z$ -bins ($0.1\le z\le 5.0$ ); fits via least-squares on $\log N\left(ϵ\right)=-D\log ϵ+C$ over 10-12 scales [58].

Figure 5. Box-counting $D_{\text{eff}}$ vs. time since the Big Bang. Monotonic, gradual decrease toward 1.0D attractor as in [88].

Figure 6. Mean density evolution $\rho \left(z\right)$ , coupled to $D_{\text{eff}}$ convergence ($\rho \propto {\left(1+z\right)}^3{D}_{\text{eff}}^{-1}$ ).

in LCDM (no self-averaging $\text{Var}\left(\xi \left(r\right)\right)\to /0$ as in [35] [90]). This evolution suggests a transition from early fractal ($D\approx 1.2\text{-}1.3$ as shown in [91]) to late near-1D organization, contrasting LCDM’s persistent $D\approx 1.5\text{-}2.0$ as demonstrated by [92].

5. Synthesis: Four Independent Lines of Evidence

The Universe exhibits stable near-1.0D topology for ≥12 Gyr: i) $P\left(k\right)$ fit by subcritical branching $m=1.054$ ($D_{\text{eff}}\approx 1.05$ as in [15]); ii) $D_{\text{eff}},\rho \to {1.0}^{+}$ attractor as explained in [35] [88]; iii) joint percolation at ${\rho }_c\approx 0.20\overline{\rho }$ as demonstrated by

Figure 7. LCDM lacks a low-density attractor, showing persistent scale-dependent decrease $\rho \propto z^{\alpha }$ , $\alpha \approx 1$ [92].

[63]; iv) shared $\gamma =2.89$ scaling ($D=3-\gamma \approx 0.11$ voids as in [85]). These constrain viable theories beyond LCDM, e.g., via modified gravity altering merger kernels (e.g., $f\left(R\right)$ models; [93]) or primordial non-Gaussianity suppressing loops as explained in [94].

6. Discussion

LCDM reproduces geometric statistics (e.g., $\xi \left(r\right)$ fractal at $D\approx 1.2\text{-}1.3$ on small scales < 100 Mpc as explained in [35] [91]), but fails topologically ($b_1^{\text{obs}}\ll b_1^{\Lambda \text{CDM}}$ as shown in [18]), with recent JWST data exacerbating tensions in early structure growth [95]. Earlier topologies (e.g., 1D fractals in 3D, $D_{\text{loc}}=1$ , $D_{\text{glob}}=3$ as in [71]; Figure 8) may resemble observations on $R\gtrsim 300$ Mpc, where homogeneity emerges as observed in [96].

Figure 8. 1.0D fractal embedded in 3D space, illustrating short-scale fractality ($D<2$ ) and large-scale homogeneity ($D=3$ as in [71]).

We establish empirical constraints ($P\left(k\right),b_i\left(R\right),{\rho }_c$ ) for future models, emphasizing topology via persistent homology as in [3] [21] and urging reconfirmation with DESI/Euclid [87] [97].

7. Conclusion

The cosmic web exhibits persistent near-1D topology ($D_{\text{eff}}\approx 1$ ) over cosmic history. Joins reveal LCDM failure ($>25\sigma$ as in [50]); percolation, scaling, and box-counting support a 1D attractor as demonstrated [88]. Reconfirmation via JWST/Euclid is needed to assess LCDM viability amid growing tensions as pointed out in [37].

AI Collaboration Statement

AI (ChatGPT-4/5) assisted with graphics and exploration, all decisions by authors.

Data Availability

Public surveys (SDSS DR12 [43] and BOSS [39]). Simulations are cited.

Author Contributions

Conceptualization: S.A.K.; Methodology: S.A.K.; Software: S.M.; Validation: S.M.; Formal analysis: S.M.; Investigation: S.A.K. and S.M.; Resources: S.A.K. and S.M.; Data curation: S.A.K. and S.M.; Writing—original draft preparation: S.A.K.; Writing—review and editing: S.M.; Visualization: S.A.K.; Supervision: S.A.K.; Project administration: S.A.K.; Funding acquisition: S.A.K. All authors have read and agreed to the published version of the manuscript.

Supplementary Material

Note A: Filament Skeleton Extraction and Parameter Robustness

To extract filamentary skeletons from the galaxy distributions, we employ the Discrete Persistent Structures Extractor (DisPerSE; [15]), following the methodology of Malavasi et al. [48] for consistency with prior topological analyses of the cosmic web. This algorithm identifies persistent topological features in the density field by computing the discrete Morse-Smale complex [11] [98], yielding a graph $G=\left(V,E\right)$ where vertices $V$ represent critical points (e.g., maxima, saddles) and edges $E$ trace ridge-like filaments connecting them. Four key parameters govern the extraction:

• Persistence Threshold (PT): Set to $3\sigma$ significance, where $\sigma$ is the noise level in the smoothed density field ${\delta }_s\left(x\right)$ . This filters features with persistence $\lambda =\left|\delta \left(x_c\right)-\delta \left(x_s\right)\right|>3\sigma$ , ensuring robustness against shot noise [59].

• Smoothing Length ($\sigma$ ): Gaussian kernel width $\sigma =2$ Mpc∙h⁻¹, as $W\left(\left|x\right|\right)={\left(2\pi {\sigma }^2\right)}^{-3/2}\text{exp}\left(-{\left|x\right|}^2/2{\sigma }^2\right)$ , balancing resolution of sub-Mpc structures with suppression of Poisson fluctuations (variance $\text{Var}\left(\delta \right)\approx 1/\overline{n}{V}_{\text{cell}}$ as in [49]).

• Minimum Segment Length (MSL): 3 Mpc∙h⁻¹, pruning short spurious edges below this threshold to focus on megaparsec-scale filaments (compare with ridge extraction in Sousbie et al. [15]).

• Linking Length (LL): 1.5 Mpc∙h⁻¹, the maximum geodesic distance for connecting critical points, analogous to a friends-of-friends parameter in halo finding [4].

These values match those in Malavasi et al. [48] for direct comparability. To test robustness, we varied each parameter by ±25% (e.g., $\sigma \in \left[1.5,2.5\right]$ Mpc∙h⁻¹; $\text{PT}\in \left[2.25\sigma ,3.75\sigma \right]$ ) and recomputed join distributions $P\left(k\right)$ and percolation thresholds ${\rho }_c$ . Results show invariance: $\langle k\rangle$ varies by <5%, and ${\chi }^2$ discrepancies with LCDM remain $>20\sigma$ (bootstrap resampling over 100 trials as shown by Newman & Barkema [73]). This confirms that statistical features are insensitive to moderate parameter choices, aligning with convergence tests in cosmic skeleton analyses [17].

Note B: Construction of RealSky Datasets from SDSS

The RealSky datasets (RealSky and RealSky2) are constructed from spectroscopic galaxy samples in the Sloan Digital Sky Survey (SDSS; [38] [43]). RealSky utilizes the Data Release 7 (DR7) North Galactic Cap main sample at $z<0.25$ , with luminosity cuts $M_r<-20.5$ (monopole) to ensure volume-limited completeness [99]. RealSky2 draws from DR12 BOSS LOWZ+CMASS regions at 0.15 < z < 0.7, incorporating color and magnitude selections: LOWZ ($r<19.6$ , 0.15 < z < 0.4) and CMASS (17.5 < i < 19.9, 0.4 < z < 0.7) for balanced redshift coverage [39] [42].

Preprocessing includes fiber collision corrections via nearest-neighbor upweighting $w_{ij}=1+{\delta }_{ij}$ (where ${\delta }_{ij}$ accounts for angular separations < 55" [100] and redshift-space distortion mitigation using the linear Kaiser model $\beta =f/b\approx 0.4$ [44] [45]). These diverse SDSS-based datasets test the robustness of filament skeletons and join counts $P\left(k\right)$ against survey footprint (angular completeness > 80%), depth (effective volume $V_{\text{eff}}~{10}^7{\text{Mpc}}^3$ ).

Note C: Fractal Dimension Estimation via 3D Box-Counting

The effective fractal dimension $D_{\text{eff}}$ is estimated using 3D box-counting on voxels classified as percolating or non-percolating structures, following standard multifractal analysis in large-scale structure [58]. For each structure, we overlay isotropic cubic boxes of varying sizes $ϵ\in \left[0.5,20\right]$ Mpc∙h⁻¹ (log-spaced, 10 - 12 sizes spanning > 1 decade) and count occupied boxes $N\left(ϵ\right)$ , where occupancy requires ≥1 galaxy or density excess $\delta >0$ . The dimension is extracted from the power-law fit:

$\log N\left(ϵ\right)=-D_{\text{eff}}\log ϵ+C,$

via least-squares regression, requiring $R^2>0.98$ for scaling validity as demonstrated in [13]. Error bars on $D_{\text{eff}}$ are computed via bootstrap resampling: for each structure, draw 100 random 80% voxel subsamples with replacement, refit the slope for each, and take the $\pm 1\sigma$ spread as uncertainty [101]. This yields ${\sigma }_D\approx 0.02$ per bin, capturing cosmic variance and finite-size effects (scaling regime $L/ϵ>10$ as in [27]). For percolating components, $D_{\text{eff}}$ converges to ≈1.01 at low $z$ , consistent with filament dominance [88] [91].

Note D: Mean Density and Standard Deviation Estimation

Mean density $\rho$ and standard deviation ${\sigma }_{\rho }$ are estimated by voxelizing the 3D galaxy distribution into cubic voxels of side (2 Mpc∙h⁻¹)³, matching the skeleton smoothing scale for consistency [49]. Voxels are classified as belonging to the percolating skeleton (giant component $S>0.5V$ ) or non-percolating (isolated) structures using the Hoshen-Kopelman algorithm with 26-connectivity [55]. For each class separately, galaxy counts per voxel $n_v$ are tallied, and statistics computed as:

$\rho =\frac{1}{N_v}{\displaystyle \sum }_v{n}_v/V_v,{\sigma }_{\rho }=\sqrt{\frac{1}{N_v-1}{\displaystyle \sum }_v{\left(n_v/V_v-\rho \right)}^2},$

where $N_v$ ∼250,000 - 300,000 for percolating (stable ${\sigma }_{\rho }/\rho <10\text{\%}$ ) and ∼30,000 - 50,000 for non-percolating structures (ensuring Poisson stability; $\lambda =\bar{n}{V}_v\gg 1$ see[42]). This approach enables direct comparison across observed and simulated datasets, revealing coupled evolution $\rho \left(z\right)\propto {\left(1+z\right)}^3{D}_{\text{eff}}^{-1}$ as in [41].

Note E: Fair Comparison of Join Distributions with ΛCDM Simulations

To ensure a rigorous comparison between observed skies and ΛCDM $N$ -body simulations with respect to join distributions $P\left(k\right)$ , we analyzed both the real galaxy distribution as in [50] and ΛCDM simulations using the identical filament-identification method as in [15], with matched number densities $\bar{n}\approx {10}^{-2}{\text{Mpc}}^{-3}$ and redshift bins 0.1 < z < 0.5. Skeletons were generated for the ΛCDM “sky” using the same procedure as [50], defining filament joins identically as intersections of three or more branches ($k=\text{deg}\left(v\right)-2\ge 1$ ) (see [15]). Our ΛCDM analysis incorporated multiple large-volume $N$ -body simulations, including MillenniumTNG [52] ($V={500}^3{\text{Mpc}}^3$ , $N_p={4320}^3$ ) and IllustrisTNG (as in [53] hydrodynamical runs TNG300), averaging over 10 realizations to suppress cosmic variance (${\sigma }_v<3\text{\%}$ [41]).

Join counts were extracted for each statistically independent sky region and compared directly to corresponding ΛCDM realizations using the ${\chi }^2$ statistic:

${\chi }^2={\displaystyle \sum }_k\frac{{\left(O_k-E_k\right)}^2}{E_k},$

where $O_k$ and $E_k$ are observed and expected counts for join order $k$ [57]. The extremely high ${\chi }^2>1700$ corresponds to a $~25\sigma$ rejection of the null hypothesis that ΛCDM reproduces the observed cosmic web join structure [73]. This strong discrepancy indicates that ΛCDM fails to account for the real Universe’s topological filament connectivity at extremely high statistical significance, persisting across parameter variations ($\text{Δ}{\sigma }_8<0.1$ [70] and feedback prescriptions [69]).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References