William Y.
Wang
*a,
Stephen J.
Thornton
a,
Bulbul
Chakraborty
b,
Anna R.
Barth
a,
Navneet
Singh
a,
Japheth
Omonira
a,
Jonathan A.
Michel
c,
Moumita
Das
c,
James P.
Sethna
a and
Itai
Cohen
ad
aDepartment of Physics, Cornell University, Ithaca, New York 14853, USA. E-mail: wyw6@cornell.edu
bDepartment of Physics, Brandeis University, Waltham, Massachusetts 02454, USA
cSchool of Physics and Astronomy, Rochester Institute of Technology, Rochester, New York 14623, USA
dDepartment of Design Technology, Cornell University, Ithaca, New York 14853, USA
First published on 26th March 2025
We study how the rigidity transition in a triangular lattice changes as a function of anisotropy by preferentially filling bonds on the lattice in one direction. We discover that the onset of rigidity in anisotropic spring networks on a regular triangular lattice arises in at least two steps, reminiscent of the two-step melting transition in two dimensional crystals. In particular, our simulations demonstrate that the percolation of stress-supporting bonds happens at different critical volume fractions along different directions. By examining each independent component of the elasticity tensor, we determine universal exponents and develop universal scaling functions to analyze isotropic rigidity percolation as a multicritical point. Our crossover scaling approach is applicable to anisotropic biological materials (e.g. cellular cytoskeletons, extracellular networks of tissues like tendons), and extensions to this analysis are important for the strain stiffening of these materials.
Rigidity percolation has been well studied in isotropic networks under different bending and stretching constraints.6–8 However, less is understood about anisotropic networks. Anisotropic networks occur naturally in many biological tissues, such as bone9,10 and tendon,11,12 both of which are composed of collagen fibrils oriented strongly in some direction. Anisotropy is expected to alter the onset of rigidity percolation, which is known to be sensitive to details of the bond distributions. For example, previous work has shown that including structural correlation within isotropic networks can result in significant changes in the critical bond occupation threshold for rigidity percolation.13,14 Furthermore, studies have shown that straining a percolated but floppy network, such as by shearing it in one direction, can drive a rigidity transition.15 For example, straining the network preferentially along the maximum extension axis activates bond stretching, which rigidifies the network. Finally, previous computational studies have also modeled anisotropic networks through an anisotropically diluted triangular lattice and found that the onset of rigidity agrees well with a simple Maxwell constraint counting argument and that the system can be approximated using an effective medium theory (see Zhang et al.2). Missing from these analyses, however, are detailed investigations of whether the critical exponents and scaling functions characterizing the rigidity transitions depend on these details of the bond distributions. Measuring the critical exponents and how they depend on the bond distributions is critical for determining whether such mechanical phase transitions are in the same universality class, which informs the relevant physics governing these transitions. Understanding this physics is an important step for learning how to control these transitions in materials ranging from biological tissues to synthetic fiber networks.
Here, we will examine the critical exponents and locations of phase boundaries in anisotropic networks. Though our model is quite similar to ones previously investigated,2 we find slightly different critical exponents for the isotropic rigidity percolation transition. Remarkably, we also discover that when we tune away from isotropy, the network exhibits two rigidity transitions. The first is associated with the component of the strain tensor for stretching along the preferential filling direction. The second rigidity transition is associated with the remaining components of the strain tensor. As such we find that the isotropic rigidity transition is a multi-critical point from which anisotropic rigidity transitions emanate (Fig. 1). It is unclear whether this intermediate, infinite anisotropy elastic phase will be present for generic lattices (see Section 4), but similar infinite anisotropy mechanical responses should arise in strain-stiffened networks (see Fig. 4).
![]() | (1) |
We apply strains εxx, εyy, and εxy to measure the elastic coefficients Cxxxx, Cyyyy, and Cxyxy directly. To measure Cxxyy, we perform a bulk compression and subtract out the energetic contributions from the independently measured Cxxxx and Cyyyy moduli. We introduce strains to the lattice by applying the proper transformation matrix to the positions of each of the nodes. For example, to stretch the network in the horizontal direction (to apply strain εxx), the following matrix is applied:
![]() | (2) |
In order to measure the energetic costs of our imposed strains in the disordered lattices, we minimize the central-force energy functional over the positions of the nodes. To capture the linear response, we truncate to leading order in the displacement of vertices:
![]() | (3) |
We then individually perform scaling analyses for each independent component of the elasticity tensor. Starting with a rigid network, we repeatedly remove bonds and minimize the energy for each strain until the network becomes sufficiently “floppy” in all directions; here, a network is considered floppy in a particular direction if the corresponding modulus falls below a threshold of Gmin ≡ 10−8, which is the simulation tolerance. The moduli as a function of p were averaged for each system size and orientation strength pair (L,r), sampled over 102–104 random seeds. This procedure allows us to perform a scaling analysis of each component of the elasticity tensor separately.
〈pc〉L − p∞c ∼ L−1/ν. | (4) |
Our analysis determines that p∞c = 0.645 ± 0.002, depicted by the black dot along the horizontal axis of Fig. 1, and ν = 1.3 ± 0.2. The location of the threshold at p∞c indicates a small deviation from the naïve Maxwell constraint counting, which states that the 2D triangular lattice to have 2 constraints per site, p = 2/3 of the lattice must be occupied. The deviation from this prediction in our measured value of p∞c is similar to what is found in other works.6 Importantly, we find that in the isotropic system (r = 1) all the moduli share the same threshold value p∞c.
Next, we perform a finite-size scaling analysis of each component of the elasticity tensor, which admits a scaling
![]() | (5) |
〈pijklc(r)〉L − p∞,ijklc(r) ∼ L−1/ν′ | (6) |
Remarkably, we find that for r > 1, the transition for the Cxxxx modulus is distinct from the threshold values for the other components of the elasticity tensor. We find that the phase boundary for Cxxxx bends towards lower values of p with increasing r (Fig. 1). The transition curves for the other moduli appear nearly identical for each system size and bend towards higher values of p with increasing r. This separation indicates a region for which the network is only rigid when strained along the preferred bond orientation direction. We verify these phase boundaries are the same at isotropy but distinct for r > 1 by measuring the separations between the histograms of pc values for different components of the elasticity tensor as a function of system size (see Appendix E). Thus, the transition from a floppy to a rigid phase occurs in two stages when the system is anisotropic.
We note that for finite system sizes, Cijkl never truly vanishes on average because if there are enough bonds remaining to create a rigid system-spanning truss , this configuration occurs with finite probability. However, we still expect that in the intermediate regime between 〈pxxxxc(r)〉L and 〈pijklc(r)〉L, finite systems will have an elasticity tensor with a principal axis in the x-direction due to our formulation of anisotropy.
Next, we test the scaling behavior of the elasticity tensor components near the critical points. In principle, including anisotropy could introduce corrections to scaling that bend the phase boundary in a trivial way, leaving all critical exponents the same as in the isotropic system. We tested this scenario, by fixing r > 1 and attempting to collapse the data near the relevant pc(r) values for each component of the elasticity tensor, keeping faniso = fiso and ν′ = ν, but found very poor collapse. This poor collapse suggests that the anisotropic phase transition is in an entirely different universality class, with different values for the critical exponents. We thus conjecture that the vicinity of r = 1 should be analyzed as a crossover scaling between two distinct critical points.
Inspired by renormalization group approaches, we analyze our data using crossover scaling functions expected to be valid in the vicinity of the isotropic critical point:
![]() | (7) |
We first estimate ζ by examining the shape of the phase boundaries in the infinite system away from isotropy. Specifically, as shown in Appendix F, because the arguments of the universal scaling function are invariant scaling combinations, this phase boundary must occur at a fixed value of X/Y1/ζ = δp/(r − rc)1/ζ (with corrections to scaling), so that the separation between the two phase boundaries in Fig. 1 scales as (r − rc)1/ζ. From this estimate based on the shape of the phase boundaries, we find ζ = 0.25 ± 0.1.
![]() | ||
Fig. 3 Crossover scaling of anisotropic rigidity percolation. Each independent elastic modulus, a function of the variables (p,r,L), collapses onto a two dimensional sheet when plotted against scaling variables X ≡ (δp)L1/ν and Y ≡ (r − rc)Lζ/ν. (left) The scaling function ![]() ![]() |
Many other quantities share this crossover scaling ansatz and allow for independent estimates of ζ. The widths of the histograms of pc values as we tune away from isotropy are also amenable to a crossover scaling analysis, with the variable Y collapsing these widths (Appendix F). From this scaling collapse, we similarly estimate ζ = 0.25 ± 0.1 and find nice collapse (Fig. 13).
This estimate for ζ and the scaling ansatz in eqn (7) can be used to collapse the elasticity tensor components for 250000 simulations consisting of anisotropy values ranging between 1.0 ≤ r ≤ 2.0, bond occupation values ranging between 0.6 ≤ p ≤ 0.68, and system sizes L ranging between 30 ≤ L ≤ 500. We show a two variable scaling collapse for the Cxxxx and Cyyyy moduli in Fig. 3. We find excellent collapse of each independent modulus onto a two dimensional sheet. The overlaid data points consist of a portion of the data used to produce the sheet and indicate various slices of constant Y: Y = 0, Y = 0.66, and Y = 1.65. The Y = 0 curve (black) shows the finite-size scaling for isotropic systems and is identical to that shown in Fig. 2. We observe similarly excellent collapse at the two higher values of Y (Fig. 15 and 16). For ease of visualization, we also include height contours projected onto the X–Y plane. The height contours for Cxxxx curve toward lower values of X as the scaling variable Y increases, reflecting the fact that the corresponding phase transition curves toward lower values of p as r increases. The other moduli, such as Cyyyy, show the opposite systematic behavior, tending towards higher values of X for increasing Y.
The critical exponents determined thus far, including ζ, are properties of the isotropic rigidity percolation critical point. We are also interested in the anisotropic critical exponents, such as the critical exponent with which each modulus vanishes with p in an infinite anisotropic system (fanisoijkl). The universal scaling function for the crossover between the isotropic and anisotropic critical points in principle contains all of these anisotropic critical exponents in its singularities in various asymptotic regimes. However, it is quite difficult to get reliable high-precision fits to generic two-variable scaling functions that include precise information about their singularities. We instead independently estimate fanisoijkl by examining the largest system sizes of simulations performed at r = 1.2, 1.5. In addition, we analyze the data for Y = 0.66, 1.65.
We find that the critical exponents for Cxxxx are distinct from those found for the isotropic system while those for the other elasticity tensor components cannot be distinguished from those found for the isotropic system:
![]() | (8) |
Exponent | Estimate |
---|---|
ν | 1.3 ± 0.2 |
f iso | 2.2 ± 0.3 |
ζ | 0.25 ± 0.1 |
ν aniso | 1.2–3.2 |
f aniso xxxx | 4.0 ± 1.0 |
f aniso ijkl (ijkl ≠ xxxx) | 2.2 ± 1.0 |
It was a surprise to us that the rigidity percolation transition for the isotropic lattice broke up into several transitions when it became anisotropic. First, obtaining multiple transitions is contrary to the naïve usage of Maxwell counting to determine the location of the rigidity transition. As shown in the insets of Fig. 1, the rigid modes become anisotropic, and span horizontally before they span vertically.
Second, our results are fundamentally different from connectivity percolation, where regardless of the value of r there can only be one transition. The left inset in Fig. 1 shows several horizontal stress-supporting chains spanning the network. The critical point at which Cxxxx first becomes non-zero presumably separates a phase where there are no stress-supporting chains from one where there are a finite density of such chains. In regular percolation, two such paths connecting the system horizontally that are separated by any finite distance will have a finite probability per unit length of being connected by bonds extending in the vertical direction. Hence, for ordinary percolation, as soon as one crosses the horizontal percolation point in an anisotropic system, it must percolate in the other directions as well. This argument suggests that lattices with bending stiffnesses and angular springs, which are believed to become rigid at the connectivity percolation threshold,5 lack this intermediate phase. However, in typical situations where bending stiffnesses are much weaker than stretching stiffnesses, a remnant of this intermediate phase should be measurable even when bending is included.
In retrospect, we should have expected separate transitions in central force rigidity percolation on regular lattices. Maxwell counting tells us when the number of zero modes can vanish in the absence of states of self stress, but does not tell us whether the zero modes couple to a given mode of deformation.18,19 Straight lines of bonds supporting stress in a large system, when connected vertically, may only contain contributions to the stress that grow quadratically (i.e. non-linearly) in the strain: the length of a beam connecting (x,y) to (x′,y + ε) grows as ε2, suggesting the corresponding linear elastic modulus is 0. A similar nonlinear response to infinitesimal strains is found to stabilize hypostatic jammed packings of ellipsoidal particles,20 in violation of simple constraint-counting arguments. The perfect square lattice has no Cxyxy shear modulus, and the perfect hexagonal lattice has no non-zero moduli except the bulk modulus – why should anisotropic random lattices not possess separate transitions?
Since generic lattices have no straight lines of bonds, it is not clear to us whether they will have intermediate phases with infinite anisotropy in linear response. (Indeed, a non-generic lattice which supports horizontal strain under tension has straight lines of bonds that buckle under compression making the effective modulus under compression still zero). Even in this case, we expect crossover scaling to be an important component of the analysis of anisotropic generic lattice networks. Our analysis suggests that the scaling of anisotropic generic networks could exhibit a crossover between distinct rigidity transition universality classes.
The difference between generic and non-generic lattices is no longer crucial after finite deformations; under tension a node connecting only two bonds will naturally straighten out to 180° and not support perpendicular stresses to linear order. Along these lines it might be interesting to investigate whether separate rigidity transitions arise in floppy isotropic random lattices under finite deformation. Preliminary results suggest strain stiffening leads to rigidity along the extension axis without rigidity under some of the other modes of deformation (see Fig. 4).
What kind of rigidity critical points do we expect? The Cxxxx transition where horizontal stress-bearing chains first arise could be self-similar (with a single diverging correlation length), but could also be self-affine (with the vertical spacing between chains diverging with a different power on the stress-supporting side than the rigid cluster lengths diverge on the floppy side, akin to directed percolation, where connectedness lengths are controlled by different critical exponents parallel and perpendicular to the preferred direction21–23). Whether the three other moduli become non-zero simultaneously or separately in this model is not numerically resolved yet, but one expects that more complicated anisotropies (say a 3D model with brick-like symmetry) will allow for separate transitions for those moduli as well.
Finally, it would be interesting to consider whether the results presented here have any bearing on biological networks. Many biological tissues, including bone, tendon, muscle, and blood vessels, have an extracellular matrix that exhibits preferential alignment, giving rise to anisotropic elastic moduli. It may be possible to understand the extracellular matrices of these tissues as anisotropic networks which have crossed the rigidity percolation transition in the stiff direction, but not the other directions. A quantitative description of these tissues may require analyzing the effects of second-order constraints, changing the contact number distribution, and exploring the crossover between this pair of zero-bending stiffness transitions and a bending-dominated regime. Moreover, cells are able to control how matrix elements are generated. Cells may generate networks with many different rigidity transitions to tune between, where the particular way matrix elements are laid down biases the thresholds of the different components of the elasticity tensor. As such, the results presented here could have profound implications for understanding more complicated networks in many biological systems.
We start by picking a random number seed, and then assigning a (uniformly chosen) random number si between 0 and 1 to each bond. The bonds are assigned and then sorted based on a key, ki, which corresponds to the value of bond occupation fraction p for which the bond would be added based on the anisotropy parameter r:
![]() | (9) |
![]() | (10) |
To handle periodic boundary conditions, we split the Hessian into two parts: Hpbc, which is computed using only the bonds that span across the network, and Hin with bonds which do not. The energy is therefore computed as:
![]() | (11) |
(Hin + Hpbc)![]() ![]() | (12) |
An affine displacement is used as an initial guess. The sparsity structure allows matrices to be stored in compressed sparse row format, reducing memory usage and improving the speed of operations. A Cholesky factorization for (Hin + Hpbc) is computed16 and used as a preconditioner. For large system sizes, we use GPUs to accelerate numerical computations, such as matrix factorizations and matrix-vector products.
We note that for systems that have under-constrained nodes, the null space of (Hin + Hpbc) has a non-zero dimension; as such we do not consider the non-affinity parameter (which sums the squared displacements from an affine transformation) as a method of analysis or for extracting critical exponents.
pijklc(L) ∼ ρijklL, ρijklL ∈ Δ[0,1]. | (13) |
Histograms of pijklc are plotted in Fig. 5. As the system size increases, the distributions become increasingly sharp and the means shift systematically. In the limit of infinite system size, we assume that each distribution converges to a delta function about p∞c.
For each system size, we compute both the means 〈pijklc〉L and standard deviations σijklL of each distribution. We expect systematic shifts in the means and standard deviations (i.e., the first and second moments) to scale as a power law with respect to L governed by a single critical exponent ν:
![]() | (14) |
We perform a joint non-linear least squares fit,24 with p∞c and ν the same for all curves (we assume they are equal for each modulus at isotropy). Fig. 6 depicts 〈pijklc〉L (left) and σijklL (right) as a function of L. The fits give an estimate of p∞c = 0.645 ± 0.002 and ν = 1.3 ± 0.2.
![]() | ||
Fig. 6 (left) Mean of pc and (right) standard deviations of the pc distribution at isotropy as a function of system size. The fitted curves match those in eqn (14) and the fitted lines on the right figure have slope −1/ν. |
Furthermore, we find a universal scaling function for the distributions ρijklL with respect to our previously defined scaling variable X ≡ (δp)L1/ν:
![]() | (15) |
We find that X collapses the density functions, with the resulting histograms shown in Fig. 7.
![]() | ||
Fig. 7 Universal scaling of rigidity distributions at isotropy for each independent modulus. The histograms all collapse when plotted against the scaling variable X. |
![]() | ||
Fig. 8 Unscaled modulus data across various system sizes at isotropy (r = 1). The smaller system sizes have a higher probability of becoming rigid at lower values of p. |
We find the exponent fiso by considering the largest available system size (L = 500) and performing a least-squares fit of the modulus, finding a range of fiso as we vary pc slightly. With fiso = 2.2, p∞c = 0.646, and ν = 1.3, we plot the data against the proposed scaling variables and find a nice collapse for all independent moduli, shown in Fig. 9. The number of samples we average over ranges from 104–102 for system sizes L ∈ [30,200] and 20 samples of L = 500.
Using the value of fiso = 1.4 ± 0.1 and ν = 1.4 ± 0.2 quoted in Broedersz et al.,6 we find best collapse with p∞c = 0.65 shown in Fig. 10, which gives good collapse for lower values of system size, but does not collapse the modulus for our largest system.
![]() | ||
Fig. 10 Universal scaling function of all independent elastic modulus using previously reported exponents.6 The data for the range L ∈ [30,200] comparable to the earlier work does give a good collapse. Having the larger system size (L = 500) explains why we find different exponents. |
Sample-to-sample, there are lattices that can support rigidity in some shear directions but not others. There are two basic scenarios. (1) In the case where including anisotropy ends up simply giving analytic corrections to scaling that bend a single phase boundary, all moduli will vanish at the same location at L = ∞, but the amplitudes of the finite-size effects may be different. Singling out the Cyyyy modulus for the sake of comparison, this means that
〈pijklc〉L − 〈pyyyyc〉L ∼ L−1/ν. | (16) |
![]() | (17) |
(2) In the case where including anisotropy leads to genuinely new critical phenomena and a pair of phase transitions, the finite-size effects of the mean and standard deviation are controlled by the finite-size scaling exponent of each anisotropic rigidity transition νaniso. If we split into a pair of phase transitions, then the spreads σijkl of each distribution will narrow with increasing L, but the separation of the means is asymptotically constant as L → ∞. This would make
![]() | (18) |
We begin by performing this measurement of the separations between the distributions of pc for all moduli at the isotropic transition, where all moduli vanish at the same location in p as L → ∞ (Fig. 12 (left)).
![]() | ||
Fig. 12 Separation of rigidity percolation threshold mean 〈pc〉 from that of the Cyyyy modulus as a function of system size at isotropy (left) and anisotropy (right). |
As expected, this measure of the separation in the means is flat for the isotropic case as a function of system size, confirming that these moduli vanish at the same asymptotic location and that the finite-size effects controlling the mean and the standard deviation have the same systematic dependence on L.
When we perform the same analysis for the anisotropic case (r = 1.5), we see systematic growth in this measure as a function of system size (Fig. 12 (right, blue)), suggesting that the locations of pc for different moduli are genuinely different in the thermodynamic limit. This is moderately strong quantitative evidence for the information that can roughly be seen by eye in Fig. 11.
Furthermore, we can consider the standard deviations of the rigidity percolation threshold distributions (as in Appendix C) for r > 1 and consider a scaling function with respect to our scaling variable Y = (r − 1)Lζ/ν
![]() | (19) |
We find that standard deviation is best collapsed with ζ = 0.25 ± 0.1, shown in Fig. 13. The determined value of the exponent ζ appears to collapse σL1/ν for all moduli except for Cxxxx at larger values of the scaling variable. We also note the nice overlap of data performed at different values of (r,L) but the same value of the scaling variable Y.
We can in principle use information from these collapse plots to make a prediction for the value of the finite-size scaling exponent close to the anisotropic phase transition νaniso. First, we note that if we fix r > 1 and send L → ∞, the spread in the distributions of pc will narrow as σ ∼ L−1/νaniso, as the finite-size effects are (at large enough system sizes) controlled by the critical exponents of the anisotropic critical point. In the scaling function for the distributional spreads of pc, this corresponds to the asymptotic limit Y → ∞. Forcing the asymptotics of the numerically determined crossover scaling function to agree with the asymptotics expected at the anisotropic critical point will give us a prediction for νaniso.
Suppose this scaling function has asymptotic behavior at large Y. Then in the limit L → ∞ with r > 1 fixed,
σL1/ν ∼ Yα ∼ Lζα/ν and σ ∼ L−1/νaniso | (20) |
![]() | (21) |
We estimate the exponent with which each modulus vanishes at their corresponding anisotropic phase transition (fanisoijkl) by considering our largest system size at two values of r away from isotropy (r = 1.2 and r = 1.5). The location of the phase transition is determined by averaging the value of p with which each modulus for each lattice becomes rigid.
To search for which value of ζ best collapses the two variable scaling function, we obtain simulations of constant Y across various system sizes. We test a value of ζ by first fixing our largest system size and then solving for the value of r as a function of system size that results in the same value of Y. Fig. 15 and 16 show a scaling collapse of all the moduli with ζ = 0.25 and with a constant value of Y = 0.66 and Y = 1.65, respectively. The scaling function for each modulus vanishes at a given value of W, which we denote as Wijklc(Y), as it is dependent on the value of Y. We note that further away from isotropy (Y = 0), there are additional corrections to scaling, resulting in a worse collapse at higher values of Y.
![]() | ||
Fig. 15 Scaling collapse of all the moduli at constant Y = 0.66. The estimate of Wc is shown with the dashed black line. |
![]() | ||
Fig. 16 Scaling collapse of all the moduli at a constant Y = 1.65. The estimate of Wc is shown with the dashed black line. |
We again obtain an estimate of the scaling exponents in Fig. 17 and 18 by plotting the rescaled moduli as a function of distance from Wc. The plots suggest crossover for the Cxxxx modulus, with fanisoxxxx governing the behavior for lower values of W − Wc and fiso for the high W − Wc regime. Furthermore, the other three moduli appear to vanish with a critical exponent indistinguishable from fiso within our estimated error bars.
![]() | ||
Fig. 17 Estimate of anisotropic scaling exponent of each modulus at a constant Y = 0.66. The simulation data is the same as found in Fig. 15. |
![]() | ||
Fig. 18 Estimate of anisotropic scaling exponent of each modulus at a constant Y = 1.65. The simulation data is the same as found in Fig. 16. |
This journal is © The Royal Society of Chemistry 2025 |