Percolation transition in entangled granular networks

TL;DR

This study investigates how entangled granular particles form networks that undergo a percolation transition, revealing the critical link threshold and the influence of particle geometry on network formation.

Contribution

It introduces a network science approach to analyze entanglement networks and proposes a continuum percolation model for rings that captures the transition.

Findings

Networks undergo a percolation transition with a critical link threshold.

Particle opening angle influences network structure and mechanical bonding.

A continuum percolation model effectively describes the transition.

Abstract

Highly nonconvex granular particles, such as staples and metal shavings, can form solid-like cohesive structures through geometric entanglement (interlocking). The network structure formed by this entanglement, however, remains largely unexplored. Here we utilize network science to investigate the entanglement networks of C-shaped granular particles under vibration through experiments and simulations. By analyzing key network properties, we demonstrate that these networks undergo a percolation transition as the number of links increases logarithmically over time; the entangled particles form a giant cluster when the number of links exceeds a critical threshold. We propose a continuum percolation model of rings that effectively describes the observed transition. Additionally, we find that particle's opening angle significantly affects mechanical bonding and, consequently, the network structure. This work highlights the potential of network-based approaches to study entangled materials, paving the way for advancements in applications ranging from mechanical metamaterials to entangled robot swarms.

Figures

a). We test nine sets of C-particles with $\theta$ ranging from $25^{\circ}$ to $135^{\circ}$, each with $N=4000$ particles. Initially, disentangled particles are poured in a conical container. The conical shape ensures a consistent geometry, regardless of the amount of particles. To induce particle entanglement, we shake the container with vertical sinusoidal vibration at 20 Hz and 1.55 mm amplitude for a duration of time $t$. Since particles’ positions and links are difficult to measure, we measure the sizes of mechanically bonded clusters. After stopping the vibration, we slowly lift the cl

a shows a snapshot of the entanglement network from a DEM simulation; the network exhibits relatively homogeneous node degrees (number of links, $k$) with no prominent hubs.

g shows that when $\langle k\rangle$ is near the threshold observed in Fig. 3a, b or 4b, both the C-particle networks and the ring CP model closely follow $n_{s}\sim s^{-2.19}$, suggesting that they fall into the standard 3D percolation universality class where $\tau=2.19$ [[43](https://arxiv.org/html/2509.00216v2#bib.bib43)]. Furthermore, $n_{s}$ of in-container clusters agrees with the ring CP model over various $\langle k\rangle$, not only at the threshold (Extended Data Fig. 6a–d). $n_{s}$ of lifted clusters also matches the ring CP model before a giant cluster forms, but deviates afterwar

b, giving $\eta_{\textrm{c}}=2.11$. This suggests that the mean degree at the percolation threshold is $\langle{k}\rangle_{\textrm{c}}=2.11$ when $N\to\infty$, substantially larger than $\langle{k}\rangle_{\textrm{c}}=1$ of ER random networks [[25](https://arxiv.org/html/2509.00216v2#bib.bib25)] and comparable to other CP models: $\langle{k}\rangle_{\textrm{c}}=2.74$ for spheres [[39](https://arxiv.org/html/2509.00216v2#bib.bib39)] and $\langle k\rangle_{\textrm{c}}=2.27$ for disks in 3D space [[40](https://arxiv.org/html/2509.00216v2#bib.bib40)]. This large $\langle{k}\rangle_{\textrm{c}}$ ca

a. To form a Hopf link, the blue ring should enclose either $P$ or $Q$, but not both. This can be conveniently seen in the $X$-$Y^{\prime}$ plane as shown in Extended Data Fig. 10b: For given $y$ and $\gamma$, $C_{2}$ can only be located within the shaded area, which is formed by two circles with radius $R$ centered at $P$ and $Q$. The shaded area can be found by subtracting the overlap area from the total area of two circles: $A_{\rm{shaded}}(y)=2\left[\pi R^{2}-A_{\rm{overlap}}(y)\right]$. The overlap area $A_{\rm{overlap}}(y)=2\left[R^{2}\theta(y)-x(y)h(y)\right]$, where $\theta(y)=\cos^{-1

a). We test nine sets of C-particles with $\theta$ ranging from $25^{\circ}$ to $135^{\circ}$, each with $N=4000$ particles. Initially, disentangled particles are poured in a conical container. The conical shape ensures a consistent geometry, regardless of the amount of particles. To induce particle entanglement, we shake the container with vertical sinusoidal vibration at 20 Hz and 1.55 mm amplitude for a duration of time $t$. Since particles’ positions and links are difficult to measure, we measure the sizes of mechanically bonded clusters. After stopping the vibration, we slowly lift the cl

a shows a snapshot of the entanglement network from a DEM simulation; the network exhibits relatively homogeneous node degrees (number of links, $k$) with no prominent hubs.

g shows that when $\langle k\rangle$ is near the threshold observed in Fig. 3a, b or 4b, both the C-particle networks and the ring CP model closely follow $n_{s}\sim s^{-2.19}$, suggesting that they fall into the standard 3D percolation universality class where $\tau=2.19$ [[43](https://arxiv.org/html/2509.00216v2#bib.bib43)]. Furthermore, $n_{s}$ of in-container clusters agrees with the ring CP model over various $\langle k\rangle$, not only at the threshold (Extended Data Fig. 6a–d). $n_{s}$ of lifted clusters also matches the ring CP model before a giant cluster forms, but deviates afterwar

b, giving $\eta_{\textrm{c}}=2.11$. This suggests that the mean degree at the percolation threshold is $\langle{k}\rangle_{\textrm{c}}=2.11$ when $N\to\infty$, substantially larger than $\langle{k}\rangle_{\textrm{c}}=1$ of ER random networks [[25](https://arxiv.org/html/2509.00216v2#bib.bib25)] and comparable to other CP models: $\langle{k}\rangle_{\textrm{c}}=2.74$ for spheres [[39](https://arxiv.org/html/2509.00216v2#bib.bib39)] and $\langle k\rangle_{\textrm{c}}=2.27$ for disks in 3D space [[40](https://arxiv.org/html/2509.00216v2#bib.bib40)]. This large $\langle{k}\rangle_{\textrm{c}}$ ca

Equations

$$P(k)=e^{-⟨k⟩}\frac{⟨k{⟩}^k}{k!},$$

$$⟨k⟩(t)=k_0+\alpha \ln (1+t/t_0),$$

$$S_i$$

$$\chi$$

$$\left(R^2-|\mathbf{b}{|}^2\right)(\mathbf{k}\cdot {\mathbf{r}}_{ij}{)}^2>|\mathbf{k}{|}^2{\left(\frac{|{\mathbf{r}}_{ij}{|}^2}{2}+\mathbf{b}\cdot {\mathbf{r}}_{ij}\right)}^2,$$

$${\tilde{v}}_{\mathrm{ex}}^{\mathrm{ring}}(\gamma )=2{\int }_0^R{A}_{\mathrm{shaded}}(y)\sin \gamma dy=\frac{32}{3}{R}^3\sin \gamma .$$

$$v_{\mathrm{ex}}^{\mathrm{ring}}={\int }_0^{2\pi }{\int }_0^{\pi }{\tilde{v}}_{\mathrm{ex}}^{\mathrm{ring}}(\gamma )\frac{\sin \gamma d\gamma d\varphi }{4\pi }=\frac{\pi }{3}{D}^3,$$

Full text

††preprint: APS/123-QED

I Introduction

Particle shape plays a crucial role in the collective behavior of granular materials. Unlike the well-studied spherical particles [1, 2], the behavior of nonspherical particles remains poorly understood. Slender nonconvex particles, in particular, exhibit distinctive mechanical properties due to their entanglement (interlocking). Entanglement manifests macroscopically as deformation-resistant geometric cohesion [3, 4, 5]. For example, entangled U-shaped [4], Z-shaped [6], star-shaped [7, 8, 9], or cross-shaped [10] particles can resist uniaxial compression and maintain freestanding columns. Moreover, S- or U-shaped particles can withstand tensile stresses [11] and be collectively lifted against gravity [12, 8]. These mechanical effects arising from entanglement have recently sparked interest across diverse fields, such as architecture [8], living systems [13, 14], soft robotics [15], and metamaterials [16, 17].

Despite this interest, the structural properties of entangled granular materials, which impact their mechanical behavior, have rarely been quantified. To address this gap, we propose a network-based approach that models entangled granular materials as networks (graphs), where nodes represent particles and edges represent topological links. While network-based approaches have been applied to force chains of convex grains [18], they have not yet been extended to the entangled structures of nonconvex grains. The network representation of entangled particles has been briefly introduced in simulation studies for microscopic particles, such as C-shaped colloids [19] and kinetoplast DNA [20, 21]. In these studies, two particles are considered entangled when they form a Hopf link—a topological configuration in which two loops are interlinked and cannot be separated without breaking. We adopt this definition to construct entanglement networks, enabling quantitative analysis using complex network theory.

Here, we investigate the entanglement networks of C-shaped granular particles (C-particles) evolving under vibration through experiments and simulations. The C shape is chosen for its adjustable entanglement capability by varying its opening angle while maintaining the same particle diameter, which allows for proper comparison between different particle batches. By lifting the entangled clusters, we demonstrate that topological links largely reflect the mechanical bonds under tension. Through measurements of key network properties, we reveal that the C-particle networks exhibit a percolation transition, which is effectively captured by our proposed continuum percolation (CP) model of rings. Additionally, we observe that the mean degree grows logarithmically with vibration time, mirroring slow relaxation in disordered systems. These results establish network-based analysis as a promising framework for elucidating the structure and dynamics of entangled granular matter.

II Experiment and Simulation Methods

We experimentally measure the cluster formation of steel C-particles under controlled vibration. Each particle has a toroidal shape with opening angle $\theta$, centerline diameter $D=9~\rm{mm}$, and thickness $d=1~\rm{mm}$ (Fig. 1a). We test nine sets of C-particles with $\theta$ ranging from $25^{\circ}$ to $135^{\circ}$, each with $N=4000$ particles. Initially, disentangled particles are poured in a conical container. The conical shape ensures a consistent geometry, regardless of the amount of particles. To induce particle entanglement, we shake the container with vertical sinusoidal vibration at 20 Hz and 1.55 mm amplitude for a duration of time $t$. Since particles’ positions and links are difficult to measure, we measure the sizes of mechanically bonded clusters. After stopping the vibration, we slowly lift the clusters one by one from the container (Fig. 1b,c) and weigh them. Because the measurement is destructive, each trial begins with a freshly disentangled sample. This procedure is repeated in 578 experimental trials at different $t$ and $\theta$ for sufficient statistics.

To measure the links between particles, we simulate the corresponding system using the discrete element method (DEM) (Fig. 1d). We test 10 sets of frictional C-particles, each comprising a chain of spheres, with $\theta$ ranging from $20^{\circ}$ to $151^{\circ}$. The particle thickness, vibration amplitude, peak acceleration, and time unit are set to match the experiment. From the positions and orientations of the particles, we determine whether they form Hopf links (Eq. 5, Extended Data Fig. 10) [19]. Although C-particles are not closed loops, this criterion well approximates their mechanical interlocking, especially when $\theta$ is small.

The entanglement networks are measured before (in-container) and after lifting the clusters in gravity. The results are explained by Monte Carlo (MC) simulations of the ring CP model. Experiment and simulation details are in Methods and Supplementary Information.

III Results

III.1 Evolution of the largest cluster

Under vibration, initially unentangled C-particles gradually entangle and form clusters, and eventually, the largest cluster dominates the system. In percolation theory, the relative size of this largest cluster, $S_{1}\equiv s_{1}/N$, serves as the order parameter, where $s_{i}$ represents the number of particles in the $i$th largest cluster, and $N$ is the total number of particles [22]. The experiments (Fig. 1e) and simulations (Fig. 1f,g) show that for small $\theta$, lifted $S_{1}(t)$ steadily grows and approaches 1, indicating percolation of entanglement. Additionally, the susceptibility, defined as the mean cluster size $\chi\equiv\sum_{i\neq 1}{s_{i}^{2}}/N$, peaks near the rapid growth of $S_{1}(t)$, i.e., the percolation threshold (critical point) as predicted by percolation theory (Fig. 1h,i) [23].

The growth rate of $S_{1}(t)$ depends on $\theta$. For $\theta\lesssim 70^{\circ}$, C-particles with larger $\theta$ entangle more easily during vibration, resulting in faster growth of $S_{1}$ at $t>0$ (Fig. 1e,f). However, for $\theta\gtrsim 70^{\circ}$, this trend reverses (Fig. 1g,j) because larger $\theta$ makes the links break more easily. When only a few tiny clusters detach from the giant cluster, $S_{1}$ remains close to 1, but when the giant cluster disintegrates during lifting, $S_{1}$ decreases dramatically. These two cases result in a bimodal $S_{1}$ distribution at large $\theta$ (Extended Data Fig. 1c,e).

The clusters become extremely fragile beyond a certain $\theta$. In experiments at $t=60$ s, giant clusters can be lifted stably for $\theta\leq 110^{\circ}$, but often disintegrate for $\theta\geq 115^{\circ}$. As a result, the ensemble-averaged $S_{1}$ at $t=60$ s plummets around $\theta=115^{\circ}$ (Extended Data Fig. 1f). This mechanical instability, emerging beyond a certain $\theta$, is also confirmed by simulations. When $\theta\gtrsim 150^{\circ}$, the clusters are too fragile to be lifted as giant clusters (Fig. 1g,j, Extended Data Fig. 1g).

This maximum $\theta$ for stability depends on particle thickness and surface friction. Ref. [24] has shown that even straight steel rods can form cohesive clusters through friction when the particle aspect ratio exceeds 100. In contrast, our C-particles are thicker, and thus their cohesion relies primarily on normal contact forces within topological links. If the particles were thinner, stronger frictional effects would raise the maximum $\theta$ for stability, reducing the fragile region in Fig. 1j.

III.2 Evolution of the degree distribution

To investigate the entangled structure, we model the granular system as an $N$-node network comprising multiple distinct clusters (connected components), where nodes represent particles and edges denote topological links. Figure 2a shows a snapshot of the entanglement network from a DEM simulation; the network exhibits relatively homogeneous node degrees (number of links, $k$) with no prominent hubs.

This degree homogeneity is evident in the degree distribution $P(k)$, a fundamental metric of network connectivity [25]. The observed $P(k)$ of the C-particle networks at various $t$ mostly agrees with the Poisson distribution (Fig. 2b,c):

[TABLE]

where $\langle k\rangle$ is the mean degree of each network. This agreement holds for all in-container networks (before lifting) (Fig. 2b) and stably lifted networks with small $\theta$ (Fig. 2c), except for lifted networks with large $\theta$ that undergo frequent cluster breakages (Extended Data Fig. 2).

The Poisson $P(k)$ is a feature of both Erdős–Rényi (ER) random networks [25] and continuum percolation (CP) models, also known as random geometric graphs (RGGs) [26]. A similar $P(k)$ is observed in kinetoplast DNA networks, quasi-2D assemblies of flexible rings [21, 20]. Despite differences in geometry and dynamics, C-particles and DNA rings share an underlying spatial randomness, which explains their similar $P(k)$ to that of CP models.

In the following sections, we show that other properties of C-particle networks differ from ER random networks but are similar to a CP model of rings.

Given that the Poisson distribution is solely determined by $\langle{k}\rangle$, knowing how $\langle{k}\rangle$ changes with $\theta$ and $t$ is crucial. Figure 2d shows that $\langle{k}\rangle(t)$ of the in-container networks can be fitted by

[TABLE]

where the fitting parameters $k_{0}$ and $t_{0}$ vary with $\theta$, while $\alpha$ is relatively independent of $\theta$ (Extended Data Fig. 3a–c). This logarithmic increase persists over a wide range of $t$, even long after $S_{1}(t)$ (Fig. 1f) has nearly saturated. The lifted networks also display logarithmic $\langle{k}\rangle(t)$, but with the curves shifting downward (Fig. 2d) because some links break during the lifting.

To our knowledge, this logarithmic evolution of entangled granular network has not been previously reported. The most relevant prior observation is the logarithmic density increase in vibrated spherical grains [27]. However, in our system, the degree exhibits a clearer and more consistent logarithmic trend than the density, which shows larger fluctuations (Extended Data Fig. 8).

Equation (2) is purely empirical, and its theoretical basis remains unclear. We speculate that insights may be gained by considering other glassy disordered systems, in which slow logarithmic relaxation (aging) is common. Examples include stress relaxation in polymers [28] and granular media [29], structural relaxation in DNA [30] and crumpled paper [31], as well as magnetization decay in spin glasses [32]. In these systems, slow relaxation typically arises from sequential transitions between metastable states separated by broadly distributed waiting times, often due to rugged energy landscapes [33, 34, 35]. These systems also exhibit history dependence (memory effects [36]. Our observation that lifted C-particle clusters retain their boundary shape (Fig. 1c) and degree distribution (Fig. 2b–d) may suggest a similar form of memory. Whether existing theories of glassy dynamics apply to entangled granular systems remains a puzzle.

III.3 Comparison with the ring percolation model

To describe the percolation transition in C-particle networks, we propose a CP model of infinitely thin rings. CP models, also known as random geometric graphs [26], have been utilized to study various disordered systems, including porous media, semiconductors, and wireless networks [37, 38]. In conventional CP models, objects such as spheres [39], disks [40], or rods [41] are randomly distributed in space and are considered connected when they overlap. In our model, randomly distributed rings in 3D are considered connected when entangled (Methods). Our MC simulations show that as ring density increases, the ring CP undergoes a percolation transition at a unique threshold and belongs to the same universality class as conventional CP models (Fig. 3,4).

Although C-particles are neither infinitely thin nor completely randomly located, their entanglement networks are well captured by the ring CP model (Fig. 3, 4). The in-container $S_{1}(\langle{k}\rangle)$ curves for various $\theta$ collapse nicely onto $S_{1}$ of the ring CP (Fig. 3a). By comparison, the lifted networks are more densely connected, i.e., have higher $\langle{k}\rangle$, than the ring CP networks with the same $S_{1}$. This deviation is stronger at larger $\theta$, because large openings facilitate the disintegration of the giant cluster, while the decrease in $\langle{k}\rangle$ remains limited owing to the tightly connected core. The percolation threshold can be roughly estimated at $S_{1}(k_{1/2})=1/2$, and the measured $k_{1/2}$ increases linearly with $\theta$ (Fig. 3a inset).

The susceptibility curves $\chi(\langle{k}\rangle)$ of the in-container networks for various $\theta$ also collapse onto that of the ring CP model (Fig. 3b). For the lifted networks, $\chi(\langle{k}\rangle)$ shifts to larger $\langle{k}\rangle$ as $\theta$ increases, but it always peaks near $k_{1/2}(\theta)$ (Extended Data Fig. 4f). This aligns with percolation theory, which predicts a susceptibility peak near the percolation threshold [23]. The peak height of $\chi$ for the lifted networks increases with $\theta$ at $\theta\leq 72^{\circ}$ (Fig. 3b inset), indicating larger non-giant clusters. This occurs because the emergence of the giant cluster, capable of absorbing the non-giant clusters, is delayed to higher $\langle{k}\rangle$ as $\theta$ increases.

The vertical length $L_{z}$ of a lifted cluster reflects the eccentricity of the top node, i.e., the maximum distance from this node to any other node [42]. Since the top particle is effectively a randomly chosen node, $L_{z}/D$ approximately measures the node-averaged eccentricity, as confirmed in Fig. 3c. Unlike the monotonically increasing $S_{1}$, $L_{z}$ of the largest cluster exhibits a peak near the percolation threshold in both simulations (Fig. 3d) and experiments (Fig. 1b, Extended Data Fig. 5a). When the $\langle{k}\rangle$ axis is normalized by $k_{1/2}$, the $L_{z}/D$ curves with different $\theta$ collapse onto the average eccentricity of the ring CP’s largest cluster (Fig. 3d), indicating that, given $k_{1/2}$, $L_{z}(\langle{k}\rangle)$ can be predicted by the ring CP model. Meanwhile, the width $W$ of the cluster monotonically increases and plateaus as $\langle{k}\rangle$ increases (Extended Data Fig. 5b). After prolonged shaking, the lifted cluster almost retains its in-container shape (Fig. 1c), which determines the saturated $L_{z}$ and $W$.

The clustering coefficient of a node $i$, $C_{i}$, quantifies the tendency of its neighbors to connect with each other, i.e., forming triangles [25]. The mean clustering coefficient $C$ is the average of $C_{i}$ over all nodes in the network. Interestingly, the $C(\langle{k}\rangle)$ curves for all in-container and lifted networks with different $\theta$ collapse onto a master curve that closely follows $C(\langle{k}\rangle)$ of the ring CP (Fig. 3e). This master curve is lower than $C(\langle{k}\rangle)$ of the sphere CP model but much higher than the vanishing $C(\langle{k}\rangle)=\langle{k}\rangle/N$ of ER random networks when $N\gg\langle{k}\rangle$ [39]. The similarity in local connectivity patterns between the C-particle networks and the ring CP model is also evident in the eigenvalue spectra of their adjacency matrices (Fig. 3f, Extended Data Fig. 9).

The cluster size distribution $n_{s}\equiv N_{s}/N$, where $N_{s}$ is the number of clusters of size $s$, is commonly used in percolation studies. According to percolation theory, at the threshold, $n_{s}\sim s^{-\tau}$, where $\tau$ is determined by the universality class [22]. Figure 3g shows that when $\langle k\rangle$ is near the threshold observed in Fig. 3a, b or 4b, both the C-particle networks and the ring CP model closely follow $n_{s}\sim s^{-2.19}$, suggesting that they fall into the standard 3D percolation universality class where $\tau=2.19$ [43]. Furthermore, $n_{s}$ of in-container clusters agrees with the ring CP model over various $\langle k\rangle$, not only at the threshold (Extended Data Fig. 6a–d). $n_{s}$ of lifted clusters also matches the ring CP model before a giant cluster forms, but deviates afterward as the giant cluster partially breaks during lifting (Extended Data Fig. 6e–h).

III.4 Finite-size scaling

Phase transitions are affected by the system size, i.e., the number of particles $N$. By scaling data from various $N$, finite-size scaling (FSS) can determine the critical point and critical exponents [44, 45]. In CP, the order parameter and susceptibility exhibit the following FSS forms near the percolation threshold $\eta_{\textrm{c}}$:

[TABLE]

Here, $\bar{\nu}\equiv d\nu$ where the space dimension $d=3$. $\beta$, $\gamma$, and $\nu$ are the critical exponents. The scaling functions $\tilde{S_{i}}$ and $\tilde{\chi}$ are independent of $N$. The reduced density $\eta\equiv nv_{\textrm{ex}}$, where $n$ is the number density, and $v_{\textrm{ex}}$ is the excluded volume [38, 39]. $v_{\textrm{ex}}$ is the average accessible volume of the center of a particle when it is connected to another fixed particle. We analytically derive $v_{\textrm{ex}}^{\textrm{ring}}=\pi D^{3}/3$ for rings with diameter $D$ (Methods, Extended Data Fig. 10). $\eta$ approaches $\langle{k}\rangle$ as $N$ increases (Fig. 4a) because they both reflect the connection probability.

We first determine $\eta_{\textrm{c}}$ of the ring CP model from MC simulations. According to equation (3), $S_{2}/S_{1}$ is independent of $N$ at $\eta_{\textrm{c}}$ [45]. This is confirmed by the unique intersection of $S_{2}(\eta)/S_{1}(\eta)$ curves for different $N$ in Fig. 4b, giving $\eta_{\textrm{c}}=2.11$. This suggests that the mean degree at the percolation threshold is $\langle{k}\rangle_{\textrm{c}}=2.11$ when $N\to\infty$, substantially larger than $\langle{k}\rangle_{\textrm{c}}=1$ of ER random networks [25] and comparable to other CP models: $\langle{k}\rangle_{\textrm{c}}=2.74$ for spheres [39] and $\langle k\rangle_{\textrm{c}}=2.27$ for disks in 3D space [40]. This large $\langle{k}\rangle_{\textrm{c}}$ can be attributed to the threshold-increasing effect of spatial constraints [46]. The threshold can also be estimated from $k_{1/2}$ (Fig. 3a), the peak of $\chi$ (Fig. 3b), or the peak of $S_{2}$ (Extended Data Fig. 7b,e), but the intersection in Fig. 4b provides the most precise value.

The ring CP model demonstrates the expected FSS behavior, conforming precisely to Eqs. (3) and (4). By using $\eta_{\textrm{c}}=2.11$ and the critical exponents for standard 3D percolation ($\beta=0.41$, $\gamma=1.80$, and $\bar{\nu}=2.64$), we obtain clear data collapses of $S_{1}$, $S_{2}$, and $\chi$ across all $N$ (Fig. 4c,d, Extended Data Fig. 7). This shows that the ring CP model belongs to the same universality class as lattice percolation and other CP models [23].

To check whether the C-particle networks also satisfy FSS, we attempt data collapses using DEM simulations for different $N$. Unlike CP models, $\langle{k}\rangle$ and cluster sizes in C-particle networks are not determined solely by $\eta$ (Extended Data Fig. 8). This is because granular particles have correlated positions and orientations due to mechanical interactions under gravity and container constraints.

Therefore, instead of $\eta$, we employ an effective reduced density, $\eta_{\textrm{eff}}(\langle{k}\rangle;N)$, defined as the value of $\eta$ in the ring CP model that yields the same $\langle k\rangle$ as the C-particle network. We compute $\eta_{\textrm{eff}}$ by applying the empirical fit for $\eta(\langle k\rangle;N)$ from the CP model (Fig. 4a), using the measured $\langle{k}\rangle$ and $N$ for the C-particles.

Using $\eta_{\textrm{eff}}$ with the same $\eta_{\textrm{c}}$ and critical exponents from Fig. 4c,d, FSS yields reasonably good data collapses for $S_{1}$, $S_{2}$, and $\chi$ (Fig. 4e,f). The robust collapse in the dilute regime ($\eta_{\textrm{eff}}<\eta_{\textrm{c}}$) indicates that the C-particle networks exhibit critical behavior similar to the standard percolation universality class. The imperfect collapse of $\chi$ in the dense regime ($\eta_{\textrm{eff}}>\eta_{\textrm{c}}$) is possibly due to strong interparticle correlations.

IV Discussion

Our experiments and simulations reveal that the entanglement networks of C-shaped granular particles under vibration exhibit a percolation transition. Striking similarities between the C-particle networks (both in-container and lifted) and our proposed ring CP model are observed across various properties, such as the order parameter and susceptibility of percolation, degree distribution, average eccentricity (vertical length), mean clustering coefficient, and adjacency spectrum. Their cluster size distribution and finite-size scaling follow the standard percolation universality class. These similarities demonstrate that continuum percolation theory is applicable to entangled granular structures.

We have also found that particle shape impacts mechanical bonding when pulled: as the opening angle increases, geometric links act less as mechanical bonds, leading to cluster disintegration and a higher percolation threshold in lifted networks. Additionally, the mean degree of the C-particle networks increases logarithmically with vibration time, exemplifying the widely observed phenomenon of logarithmic aging in disordered systems.

These findings highlight the potential of network-based approaches in studying entangled granular materials. Future studies may incorporate additional complexities, such as the particles’ relative positions, orientations, velocities, and interactions, by assigning features to nodes and edges in the network model. The spatial data can be measured experimentally using X-ray tomography [47]. Furthermore, graph-based machine learning models can be developed to predict network dynamics and thus the mechanical response of entangled materials [48]. This network framework can be extended to various other entangled materials, such as mechanically interlocked molecules [49, 50], polymer chains [51], kinetoplast DNA [20, 21], organisms with branching structures [14], and synthetic nonconvex particles from micro [52, 53, 54] to macro scales [8]. Such advancements will help predict and control their mechanical behavior for diverse applications, including molecular machines [49, 50], entangled robots [15], and granular metamaterials [17, 16, 8].

V Methods

V.1 Experiment methods

We used C-particles made of type 304 stainless steel, which are strong enough to lift thousands of particles with minimal deformation. To ensure low and uniform surface roughness, we ground freshly made C-particles in a vibrating container for 4 hours until their surfaces became shiny. This duration is sufficient because, after the first hour of shaking, the metal powder generated by grinding was considerably reduced.

To prepare consistent and well-disentangled initial states, the C-particles are trickled through a vibrating funnel with a hole of internal diameter $2.6~\rm{cm}$ (Extended Data Fig. 1a and Supplementary Video 1). The dripping particles are collected by a second funnel, whose bottom hole (diameter: $3.4~\rm{cm}$) is blocked by a tape. The second funnel remains stationary during this process. Once all dripping particles are collected, the second funnel is placed on the electrodynamic shaker (ECON EDS-300 model). Starting at $t=0$, a vertical displacement of $z(t)=A\sin(2\pi ft)$ is applied. We set $f=20~\rm{Hz}$ and $A=1.55~\rm{mm}$, and the corresponding peak acceleration is $A\cdot(2\pi f)^{2}=2.5g$, where $g=9.8~\rm{m/s^{2}}$ (Extended Data Fig. 1a and Fig. 1a).

After vibration, we use a hook made of a paper clip to gently hook and lift the C-particle clusters from the surface one by one (Fig. 1b,c and Supplementary Video 2). The sizes of the clusters are measured by their weights. We neglect very small clusters ($<3.5$ grams) because their sizes are sensitive to the way they are hooked up and do not significantly affect the $S_{1}$ and $\chi$ behaviors. After measurement, all clusters are poured back into the first vibrating funnel (Extended Data Fig. 1a) to be disentangled in preparation for the next trial. To obtain sufficient statistics, we repeat this procedure multiple times for each vibration $t$, as listed in Supplementary Table S3-S5.

V.2 Discrete element method (DEM) simulation

In DEM simulations, the geometries of the funnel container and C-particles are similar to those in our experiment. We use 10 different sets of rigid C-particles with opening angles of $20^{\circ}$, $33^{\circ}$, $46^{\circ}$, $59^{\circ}$, $72^{\circ}$, $98^{\circ}$, $112^{\circ}$, $125^{\circ}$, $138^{\circ}$, and $151^{\circ}$ which are composed of 26, 25, 24, 23, 22, 20, 19, 18, 17, and 16 identical small hard spheres, respectively (Supplementary Fig. S2, Table S1). The diameter $D$ of the C-particles is 8.8 times the diameter $d$ of the small spheres. This ratio is similar to the C-particles with $D=9$ mm and $d=1$ mm in the experiment. The total number of particles $N=4000$ is the same as in the experiment. For particles with $\theta=20^{\circ}$, $N=500$ and 1400 are also simulated for the finite-size scaling. The friction coefficient is set to 0.4 for both particle-particle and particle-container interfaces. Details on the DEM contact model are provided in Supplementary Information.

To prepare a disentangled initial state, C-particles are randomly created with a volume fraction of 0.05 in a vertical cylindrical region with diameter $9D$ and height $22D$ right above the container. The particles then fall into the container under gravity (Supplementary Video 3).

In the entanglement process, the particles are shaken with an amplitude of $A=0.172D$ (Supplementary Video 4) and a peak acceleration of $2.5g$ to replicate the experiment.

After being vibrated for a desirable time, individual clusters near the surface are lifted one by one (Supplementary Video 5-7). The lifting speed $0.5\sqrt{gD}$ is equivalent to the speed of a particle that has fallen a relatively short distance, $D/8$. This speed is chosen to gently lift the cluster while minimizing the overall simulation time. If the bottom of the lifted cluster is more than $1.1d$ above the top particle remaining in the container, we consider the cluster fully lifted and remove it. This process is repeated until all particles have been removed. The number of trials are listed in Supplementary Table S6,S7.

V.3 Monte Carlo (MC) simulation

The continuum percolation (CP) model consisting of closed, infinitely thin rings is simulated using Monte Carlo methods. The ring centers, $\bm{r}_{i}$ for $i=1,2,...,N$, are randomly generated within a cubic volume $V$ with non-periodic boundaries. The normal orientation of the $i$’th ring is described by the unit vector $\hat{\bm{n}}_{i}=\left[\sin(\theta_{i})\cos(\phi_{i}),\,\sin(\theta_{i})\sin(\phi_{i}),\,\cos(\theta_{i})\right]$ where the two angles in the spherical coordinate $\theta_{i}=\cos^{-1}(2x_{i}-1)$ and $\phi_{i}=2\pi x_{i}$. Here, $x_{i}$ is a random variable with a uniform distribution in $[0,1]$. Snapshots at different reduced densities $\eta$ are shown in Supplementary Fig. S1. The ensemble averages of the relevant variables are obtained at each $\eta$. The numbers of trials are listed in Supplementary Table S8,S9.

In contrast to our newly proposed ring CP model, the sphere CP model has already been intensively studied as the simplest CP model in 3D [38, 39]. In each MC simulation of sphere CP model, $N=4000$ monodispersed spheres are randomly distributed within a cube with non-periodic boundaries. Two spheres are considered connected if they overlap, i.e., their distance is smaller than their diameter. The number of trials used in Fig. 3e is listed in Supplementary Table S8. The eigenvalue spectra (Extended Data Fig. 9g-i) are averaged over 300 trials.

In addition to ring and sphere CP models, we also simulate ER random networks to compare their spectra with the C-particle network. To achieve a desired mean degree of $\langle{k}\rangle$, $N\langle{k}\rangle/2$ connections are randomly attached to $N=4000$ nodes with equal probability. The eigenvalue spectra of the networks (Fig. 3f and Extended Data Fig. 9j-l) are averaged over 300 trials.

V.4 Entanglement criterion

Entanglement generally refers to geometric constraints between elongated or non-convex bodies that cause mechanical coupling. It can be quantified in various ways. For example, in filamentous materials such as polymers and rods, entanglement can be statistically estimated through tube models [51], or measured by the average crossing number over all directions [24, 55].

For C-particles, we adopt a topological entanglement criterion from Ref. [19]: two particles are defined as entangled when the circles passing through their centerlines form a Hopf link (Extended Data Fig. 10). This definition is implemented in our DEM and MC simulations. Mathematically, two circles $i$ and $j$ form a Hopf link if and only if

[TABLE]

where $R$ is their radius, $\bm{r}_{i}$ and $\bm{r}_{j}$ are the positions of their centers, and $\hat{\bm{n}}_{i}$ and $\hat{\bm{n}}_{j}$ are the unit vectors normal to their planes. $\bm{k}\equiv\hat{\bm{n}}_{i}\times\hat{\bm{n}}_{j}$, $\bm{r}_{ij}\equiv\bm{r}_{j}-\bm{r}_{i}$, and $\bm{b}\equiv({\hat{\bm{n}}_{j}\cdot\bm{r}_{ij}}/{|\bm{k}|^{2}})\left[\left(\hat{\bm{n}}_{i}\cdot\hat{\bm{n}}_{j}\right)\hat{\bm{n}}_{i}-\hat{\bm{n}}_{j}\right]$.

V.5 Mean clustering coefficient

In Fig. 3e, the mean clustering coefficient $C\equiv\sum_{i=1}^{N}{C_{i}}/N$, where $N$ is the total number of nodes in the network, and $C_{i}$ is the local clustering coefficient of node $i$ [25]. $C_{i}\equiv 2t_{i}/(k_{i}^{2}-k_{i})$, where $t_{i}$ is the number of triangles (loops of length 3) attached to node $i$, and $k_{i}$ is the degree of node $i$. $C_{i}=0$ when $k_{i}$ is 0 (isolated node) or 1 (leaf node).

V.6 Calculation of the excluded volume of rings

In CP models, the connection probability of two objects is determined by the excluded volume $v_{\rm{ex}}$: when edge effects are negligible, $p=v_{\rm{ex}}/V$, where $V$ is the total volume of the system. $v_{\rm{ex}}$ is generally defined as the volume around an impenetrable object within which the center of another identical object cannot access because of the presence of the first object. For example, $v_{\rm{ex}}=4\pi D^{3}/3$ for spheres with diameter $D$ [39]. In the CP model of infinitely thin rings, $v_{\rm{ex}}$ is equivalent to the volume accessible to the center of a ring when it is connected, i.e., forming a Hopf link, to another fixed ring.

Inspired by Onsager’s calculation of the excluded volume of cylinders [56, 57], we analytically derive the excluded volume of infinitely thin rings, $v_{\rm{ex}}^{\rm{ring}}$, as follows: Consider two connected rings of radius $R$ whose centers are $C_{1}$ and $C_{2}$, as shown in Extended Data Fig. 10a. To form a Hopf link, the blue ring should enclose either $P$ or $Q$, but not both. This can be conveniently seen in the $X$-$Y^{\prime}$ plane as shown in Extended Data Fig. 10b: For given $y$ and $\gamma$, $C_{2}$ can only be located within the shaded area, which is formed by two circles with radius $R$ centered at $P$ and $Q$. The shaded area can be found by subtracting the overlap area from the total area of two circles: $A_{\rm{shaded}}(y)=2\left[\pi R^{2}-A_{\rm{overlap}}(y)\right]$. The overlap area $A_{\rm{overlap}}(y)=2\left[R^{2}\theta(y)-x(y)h(y)\right]$, where $\theta(y)=\cos^{-1}\left({x(y)}/{R}\right)$, $x(y)=\sqrt{R^{2}-y^{2}}$, and $h(y)=y$. Therefore, the excluded volume for a fixed $\gamma$ can be calculated as

[TABLE]

Given that the rings’ orientations are random, we take the average over all orientations to obtain

[TABLE]

where $D=2R$. We further verified this result using MC simulation.

The excluded volume of disks (circular plates) of diameter $D$ is $v_{\textrm{ex}}^{\textrm{disk}}=\pi^{2}D^{3}/8$ [56], slightly larger than $v_{\rm{ex}}^{\rm{ring}}$, as expected. Thus the reduced density (equal to the mean degree as $N\to\infty$) of disks is $\eta={\pi^{2}}nD^{3}/8$. Using this and the critical number density $\pi n_{\textrm{c}}D^{3}/6=0.9614$ from Ref. [40], we obtain $\langle{k}\rangle_{\textrm{c}}\equiv\eta_{\textrm{c}}=2.27$ for disks.