Percolation transition in entangled granular networks

Domain: physics
Published: 2025-12-10
Editorial outcome: minor_revision

Abstract
Highly nonconvex granular particles, such as staples and metal shavings, can form solid-like cohesive structures through geometric entanglement (interlocking). However, the network structure formed by this entanglement remains largely unexplored. Here, we employ network science to investigate the entanglement networks of C-shaped granular particles under vibration in experiments and simulations. Analysis of key network properties reveals 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 this observed transition. Furthermore, we find that the particles’ opening angle significantly affects mechanical bonding and, consequently, the network structure. This work demonstrates the promise of network-based approaches for studying entangled materials, with potential applications from mechanical metamaterials to entangled robot swarms.

Introduction
Particles shape plays a crucial role in the collective behavior of granular materials. Unlike the well-studied spherical particles1,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 cohesion3–5. For example, entangled U-shaped4, Z-shaped6, star-shaped7–9, or cross-shaped10 particles can resist uniaxial compression and maintain freestanding columns. Moreover, S- or U-shaped particles can withstand tensile stresses11,12 and be collectively lifted against gravity8,13. These mechanical effects arising from entanglement have recently sparked interest across diverse fields, such as architecture8, living systems14,15, soft robotics16, and metamaterials17,18.

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 grains19, 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 colloids20 and kinetoplast DNA21,22. 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 network analysis.

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.

Results
We experimentally measure the cluster formation of steel C-particles under controlled vibration. Each particle has a toroidal shape with opening angle θ, centerline diameter D = 9 mm, and thickness d = 1 mm (Fig. 1a). We test nine sets of C-particles with θ ranging from 25° to 135°, 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 θ 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 θ ranging from 20° to 151°. 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, Supplementary Fig. 10)20. Although C-particles are not closed loops, this criterion well approximates their mechanical interlocking, especially when θ 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.

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, S1 ≡ s1/N, serves as the order parameter, where s1 represents the number of particles in the ith largest cluster, and N is the total number of particles23. The experiments (Fig. 1e) and simulations (Fig. 1f, g) show that for small θ, lifted S1(t) steadily grows and approaches 1, indicating percolation of entanglement P. Additionally, the susceptibility, defined as the mean cluster size χ ≡ Σi≠1 si2/N, peaks near the rapid growth of S1(t), i.e., the percolation threshold (critical point), as predicted by percolation theory (Fig. 1h, i)24.

For θ ≲ 70°, C-particles with θ show steady growth of S1 at t > 0 (Fig. 1e, f). However, for θ ≳ 70°, this trend reverses (Fig. 1g, j) because larger θ makes the links break more easily. When only a few tiny clusters detach from the giant cluster, S1 remains close to 1, but when the giant cluster disintegrates during lifting, S1 decreases dramatically. These two cases result in a bimodal S1 distribution at large θ (Supplementary Fig. 1c, e).

The clusters become extremely fragile beyond a certain θ. In experiments at t = 60 s, giant clusters can be lifted stably for θ ≤ 110°, but often disintegrate for θ ≥ 115°. As a result, the ensemble-averaged S1 at t = 60 s plummets around θ = 115° (Supplementary Fig. 1f). This mechanical instability, emerging beyond a certain θ, is also confirmed by simulations. When θ ≳ 150°, the clusters are too fragile to be lifted as giant clusters (Fig. 1g, j, Supplementary Fig. 1g).

This maximum θ for stability depends on particle thickness and surface friction. Ref. 25 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. Reducing particle thickness or increasing the friction coefficient would enhance cohesion, there by raising the maximum stable θ and shrinking the fragile region in Fig. 1j.

Evolution of the degree distribution
To investigate the entangled structure, we model the granular system as an N-node network comprising multiple disconnected clusters, 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 connectivity26. The observed P(k) of the C-particle networks at various t mostly agrees with the Poisson distribution (Fig. 2b, c):

P(k) = e−〈k〉 〈k〉k / k!, (1)

where 〈k〉 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 θ (Fig. 2c), except for lifted networks with large θ that undergo frequent cluster breakages (Supplementary Fig. 2).

The Poisson P(k) is a feature of both Erdös-Rényi (ER) random networks26 and continuum percolation (CP) models, also known as random geometric graphs27. A similar P(k) is observed in kinetoplast DNA networks, which are quasi-2D assemblies of flexible loops21,22. These results indicate that an underlying spatial randomness is shared by these systems, despite differences in geometry and dynamics. 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 〈k〉, knowing how 〈k〉 changes with θ and t is crucial. Figure 2d shows that 〈k〉(t) of the in-container networks can be fitted by

〈k〉(t) = k0 + α ln(1 + t/t0), (2)

where the fitting parameters k0 and t0 vary with θ, while α is relatively independent of θ (Supplementary Fig. 3a–c). This logarithmic increase persists over a wider range of t, even long after S1(t) (Fig. 1f) has nearly saturated. The lifted networks also display logarithmic 〈k〉(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 grains28. However, in our system, the mean degree exhibits a clearer logarithmic trend, whereas the density shows significantly larger fluctuations (Supplementary 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 polymers29 and granular media30, structural relaxation in DNA31 and crumpled paper32, as well as magnetization decay in spin glasses33. In these systems, slow relaxation typically arises from sequential transitions between metastable states separated by broadly distributed waiting times, often due to rugged energy landscapes34–36. These systems also exhibit history dependence memory effects37. 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.

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 have been utilized to study various disordered systems, including porous media, semiconductors, and wireless networks38,39. In conventional CP models, objects such as spheres40, disks41, or rods42 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 (Figs. 3, 4).

Although C-particles are neither infinitely thin nor completely randomly located, their entanglement networks are well captured by the ring CP model (Figs. 3, 4). The in-container S1(〈k〉) curves for various θ collapse onto the ring CP (Fig. 3a), but the lifted networks deviate from the ring CP (Fig. 3a), i.e., by comparing 〈k〉, they are lower than the ring CP networks with the same S1. This deviation is stronger at larger θ, because large openings facilitate the disintegration of the giant cluster, while the decrease in 〈k〉 remains limited owing to the tightly connected core. The percolation threshold can be roughly estimated at S1(k1/2) = 1/2, and the measured k1/2 increases linearly with θ (Fig. 3a inset).

The susceptibility curves χ(〈k〉) of the in-container networks for various θ also collapse onto that of the ring CP model (Fig. 3b). For the lifted networks, χ(〈k〉) shifts to larger 〈k〉 as θ increases, but it always peaks near k1/2(θ) (Supplementary Fig. 4f). This aligns with percolation theory, which predicts a susceptibility peak near the percolation threshold24. The peak height of χ for the lifted networks increases with θ at θ ≤ 72° (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 〈k〉 as θ increases.

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

The clustering coefficient of a node i, Ci, quantifies the tendency of its neighbors to connect with each other, i.e., forming triangles26. The mean clustering coefficient C is the average of Ci over all nodes in the network. Interestingly, the C(〈k〉) curves for all in-container and lifted networks with different θ collapse onto a master curve that closely follows C(〈k〉) of the ring CP (Fig. 3e). This master curve is lower than C(〈k〉) of the sphere CP model but much higher than the vanishing C(〈k〉) = 〈k〉/N of ER random networks when N ≫ 〈k〉40. 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, Supplementary Fig. 9).

The cluster size distribution ns ≡ Ns/N, where Ns is the number of clusters of size s, is commonly used in percolation studies. According to percolation theory, at the threshold, ns ~ s−τ, where τ is determined by the universality class23. Figure 3g shows that when 〈k〉 is near the threshold observed in Figs. 3a, b or 4b, both the C-particle networks and the ring CP model closely follow ns ~ s−2.19, suggesting that they fall into the standard 3D percolation universality class where τ = 2.1944. Furthermore, ns of in-container clusters agrees with the ring CP model over various 〈k〉, not only at the threshold (Supplementary Fig. 6a–d). ns 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 (Supplementary Fig. 6e–h).

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 exponents45,46. In CP, the order parameter and susceptibility exhibit the following FSS forms near the percolation threshold ηc:

S1 Nβ/ν = S~((η − ηc)N1/ν), (3)

χ N−γ/ν = χ~((η − ηc)N1/ν). (4)

Here, ν ≡ dν where the space dimension d = 3. β, γ, and ν are the critical exponents. The scaling functions S~ and χ~ are independent of N. The reduced density η ≡ nvex, where n is the number density, and vex is the excluded volume39,40. vex is the average accessible volume of the center of a particle when it is connected to another fixed particle. We analytically derive vexring = πD3/3 for rings with diameter D (“Methods”, Supplementary Fig. 10). η approaches 〈k〉 as N increases (Fig. 4a) because they both reflect the connection probability.

We first determine ηc of the ring CP model from MC simulations. According to Eq. (3), S1/S2 is independent of N at ηc46. This is confirmed by the unique intersection of S1(η)/S2(η) curves for different N in Fig. 4b, giving ηc = 2.11. This suggests that the mean degree at the percolation threshold is 〈k〉c = 2.11 when N → ∞, substantially larger than 〈k〉c = 1 of ER random networks26 and comparable to other CP models: 〈k〉c = 2.74 for spheres40 and 〈k〉c = 2.27 for disks in 3D space41. This large 〈k〉c can be attributed to the threshold-increasing effect of spatial constraints47. The threshold can also be estimated from k1/2 (Fig. 3a), the peak of χ (Fig. 3b), or the peak of S2 (Supplementary 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 ηc = 2.11 and the critical exponents for standard 3D percolation (β = 0.41, γ = 1.80, and ν ≡ 2.64), we obtain clear data collapses of S1, S2, and χ across all N (Fig. 4c, d, Supplementary Fig. 7). This shows that the ring CP model belongs to the same universality class as lattice percolation and other CP models24.

To check whether the C-particle networks also satisfy FSS, we attempt data collapses using DEM simulations for different N. Unlike CP models, 〈k〉 and cluster sizes in C-particle networks are not determined solely by η (Supplementary Fig. 8). This is because granular particles have correlated positions and orientations due to mechanical interactions under gravity and container constraints. Therefore, instead of η, we employ an effective reduced density, ηeff(〈k〉; N), defined as the value of η in the ring CP model that yields the same 〈k〉 as the C-particle network. We compute ηeff by applying the empirical fit for η(〈k〉; N) from the CP model (Fig. 4a), using the measured 〈k〉 and N for the C-particles. Using ηeff with the same ηc and critical exponents from Fig. 4c, d, FSS yields reasonably good data collapses for S1, S2, and χ (Fig. 4e, f). The robust collapse in the dilute regime (ηeff < ηc) indicates that the C-particle networks exhibit critical behavior similar to the standard percolation universality class. The imperfect collapse of χ in the dense regime (ηeff > ηc) is possibly due to strong interparticle correlations.

Methods
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 cm (Supplementary Fig. 1a and Supplementary Movie 1). The dripping particles are collected by a second funnel, whose bottom hole (diameter: 3.4 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πft) is applied. We set f = 20 Hz and A = 1.55 mm, and the corresponding peak acceleration is A ⋅ (2πf)2 = 2.5g, where g = 9.8 m/s2 (Supplementary 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 Movie 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 S1 and χ behaviors. After measurement, all clusters are poured back into the first vibrating funnel (Supplementary 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 Tables 3–5.

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°, 33°, 46°, 59°, 72°, 98°, 112°, 125°, 138°, and 151° which are composed of 26, 25, 24, 23, 22, 20, 19, 18, 17, and 16 identical small hard spheres, respectively (Supplementary Fig. 12, Table 1). 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 θ = 20°, 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 Movie 3). In the entanglement process, the particles are shaken with an amplitude of A = 0.172D (Supplementary Movie 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 Movie 5–7). The lifting speed 0.5 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 6, 7.

Monte Carlo (MC) simulation
The continuum percolation (CP) model consisting of closed, infinitely thin rings is simulated using Monte Carlo methods. The ring centers, ri 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 n^i = (sin(θi)cos(ϕi), sin(θi)sin(ϕi), cos(θi)) where the two angles in the spherical coordinate θi = cos−1(2xi − 1) and ϕi = 2πxi. Here, xi is a random variable with a uniform distribution in [0, 1]. Snapshots at different reduced densities η are shown in Supplementary Fig. 11. The ensemble averages of the relevant variables are obtained at each η. The numbers of trials are listed in Supplementary Table 8, 9.

In contrast to our newly proposed ring CP model, the sphere CP model has already been intensively studied as the simplest CP model in 3D39,40. 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 8. The eigenvalue spectra (Supplementary 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 〈k〉, N〈k〉/2 connections are randomly attached to N = 4000 nodes with equal probability. The eigenvalue spectra of the networks (Fig. 3f and Supplementary Fig. 9j–l) are averaged over 300 trials.

Entanglement criterion
Entanglement generally refers to geometric constraints between elongated or non-convex bodies that cause mechanical bonding. It can be quantified in various ways. For example, in filamentous materials such as polymers and rods, entanglement can be statistically estimated through tube models52, or measured by the average crossing number over all directions25,56.

For C-particles, we adopt a topological entanglement criterion from Ref. 20: two particles are defined as entangled when the circles passing through their centerlines form a Hopf link (Supplementary 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

|R2 − |rj − ri|2|2 > |k|2 |rj − ri|2 + b ⋅ (rj − ri),

where R is their radius, ri and rj are the positions of their centers, and n^i and n^j are the unit vectors normal to their planes. k ≡ n^i × n^j, rij ≡ rj − ri and b ≡ ((n^i ⋅ rij / |k|2) n^j − (n^j ⋅ rij / |k|2) n^i).

Mean clustering coefficient
In Fig. 3e, the mean clustering coefficient C ≡ Σi=1N Ci/N, where N is the total number of nodes in the network, and Ci is the local clustering coefficient of node i26. Ci ≡ 2ti / (ki(ki − 1)), where ti is the number of triangles (loops of length 3) attached to node i, and ki is the degree of node i. Ci = 0 when ki is 0 (isolated node) or 1 (leaf node).

Calculation of the excluded volume of rings
In CP models, the connection probability of two objects is determined by the excluded volume vex: when edge effects are negligible, p = vex/V, where V is the total volume of the system. vex 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, vex = 4πD3/3 for spheres with diameter D40. In the CP model of infinitely thin rings, vex 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 cylinders57,58, we analytically derive the excluded volume of infinitely thin rings, vexring, as follows: Consider two connected rings of radius R whose centers are C1 and C2, as shown in Supplementary 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-Y0 plane as shown in Supplementary Fig. 10b: For given y and γ, C2 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: Ashaded(y) = 2πR2 − Aoverlap(y). The overlap area Aoverlap(y) = 2R2θ(y) − x(y)h(y), where θ(y) = cos−1(x(y)/R), x(y) = R2 − y2, and h(y) = y. Therefore, the excluded volume for a fixed γ can be calculated as

vexring(γ) = 2 ∫0R Ashaded(y) sinγ dy = 32/3 R3 sinγ.

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

vexring = 1/4π ∫02π ∫0π vexring(γ) sinγ dγ dϕ = π2D3/3,

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

The excluded volume of disks (circular plates) of diameter D is vexdisk = π2D3/857, slightly larger than vexring, as expected. Thus the reduced density (equal to the mean degree as N → ∞) of disks is η = π2nD3/8. Using this and the critical number density πnD3/6 = 0.9614 from Ref. 41, we obtain 〈k〉c ≡ ηc = 2.27 for disks.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is a manuscript on the percolation of entangled granular networks which has a certain charme since it is easy to understand and the corresponding experiment is macroscopic. Still the physics behind the entanglements involving a percolation transition is interesting and fascinating. After having read the paper, I think that the authors did a careful job in posing the problem and performing the experiment. It is always good to have a computer simulation in parallel. The interdisciplinary topic (granules, statistical physics, clusters etc) is ideal for Nature Communications. But before publication the following points need to be addressed by the authors:\n1) The depiction of a node network introduced in line 146 is much too brief. In the main text, in Figure 2a, it is sketched but should be discussed in more detail. What are the dots surrounding the cluster here and what is their distance? What is the physical significance of this representation?\n2) I worry about the protocol dependence of the results. The duration of shaking is one important parameter in preparing the system. It is mentioned that 578 experiments were performed but under which conditions? For example: Has a shaken system been lifted , put on the cone again and then shaken again in order to make sure that the result does not depend on the preparation and initial history of the sample? And have initially completely disentangled C-particles been thrown into the cone and did then the shaking was yield the same results as for an initially entangled configuration? Since the time to produce enanglements can be huge this should be carefully checked in order to establish the full protocol independence of the data. Sure the same should be done in the MC simulation but I would in particular see this confirmed in the experiments.", "reviewer_3": "The authors combine experiment, discrete element (DEM) simulations and Monte Carlo simulations of a continuum model to study the entanglement (linking) properties of systems of C-shaped particles with different opening angles theta. Both in experiments and in the DEM simulations, entanglement between the particles is achieved by subjecting the system to sinusoidal vibration.\nThe authors analyze the resulting entanglement network (in which each node represents a particle and each edge a link/entanglement) by measuring several graph theoretical properties, such as the size of the largest cluster, the average degree (connectivity) of the nodes, and the cluster size distribution. The authors additionally compare their results to known theoretical predictions for the 3D percolation model.\nThe main results can be summarized as follows:\n- Clusters of particles with larger opening angle theta entangle more easily than those of particles with small theta, but there is a critical theta (~70 degrees) after which this trend reverses, and the clusters entangle less and break more easily. The opening angle thus affects the topological and mechanical properties of the network.\n- Most properties of the entanglement network (such as degree distribution and cluster size distribution) are well described by the 3D continuum percolation model. The mean degree of the nodes increases logarithmically with vibration time.\nThe experimental methods employed seem sound, as well as the simulations and the analyses performed. I find the results and the methods clearly explained. However, I have to say I do not find the main results as presented to be very surprising (from the Abstract: “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.” [...] “Additionally, we find that particle’s opening angle significantly affects mechanical bonding and, consequently, the network structure”). The logarithmic increase of the number of links over time is interesting, but it is not motivated theoretically and the connection to other systems is a bit vague, as detailed below. I thus believe that the paper would benefit from a more detailed comparison with existing experimental and simulation results for similar systems.\nIn particular, a very similar work, studying rod-like particles instead of C-shaped particles via experiments and simulations, was recently published [1]. It would be interesting to compare the results of the present work with those of [1]. In particular, the authors of [1] also find a percolation of the entanglement network which results into mechanical stabilization, despite their particles being straight rods, which I believe for the C-shaped particles would roughly correspond to a particle with very large diameter/thickness (D/d) and opening angle approaching 360 degrees. This seems to partially contradict the observation that mechanical stability becomes compromised for opening angles larger than 150 degrees (pag. 2: “When θ ≳ 150, the clusters are too fragile to be lifted as giant clusters”), so this is a point that I would like to see addressed.\nAn even more similar system to the one studied by the authors, consisting of U-shaped particles, has been studied experimentally and in simulations by Gravish et al. [2] (Ref. 3 in the present manuscript). How do the results of this study differ from those of [2]? Is there any fundamental qualitative difference between the two systems and their behaviors? If so, these differences should be clearly explained, because the two systems look very similar to me, which makes me question the novelty of this study.\nFinally, since the authors mention kinetoplast DNA, I think it would be relevant to compare more in detail the results of this paper with those obtained from simulations of kinetoplast. I’m aware that the latter is a quasi-2D system whereas the system studied by the authors is 3D, but a comparison might still be informative.\nPlease find below some comments and questions that I believe should be addressed (in addition to those reported above) before this paper can be considered for publication.\n- I believe that the operative definition of entanglement used in the paper (two particles are considered entangled when the two circles passing through their centerlines form a Hopf link) should be reported in the main text for clarity (currently, the main text reports that “Two particles are defined as linked when the circles along their centerlines are intertwined”, but “intertwined” is not a well-defined concept).\n- Related to my previous point, I think that it might be useful to spend a few words on the meaning of entanglement. Due to the particular shape of the particles studied in this work, their entanglement can be straightforwardly defined using the topological concept of linking. However, this is not the case in other systems, in particular those made of slender objects such as polymers, living worms or rod-like particles (see [1]). I feel that the term “entanglement” tends to get used in a loose way to indicate rather different physical situations, and thus some clarification might help the readers who are not familiar with all the different definitions.\n- Is there any theoretical basis for the functional form reported in Eq. (2), or is this a purely empirical fit? Below this equation, the authors state that “similar logarithmic time evolution has been observed in various disordered systems”, reporting several references. However, from what I could see after a cursory look at these references, none of these study the time evolution of an entanglement network. What is, then, the justification to compare this logarithmic behavior to those observed in these papers, that concern very different quantities? I find the connection a bit vague and not justified well enough. I would invite the authors to compare their results to the time evolution of the entanglement network reported in [1] and [2], for example. What comes closest to a microscopic theoretical explanation seems to be Ref. 32, although I am not sure of how to compare its predictions with the results here presented.\n- The authors find using finite-size scaling and other methods that the mean node degree at percolation is <k_c>=2.11 for infinite system size. I find this results to make intuitive sense, as the network cannot percolate if the degree of each node is 2 or less, and a certain number of nodes of degree 3 or more is needed. Have the authors compared this predictions to analytical models for similar systems? For example, it is known that in systems of particles with a limited number of interaction sites (limited valence or “patchy” particles), the phase diagram of the system develops a phase-separated region as soon as the mean particle valence is above 2 – see for example [3]. This result can be explained analytically by the Flory-Stockmeyer combinatorial theory [4]. I wonder if there is any connection to draw between this result and the one obtained here.\n- I could not understand how the effective reduced density eta_{eff} (pag. 6) is defined/calculated. The definition of this quantity and how it is calculated should be made clearer in the main text.\n[1] Jung et al., Entanglement transition in random rod packings, PNAS 2025 Vol. 122 No. 8\n[2] Gravish et al., Entangled Granular Media, PRL 2012 108, 208001\n[3] Sciortino and Zaccarelli, Reversible gels of patchy particles, Current Opinion in Solid State and Materials Science 15 (2011) 246–253\n[4] Stockmayer W. Theory of molecular size distribution and gel formation in branched-chain polymers. J Chem Phys 1943, 11:45."}, "author_response": "Responses to Reviewers\nManuscript Title: Percolation Transition in Entangled Granular Networks\nManuscript ID: NCOMMS-25-28470\nSeongmin Kim, Daihui Wu, and Yilong Han*\n*yilong@ust.hk\nSeptember 4, 2025\nWe appreciate the opportunity to revise our manuscript and thank the referees for their con- structivefeedback. Below, weaddresstheirconcernspointbypoint. Changesinthemanuscriptare highlightededinblue. Referencenumbersinbrackets[]correspondtothoseintherevisedmanuscript.\nReviewer #1\nReviewer: This is a manuscript on the percolation of entangled granular networks which has a certain charm since it is easy to understand and the corresponding experiment is macroscopic. Still the physics behind the entanglements involving a percolation transition is interesting and fascinating. After having read the paper, I think that the authors did a careful job in posing the problem and performing the experiment. It is always good to have a computer simulation in par- allel. The interdisciplinary topic (granules, statistical physics, clusters etc) is ideal for Nature Communications. But before publication the following points need to be addressed by the au- thors:\nReviewer: The depiction of a node network introduced in line 146 is much too brief. In the main text, in Figure 2a, it is sketched but should be discussed in more detail. What are the dots surrounding the cluster here and what is their distance? What is the physical significance of this representation?\nResponse: Thank you for the comment. We agree that additional clarification for Fig. 2a would strengthen the discussion. The main purpose of Fig. 2a is to visualize the entanglement network. Also, the figure shows that the whole network contains multiple clusters (i.e., connected compo- nents). Moreover, the relatively homogeneous degrees represented by the color and size of the nodes show that this network is clearly different from the scale-free networks, which features the power-law degree distribution with a few extremely high-degree hubs.\nThe node placement in Fig. 2a is schematic (using the ‘force-directed layout (Rruchterman- Reingold layout)’ in MATLAB, which applies attractive forces between adjacent nodes and repulsive forces between distant nodes) and does not reflect geographic positions or any other physical quantities. We have clarified this in the caption of Fig. 2a.\nTo clarify these, we have added\n• Lines 168-169 (Fig. 2a caption):\nNode positions are set by a force-directed layout for visual clarity and do not re- flect actual spatial coordinates. Isolated dots represent disconnected single C- particles.\n• Lines 175–184:\nTo investigate the structural basis of the observed percolation transitions, we model the granular system as an N-node network comprising multiple discon- nected clusters, where nodes represent particles and edges denote topological links. Figure 2a shows a snapshot of the entanglement network from a DEM sim- ulation; the network exhibits relatively homogeneous node degrees (number of links, k) with no prominent hubs.\nThis degree homogeneity is evident in the degree distribution P(k), a fundamen- tal metric of network connectivity [25].\nReviewer: I worry about the protocol dependence of the results. The duration of shaking is one important parameter in preparing the system. It is mentioned that 578 experiments were performed but under which conditions? For example: Has a shaken system been lifted, put on the cone again and then shaken again in order to make sure that the result does not depend on the preparation and initial history of the sample? And have initially completely disentangled C-particles been thrown into the cone and did then the shaking was yield the same results as for an initially entangled configuration? Since the time to produce entanglements can be huge this should be carefully checked in order to establish the full protocol independence of the data. Sure the same should be done in the MC simulation but I would in particular see this confirmed in the experiments.\nResponse: Thank you for the comment. We agree that establishing full protocol independence is crucial in our experiment because the entanglement network exhibits strong history dependence (memory effects).\nSince our measurements are destructive, we reset the sample (i.e., erase its memory) after each trial using the disentanglement process, where entangled clusters are trickled through a vibrat- ing funnel to become disentangled (Extended Data Fig. 1a). Each of the 578 trials begins with a fresh sample for the vibration at t = 0, avoiding reuse of previously shaken and measured clusters. Thus, each data point in Extended Data Fig. 1b,c corresponds to a full cy- cle of disentangling → shaking → lifting all clusters → measuring → disentangling. The same protocol applies to DEM simulations, except the disentanglement process is unnecessary since disentangled particles can be generated directly. The exact number of trials for each condition is provided in Supplementary Tables S3–S7.\nTo clarify this, we have added\n2\n• Lines 87–89:\nBecause the measurement is destructive, each trial begins with a freshly disentan- gled sample.\n• Line 700–702 (‘Experiment methods’ section in Methods):\nAfter measurement, all clusters are poured back into the first vibrating funnel (Extended Data Fig. 1a) to be disentangled in preparation for the next trial.\nReviewer #3\nReviewer: The authors combine experiment, discrete element (DEM) simulations and Monte Carlo simulations of a continuum model to study the entanglement (linking) properties of sys- tems of C-shaped particles with different opening angles theta. Both in experiments and in the DEM simulations, entanglement between the particles is achieved by subjecting the system to sinusoidal vibration.\nThe authors analyze the resulting entanglement network (in which each node represents a parti- cle and each edge a link/entanglement) by measuring several graph theoretical properties, such as the size of the largest cluster, the average degree (connectivity) of the nodes, and the cluster size distribution. The authors additionally compare their results to known theoretical predictions for the 3D percolation model.\nResponse: We appreciate the reviewer’s thoughtful summary of our key findings. We would like to add that (1) we also analyzed the adjacency-matrix eigenvalue spectra, an approach that, to our knowledge, has not previously been applied to geometrically entangled systems; and (2) while continuum percolation models of overlapping objects (e.g., spheres, disks) are well established, our work is the first to propose and simulate continuum percolation of entangled (topologically linked) objects, where connections are defined by topological links. We derived, for the first time, a closed-form expression for the excluded volume of infinitely thin rings: vring = πD3/3, and ex identified the critical point η = 2.11—both of which are mathematically meaningful results.\nc\nThe main results can be summarized as follows:\n- Clusters of particles with larger opening angle theta entangle more easily than those of particles with small theta, but there is a critical theta ( 70 degrees) after which this trend reverses, and the clusters entangle less and break more easily. The opening angle thus affects the topological and mechanical properties of the network.\n- Most properties of the entanglement network (such as degree distribution and cluster size dis- tribution) are well described by the 3D continuum percolation model. The mean degree of the nodes increases logarithmically with vibration time.\nThe experimental methods employed seem sound, as well as the simulations and the analyses performed. I find the results and the methods clearly explained. However, I have to say I do not find the main results as presented to be very surprising (from the Abstract: “we demonstrate 3 that these networks undergo a percolation transition as the number of links increases logarithmi- cally over time; the entangled particles form a giant cluster when the number of links exceeds a critical threshold.” [...] “Additionally, we find that particle’s opening angle significantly affects mechanical bonding and, consequently, the network structure”). The logarithmic increase of the number of links over time is interesting, but it is not motivated theoretically and the connection to other systems is a bit vague, as detailed below. I thus believe that the paper would benefit from a more detailed comparison with existing experimental and simulation results for similar systems.\nResponse: Besides the logarithmic aging, the following results are somewhat surprising:\n(1) There is no a priori reason to expect the entanglement networks of C-shaped grains to align with those of a continuum percolation model, since particle positions and orientations are gov- erned by mechanical interactions, gravity, and container constraints, rather than by perfect ho- mogeneous randomness. Nevertheless, the measured network properties are surprisingly similar.\n(2) Lifted clusters become highly elongated at the percolation threshold (Fig. 1b), and their ver- tical lengths can be predicted from the average eccentricity of the ring CP model at the corre- sponding degree (Fig. 3d).\nWe agree that comparisons with existing results for similar systems could highlight the novelty of our work.\nReviewer: In particular, a very similar work, studying rod-like particles instead of C-shaped particles via experiments and simulations, was recently published (Jung et al. [24]). It would be interesting to compare the results of the present work with those of [24]. In particular, the au- thors of [24] also find a percolation of the entanglement network which results into mechanical stabilization, despite their particles being straight rods, which I believe for the C-shaped particles would roughly correspond to a particle with very large diameter/thickness (D/d) and opening angle approaching 360 degrees. This seems to partially contradict the observation that mechan- ical stability becomes compromised for opening angles larger than 150 degrees (pag. 2: “When θ ≳ 150, the clusters are too fragile to be lifted as giant clusters”), so this is a point that I would like to see addressed.\nResponse: Thank you for the comment and for bringing the relevant work by Jung et al. (now cited as [24]) to our attention. While both works investigate the cohesion of granular materials through entanglement, there are several important differences:\n(1) The key distinction is the mechanism of cohesion. Jung et al. demonstrate that for extremely thin rods (aspect ratio α > 100), cohesion arises from friction between rods. This is directly ev- idenced by their Supplementary Movie S10, in which the removal of friction causes a complete loss of cohesion. In contrast, our C-particles are thicker, and thus their cohesion relies primarily on normal contact forces within topological links. If the C-particles were thinner, stronger fric- tional effects would raise the maximum θ for stability, reducing the fragile region in Fig. 1j.\n(2) Owing to this difference in cohesion mechanism, straight rods lack the ability to form a long, stable chain like the one shown in our Fig. 1b.\n(3) The very definition of “entanglement” differs, tailored to each mechanism: Jung et al. use the average crossing number (ACN), a metric typically applied to filamentous materials. The ACN quantifies the number of particle crossings when projected onto a plane, averaged over all ori- entations. While practical for filaments, ACN cannot directly capture topological links between 4 rings or C-particles. For example, the ACN of a dense pile of unentangled rings or C-particles (i.e., those not forming Hopf links) would still be nonzero, as their projected images contain many crossings. Similarly, for two overlapping rings lying in a plane, the ACN is similar whether they form a Hopf link or not. In contrast, our Hopf-link-based definition distinguishes these cases and serves as a reliable predictor of mechanical bonding, though it is not applicable to rods.\n(4) The methodology for assessing stability differs. Jung et al. begin lifting with a hook from the bottom of the material, which relies less on tensile stress. In contrast, our method of lifting a cluster from its top surface directly tests its tensile strength. If the rod cluster were instead lifted from the surface, it would likely break more easily.\nIn summary, the two studies are complementary rather than contradictory.\nTo clarify this, we have added\n• Line 154–163: We cite Jung et al. [24] to point out the impact of friction on mechanical stability.\nThis maximum θ for stability depends on particle thickness and friction coeffi- cient. 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. Reducing particle thickness or increasing the fric- tion coefficient would enhance cohesion, thereby raising the maximum stable θ and shrinking the fragile region in Fig. 1j.\n• Lines 778–784: We cite Jung et al. [24] as an example of a different entanglement criterion.\nEntanglement generally refers to geometric constraints between elongated or non- convex bodies that cause mechanical bonding. It can be quantified in various ways. For example, in filamentous materials such as polymers and rods, entan- glement can be statistically estimated through tube models [51], or measured by the average crossing number over all directions [24, 55].\nReviewer: An even more similar system to the one studied by the authors, consisting of U- shaped particles, has been studied experimentally and in simulations by Gravish et al. [4] (Ref. 3 in the present manuscript). How do the results of this study differ from those of Ref. [4]? Is there any fundamental qualitative difference between the two systems and their behaviors? If so, these differences should be clearly explained, because the two systems look very similar to me, which makes me question the novelty of this study.\nResponse: Thank you for this important point. While the systems appear similar, our study differs from Gravish et al. [4] in several key aspects that define its novelty:\n(1) A major advantage of the C-shape is that when we tune the entanglement capability through the opening angle, the particle diameter remains constant. This minimizes changes to the overall shape and size, enabling proper comparisons across opening angles. In contrast, varying the arm lengths of U-shaped particles significantly alters their overall shape and size, complicating direct comparison.\n5\n(2) We investigated cluster formation under vibration and the mechanical response to lifting, in contrast to Gravish et al. [4], who studied column collapse under vibration. This cluster-lifting approach provides a direct demonstration of tensile load-bearing capacity than column collapse experiments.\n(3) We employed a network-based framework and tested with percolation theory, which are ab- sent in Gravish et al. [4]. This framework allows us to reveal that entangled granular networks are explained by continuum percolation and display slow logarithmic growth in degree.\n(4) We performed Discrete Element Method (DEM) simulations, which provide detailed infor- mation on the internal structure arising from mechanical interactions and particle dynamics. In contrast, the Monte Carlo simulation employed by Gravish et al. is purely mathematical, lacking mechanical interactions or dynamics.\nWe anticipate qualitatively similar mechanical behavior for analogous shapes like U-shape. This expected similarity highlights the value of our framework; it is not specific only to the C-shape but provides a generalizable model and powerful analytical tools that is applicable to a broad class of entangled materials.\nTo clarify this, we have added\n• Lines 68–70: We emphasize our application of network-based framework.\nThese results establish network-based analysis as a promising framework for elu- cidating the structure and dynamics of entangled granular matter.\n• Lines 65–68: We emphasize our observation of logarithmic degree evolution.\nAdditionally, we observe that the mean degree grows logarithmically with vibra- tion time, mirroring slow relaxation in disordered systems.\n• Lines 55–59: We point out the benefit of using C-shape.\nThe 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.\nReviewer: Finally, since the authors mention kinetoplast DNA, I think it would be relevant to compare more in detail the results of this paper with those obtained from simulations of kineto- plast. I’m aware that the latter is a quasi-2D system whereas the system studied by the authors is 3D, but a comparison might still be informative.\nResponse:\nWe agree. Although kinetoplast DNAs are soft loops in a quasi-2D space, their degree distribu- tion is also Poisson-like (Fig. 6g in Ref. [21] and Fig. 2d in Ref. [20]). According to Ref. [21], this is because the entanglement is almost random, although the details differ from those of C- particles.\nTo clarify this, we have added\n6\n• Lines 196–200:\nA similar P(k) is observed in kinetoplast DNA networks, which are quasi-2D as- semblies of flexible loops [20, 21]. These results indicate that an underlying spa- tial randomness is shared by these systems, despite differences in geometry and dynamics.\nReviewer: Please find below some comments and questions that I believe should be addressed (in addition to those reported above) before this paper can be considered for publication.\nI believe that the operative definition of entanglement used in the paper (two particles are con- sidered entangled when the two circles passing through their centerlines form a Hopf link) should be reported in the main text for clarity (currently, the main text reports that “Two particles are defined as linked when the circles along their centerlines are intertwined”, but “intertwined” is not a well-defined concept).\nRelated to my previous point, I think that it might be useful to spend a few words on the mean- ing of entanglement. Due to the particular shape of the particles studied in this work, their en- tanglement can be straightforwardly defined using the topological concept of linking. However, this is not the case in other systems, in particular those made of slender objects such as poly- mers, living worms or rod-like particles (see [24]). I feel that the term “entanglement” tends to get used in a loose way to indicate rather different physical situations, and thus some clarification might help the readers who are not familiar with all the different definitions.\nResponse: Thank you for the comment. We agree that clarifying the definition of entanglement in the main text would be helpful for the readers.\nTo clarify this, we have added\n• Line 24–25: We cite three papers that explain “geometric cohesion”, another term for en- tanglement in mechanical contexts.\nEntanglement manifest macroscopically as deformation-resistant geometric cohe- sion [3–5]\n• Lines 47–50:\nIn 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.\n• Lines 97–102:\nFrom 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 θ is small.\n7\n• Lines 778–784 (‘Entanglement criterion’ in Methods):\nEntanglement generally refers to geometric constraints between elongated or non- convex bodies that cause mechanical bonding. It can be quantified in various ways. For example, in filamentous materials such as polymers and rods, entan- glement can be statistically estimated through tube models [51], or measured by the average crossing number over all directions [24, 55].\nReviewer: Is there any theoretical basis for the functional form reported in Eq. (2), or is this a purely empirical fit? Below this equation, the authors state that “similar logarithmic time evo- lution has been observed in various disordered systems”, reporting several references. However, from what I could see after a cursory look at these references, none of these study the time evo- lution of an entanglement network. What is, then, the justification to compare this logarithmic behavior to those observed in these papers, that concern very different quantities? I find the con- nection a bit vague and not justified well enough. I would invite the authors to compare their results to the time evolution of the entanglement network reported in Jung et al. [24] and Grav- ish et al. [4], for example. What comes closest to a microscopic theoretical explanation seems to be Ref. 32, although I am not sure of how to compare its predictions with the results here pre- sented.\nResponse: Thank you for your comment. We acknowledge that Eq. (2) is a purely empirical fit to our data, and its theoretical foundation remains an important open question.\nIn the absence of prior reports of logarithmic network growth that exactly match our observa- tions, we cited the closest available examples in Refs. [27–35]. These references are all glassy dis- ordered systems (polymers, granular media, DNA, crumpled paper and spin glass) characterized by slow logarithmic relaxation. This logarithmic evolution arises from sequential transitions be- tween metastable states separated by broadly distributed waiting times—often due to rugged en- ergy landscapes.\nAmong these, the most relevant observation is the logarithmic increase in density for vibrated spherical grains in Richard et al. [27], but our logarithmic increase is about the mean degree.\nNotably, these systems exhibit history dependence (memory effects). 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.\nWe note that the packing fraction of U-particles (Fig. 2a in Ref. [4], Gravish et al.) and the num- ber of contacts in rods (Fig. 4(c) in Ref. [24] Jung et al.) exhibit qualitatively similar time evolu- tion in that the slope decreases over time. However, their curves saturate to a constant within 20 seconds (Ref. [4]) and a few seconds (Ref. [24]), while we observe a continued logarithmic growth of the mean degree over much longer timescales. This discrepancy may be due to differences in the measured quantities and particle shapes.\nWe have added the following text in the main text:\n• Lines 216–239:\nTo our knowledge, this logarithmic evolution of entangled granular network has not been previously reported. The most relevant prior observation is the loga- rithmic density increase in vibrated spherical grains [27]. However, in our system, 8 the mean degree exhibits a clearer logarithmic trend whereas the density shows significantly larger fluctuations (Extended Data Fig. 8).\nEquation (2) is purely empirical, and its theoretical basis remains unclear. We speculate that insights may be gained by considering other glassy disordered sys- tems, 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–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 dy- namics apply to entangled granular systems remains a puzzle.\nReviewer: The authors find using finite-size scaling and other methods that the mean node de- gree at percolation is ⟨k ⟩ = 2.11 for infinite system size. I find this results to make intuitive c sense, as the network cannot percolate if the degree of each node is 2 or less, and a certain num- ber of nodes of degree 3 or more is needed. Have the authors compared this predictions to ana- lytical models for similar systems? For example, it is known that in systems of particles with a limited number of interaction sites (limited valence or “patchy” particles), the phase diagram of the system develops a phase-separated region as soon as the mean particle valence is above 2 – see for example (Sciortino and Zaccarelli, 2011). This result can be explained analytically by the Flory-Stockmeyer combinatorial theory (Stockmayer, 1943). I wonder if there is any connection to draw between this result and the one obtained here.\nResponse:\nWe appreciate the reviewer’s insightful comment, which allows us to clarify several important points. First, the term “network” refers to the entire system, as shown in Fig. 2a, which com- prises many disconnected clusters. The global mean degree ⟨k⟩ is averaged over all nodes, includ- ing a non-negligible population of isolated nodes with k = 0. Second, the percolation threshold is theoretically defined as the birth point of the largest cluster, where its fraction is still zero in the infinite-size limit [22]: S = 0 for ⟨k⟩ ≤ ⟨k⟩ and S ∝ (⟨k⟩−⟨k⟩ )β for ⟨k⟩ > ⟨k⟩ . At the 1 c 1 c c threshold, the network is still far from a fully connected graph and contains many isolated nodes. This allows the critical mean degree to be lower than 2, as seen in ER random networks where ⟨k⟩ = 1. Thus, the intuition that “a network cannot percolate if the average degree is less than c or equal to 2” does not apply to the theoretically defined ⟨k⟩ ; there is no general requirement c that ⟨k⟩ ≥ 2.\nc\nIn our ring CP model, we obtained ⟨k⟩ = 2.11, which is much larger than ⟨k⟩ = 1 for ER c c random networks. This reflects the threshold-increasing effect of spatial constraints in geomet- ric graphs [46]. Unlike ER networks, which live in an abstract, non-metric space where any pair of nodes can link, CP models are embedded in Euclidean space and connections are limited to nearby objects. As a result, more links per node are required to span the system in CP than in ER. Consistent with this, other CP models also exhibit higher percolation thresholds than ER: ⟨k⟩ = 2.74 for spheres [39] and ⟨k⟩ = 2.27 for disks in 3D space [40].\nc c\nTo clarify this, we have added:\n• Lines 175–181:\nTo investigate the entangled structure, we model the granular system as an N- node network comprising multiple disconnected clusters, 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.\n• Lines 387–393:\nThis suggests that the mean degree at the percolation threshold is ⟨k⟩ = 2.11 c when N → ∞, substantially larger than ⟨k⟩ = 1 of ER random networks [25] c and comparable to other CP models: ⟨k⟩ = 2.74 for spheres [39] and ⟨k⟩ = c c 2.27 for disks in 3D space [40]. This large ⟨k⟩ can be attributed to the threshold- c increasing effect of spatial constraints [46].\n• Lines 840–845 (‘Calculation of the excluded volume of rings’ section in Methods): We pro- vide the excluded volume of disks in 3D space, which is used to obtain ⟨k⟩ = 2.27 from c Ref. [40].\nThe excluded volume of disks (circular plates) of diameter D is vdisk = π2D3/8 [56], ex slightly larger than vring, as expected. Thus the reduced density (equal to the ex mean degree as N → ∞) of disks is η = π2nD3/8. Using this and the critical number density πn D3/6 = 0.9614 from Ref. [40], we obtain ⟨k⟩ ≡ η = 2.27 for c c c disks.\nReviewer: I could not understand how the effective reduced density η (pag. 6) is de- eff fined/calculated. The definition of this quantity and how it is calculated should be made clearer in the main text.\nResponse: We are grateful to the reviewer for this comment, which provides us with the oppor- tunity to clarify the definition and computation of η .\neff\nThe effective reduced density, η (⟨k⟩;N), is defined as the value of η in the ring CP model that eff yields the same ⟨k⟩ as the C-particle network. We use this value because in granular systems, the true reduced density η ≡ nv is not a good indicator of connectivity, as pointed out in line 408– ex 413:\nWe compute η by applying the empirical fit for η(⟨k⟩;N) from the CP model (Fig. 4a), using eff the measured ⟨k⟩ and N for the C-particles.\nIn Fig. 4a, the lines represent empirical fits of the form η = p 1 ⟨k⟩ + p 2 ⟨k⟩p3 (p 1 , p 2 , and p 3 are in Supplementary Table 2), used to calculate η (⟨k⟩;N). Thus, the effective reduced density is eff explicitly given by\nη (⟨k⟩) = p ⟨k⟩+p ⟨k⟩p3\neff 1 2\nwhere ⟨k⟩ is the measured mean degree in the C-particle network.\nTo clarify this, we have added:\n• Lines 350–352 (Fig. 4a caption): We provide the functional form of the empirical fitting curves, which are used to compute η .\neff\nLines represent empirical fits of the form η = p 1 ⟨k⟩+p 2 ⟨k⟩p3 (p 1 , p 2 , and p 3 are in Supplementary Table 2), used to calculate η (⟨k⟩;N) in e,f.\neff\n• Lines 413–418: We clarify the definition and computation methods of η .\neff\nTherefore, instead of η, we employ an effective reduced density, η (⟨k⟩;N), de- eff fined as the value of η in the ring CP model that yields the same ⟨k⟩ as the C- particle network. We compute η by applying the empirical fit for η(⟨k⟩;N) from eff the CP model (Fig. 4a), using the measured ⟨k⟩ and N for the C-particles.\nAdditional Changes\n• In the Introduction, we reordered of sentences in the 2nd paragraph to improve the logical flow with added sentences.\n• Lines 370–372: Removed “and” between equations to save space.\n• Lines 383–384, 398–403: To include the comparison of ⟨k⟩ with other CP models while pre- c servingreadability,wedividedtheparagraphintotwoandmodifiedthesentencesaccordingly.\n• Lines438–440: Weaddedasentencetoclarifytheresults’implicationsandsplittheparagraph into two for better readability.\nThese similarities demonstrate that continuum percolation theory applies to entan- gled granular structures.\nReferences\nSciortino, F. and Zaccarelli, E. (2011). Reversible gels of patchy particles. Current Opinion in Solid State and Materials Science, 15(6):246–253.\nStockmayer, W. H. (1943). Theory of Molecular Size Distribution and Gel Formation in Branched- Chain Polymers. The Journal of Chemical Physics, 11(2):45–55.\n12"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have prepared a careful reply and have reacted to all possible concerns raised by the two previous referees.\nThey have answered in a competent and convincing way. The recent paper by Jung, Plumb-Reyes, Lin and Mahadevan is now included in the discussion and I agree that although the topic is similar the samples used in the present manuscript are different: Jung et al use straight rods which do not entangle a priori while the present paper is on C-particles. And there are no contradictions between the two papers.", "reviewer_3": "The authors addressed all the concerns that I had raised, and modified the manuscript accordingly, clarifying the methodology and strengthening the connections with the existing literature.\nI find this work appealing because it introduces a simple experimental system which, complemented by theory and simulations, provides some interesting insights on the physics of entanglement in granular media.\nDespite somewhat similar systems of C-shaped particles having already been studied in simulations and experiments (see refs. 4 and 19 in the paper), this work is the first one to combine experiments, discrete elements simulations, network analysis and continuum percolation theory to study a system where the degree of entanglement can be experimentally tuned (via the aperture angle) in a controlled way, minimizing changes to the particle shape and size.\nSome of the findings, like the prediction of the percolation threshold (critical density), the development of a continuum percolation theory for entangled (linked) rings, and the observation of a logarithmic increase of the number of links over time, are interesting and could bear important connections with other physical systems of entangled and disordered materials.\nThus, I believe that this work will be of general interest to researchers working on granular and disordered materials, and to the broadly defined soft matter community. For these reasons, I recommend it for publication in Nature Communications.", "reviewer_4": "I think the paper is quite interesting, carefully written and the responses to the reviewers are excellent and well thought out. I advocate publication."}, "author_response": ""}]

Decision letter