Frictional rigidity percolation: A new universality class and its superuniversal connections through minimal rigidity proliferation.

We introduce two new complementary concepts, frictional rigidity percolation and minimal rigidity proliferation, to help identify the nature of the frictional jamming transition as well as significantly broaden the scope of rigidity percolation. To probe frictional rigidity percolation, we construct rigid clusters using a (3,3) pebble game for sliding and frictional contacts first on a honeycomb lattice with next-nearest neighbors, and second on a hierarchical lattice. For both lattices, we find a continuous rigidity transition. Our numerically obtained transition exponents for frictional rigidity percolation on the honeycomb lattice are distinct from those of central-force rigidity percolation. We propose that localized motifs, such as hinges connecting rigid clusters that are allowed only with friction, could give rise to this new frictional universality class. And yet, the distinction between the exponents characterizing the spanning rigid cluster for frictional and central-force rigidity percolation is small, motivating us to look for a limit in which they are identical, i.e., a search for mechanisms of superuniversality. To achieve this goal, we construct a minimally rigid cluster generating algorithm invoking generalized Henneberg moves, dubbed minimal rigidity proliferation. For both frictional and central-force rigidity percolation, these clusters appear to be in the same universality class as connectivity percolation, suggesting superuniversality between all three transitions for such minimally rigid clusters. These combined results allow us, for the first time in rigidity percolation, to directly compare two universality classes on the same lattice and to highlight unifying and distinguishing concepts of rigidity transitions in disordered systems.


I. INTRODUCTION
At the heart of every rigidity transition is the emergence of a spanning rigid cluster-an entity of interconnected bonds that are rigid with respect to each other. For disordered systems, the starting point of choice has become randomly diluted spring networks with central-force interactions [1][2][3][4]. As bonds are randomly diluted from a triangular lattice, either a regular one or with slightly randomized lattice points (a generic lattice), the system goes from rigid with a nonzero shear modulus to floppy without this feature [5][6][7][8]. Underlying this mechanical phase change is the transition from a system with a spanning rigid cluster to a system without one, as identified by the combinatorial (2,3) pebble game [9]. The location of this rigidity transition occurs approximately at isostaticity where the number of degrees of freedom (d.o.f.) are frozen out by the number of force-balance equation constraints, known as Maxwell constraint counting [10].
The rigidity transition in the central-force, randomly bond-diluted triangular lattice was numerically found to be a continuous one with a correlation length exponent, ν ¼ 1.21 AE 0.06, an order parameter exponent β ¼ 0.18 AE 0.02, and a fractal dimension of the spanning rigid cluster, d f ¼ 1.86 AE 0.02 [6,8,11]. These exponents differ slightly from two-dimensional connectivity percolation where ν ¼ 4=3, β ¼ 5=36, and d f ¼ 91=48 [12]. Despite the small difference in exponents, it was eventually argued that centralforce rigidity percolation (RP) is in a separate universality class since there are nonlocal effects in terms of how rigid clusters grow that differ from connectivity percolation [6,8]. Meanwhile, Bethe lattices with no loops are amenable to analytical treatment and demonstrate that the spanning rigid cluster at the transition is not fractal [13,14]. To add to the complexity, numerical simulations of three-dimensional lattices with central-force interactions indicate a discontinuous rigidity transition as well, in contrast to the two-dimensional case [15]. Finally, central-force models with next-neighbor springs can exhibit hybrid rigidity transitions with both continuous and discontinuous features [16]. With this rather varied set of phase transitions when changing just the type of lattice, the general solution to the central-force RP problem is far from clear, if it is even possible.
Rigidity percolation with bond-bending forces adds another "dimension" to the problem [17][18][19][20][21]. Numerical simulations of two-dimensional systems measuring elastic properties suggest that bond-bending forces drive the transition into a different universality class [22]. However, since there is currently no pebble game for bond-bending forces even in two dimensions, a direct comparison to central-force rigidity percolation in terms of ν, β, and d f has yet to be made.
Particle packings also undergo a rigidity transition as a function of the packing fraction [23][24][25] and additionally feature contact network rearrangements, unlike the randomly diluted spring networks. The rigidity transition in such systems has been labeled the jamming transition, where the system moves from a zero to a nonzero bulk and shear modulus with increasing packing fraction, suggesting the emergence of a spanning rigid cluster [24]. This suggestion was made explicit by first extracting the contact network of a two-dimensional frictionless (i.e., central forces only) particle packing at jamming. Second, constructing the rigid clusters from this network via the (2,3) pebble game shows that at the onset of rigidity or jamming every particle participating in the contact network is part of one rigid cluster; i.e., the spanning rigid cluster is bulky at the transition [26].
In frictional particle packings, jamming is even less understood. Recently, two of us [27] developed a new pebble game algorithm incorporating both particle translations and rotations to compute rigid clusters for networks abstracted from two-dimensional frictional particle packings. We applied it to molecular dynamics simulations of frictional particle packings at a fixed packing fraction experiencing slow shear, and determined rigid clusters at constant strain intervals as the packing repeatedly goes through jamming or unjamming. The size of the largest rigid cluster indicates a continuous transition, and so does the observation of a roughly power-law distribution of rigid cluster sizes near the jamming transition [27]. Within the spanning rigid clusters, we found regions of floppiness, showing partial rigidity. The floppy regions are also physically relevant as the nonaffine motion of the particles is smaller inside the rigid cluster compared to outside, and the pressure is higher. We present four such rigid cluster images close to the transition for four different values of the friction coefficient μ in Fig. 1. The major open question arising from this recent study of frictional rigid clusters is whether or not the rigidity transition is actually continuous or not, as a continuous rigidity transition would be very different from the frictionless case. Unfortunately, the molecular dynamics simulation could not tune the system to be arbitrarily close to the rigidity transition, and so it remains difficult to assess the nature of the transition using such simulations.
Here, we ascertain the nature of the rigidity transition in frictional systems with the introduction of frictional RP, that is, the study of rigid clusters constructed via the frictional (3,3) pebble game, on randomly diluted lattices with bonds denoted as either frictional or sliding. We study two types of lattices: the honeycomb lattice with nextnearest neighbors (NNN) and a hierarchical lattice. The former can be studied numerically, while the latter is amenable to analytical calculations. The results for both lattices can then be compared to results for central-force RP to arrive at the first direct comparison ever between different types of forces on the same lattice. In generalizing (a) (b) (c) (d) FIG. 1. Rigid clusters in simulated frictional packings under slow shear. Four snapshots from the molecular dynamics simulation showing partially rigid systems close to the frictional rigidity transition. In black is the largest rigid cluster, while floppy regions are colored gray. Blue, green, and red disks correspond to three, two, and one leftover pebble, respectively. The range of partial rigidity decreases with increasing friction coefficient μ, or equivalently, increasing q. (a) μ ¼ 0.2, with a transition at q ¼ 0.78 and average coordination number z ¼ 3.35, (b) μ ¼ 0.3, with a transition at q ¼ 0.86 and z ¼ 3.15, (c) μ ¼ 0.5, with a transition at q ¼ 0.95 and z ¼ 3.0, and finally (d) μ ¼ 10 with a transition at q ¼ 1.0 and z ¼ 2.8. The last value is due to the large number of contactless particles (rattlers) in the packing, visible in blue.
rigidity percolation, we can now take a step back to look at a broader set of problems to be more readily able to identify special cases in a forest of rather eclectic trees. Our findings can be summarized as follows.
(1) Both our numerical results on the honeycomb lattice with next-nearest neighbors and our analytical results on a hierarchical lattice exhibit a continuous rigidity transition with a distinct set of clustermerging mechanisms not available to central forces to lead to a new universality class distinct from central-force rigidity percolation.
(2) Since now all three universality classes (connectivity percolation, central-force rigidity percolation, and frictional rigidity percolation) have only slightly distinct order parameter exponents and fractal dimension of each spanning rigid cluster at the transition, we hypothesize that within a subset of minimally rigid clusters, there may be superuniversality across these transitions. We construct a new algorithm for both types of rigidity percolation, dubbed minimal rigidity proliferation (MRP), to directly create minimally rigid spanning clusters whose properties are indeed similar to clusters in connectivity percolation. More specifically, our numerical finite-size scaling results show that while the correlation length exponent ν and the nonspanning rigid cluster size exponent τ are rather different between the two types of forces, the order parameter exponent β and the fractal dimension of the spanning rigid cluster at the transition d f are rather similar. We conclude that central-force RP and frictional RP represent two distinct universality classes based on the rather different values of ν and τ. We analyze the mechanics of merging rigid clusters and find, for example, that hinges between rigid clusters constitute rigid connectors for frictional RP, whereas they are not for central-force RP. The different rigid cluster merging mechanisms between the two cases, along with different types of rigid clusters, potentially account for the difference in ν and τ, as evidenced by our analytical calculations on a hierarchical lattice. However, the nearness of the other two exponents pertaining to the spanning rigid cluster (β and d f ) motivates us to search for the possibility of superuniversality in which the two universality classes collapse onto the universality class of connectivity percolation. Superuniversality occurs when, despite a symmetry difference, two different models are in the same universality class, i.e., have the same exponents. An example of superuniversality is the three-state Potts model and the four-state Potts model in the presence of quenched disorder belonging to the same universality class, despite having different spin symmetries [28].
We explore the superuniversal possibility only within the subset of rigid cluster configurations that are minimally rigid. By (1) using concepts from invasion percolation [29] used to build only spanning clusters in connectivity percolation and (2) extending the Henneberg moves [30] to grow a minimally rigid network with central forces only, we construct a new algorithm to directly grow a minimally rigid network with frictional forces, which we dub minimal rigidity proliferation. We argue that the order parameter exponent in minimal rigidity proliferation is indeed the same across connectivity percolation, central-force RP, and frictional RP. This search for superuniversality within a subset of rigid cluster configurations will motivate searches for other subsets in which the same phenomenon arises, and then determine how to navigate between those subsets that are distinct and those that are not. This more nuanced approach moves beyond the standard paradigm of centralforce rigidity percolation and connectivity percolation as being completely distinct.
The structure of this paper is as follows. We first describe frictional rigidity percolation on the honeycomb lattice with next-nearest neighbors and present our results. We then describe frictional rigidity percolation on a hierarchical lattice and discuss how our analytical calculations can help bolster the numerical results on the honeycomb lattice. Moving on, we present minimal rigidity proliferation and argue how its strategic bond occupation method is different enough from our earlier random bond occupation to lead to superuniversality amongst the different models. We then conclude with a discussion of the implications of our results. A graphical abstract is provided in Fig. 2 and serves as a roadmap for the paper.

II. FRICTIONAL RIGIDITY PERCOLATION
A. Honeycomb lattice with next-nearest neighbors

Model
To motivate the model we begin by reviewing Maxwell constraint counting in two-dimensional frictionless packings (or central-force spring networks) with N particles (vertices) and average coordination number z [10]. The total number of d.o.f. is 2N, while there are ½ðNzÞ=2 forcebalance constraints. When the number of d.o.f. is equal to the number of force-balance constraints, the system is minimally rigid, i.e., where the number of global rigid body translations and rotations has been subtracted out since they are trivial. This equation yields the critical average coordination number z c ¼ 4 (as N → ∞) for the onset of rigidity. Much numerical work with frictionless particle packings has shown that this counting is an extremely good method to determine the rigidity transition point [31][32][33]. No states of self-stress are observed in particle packings at the transition, since the values of the purely repulsive forces are uniquely determined by the boundary conditions at this point, such that a more involved constraint counting approach is not needed [34]. For two-dimensional frictional particle packings, some contacts are below the Coulomb threshold with the magnitude of the tangential force less than the magnitude of the repulsive, central force times the friction coefficient. At such contacts, two particles can only rotate and translate with respect to one another just as a gear does, and these are denoted as frictional contacts. There are also contacts at the Coulomb threshold in which two particles slide with respect to each other. For these sliding contacts, the magnitude of the tangential force is set by the magnitude of the repulsive, central force; i.e., there is only one constraint. We distinguish between these two types of contacts by denoting q to be the probability of having a frictional contact, with 1 − q denoting the probability of having a sliding constact; i.e., if q ¼ 1, all contacts are frictional. Then, performing the Maxwell constraint counting as above, since each particle has 3 translational and rotational d.o.f., there are 3ðN − 1Þ total d.o.f. (subtracting out the trivial global d.o.f. in which there is no relative motion between the particles). Moreover, the interparticle forces yield f½zð1 þ qÞN=2g total constraints. We, therefore, arrive at the minimal rigidity criterion, or where q denotes the probability of having a frictional bond. For N → ∞ and q ¼ 1, all bonds are frictional and z c ¼ 3.
(2) describes a line of transition points interpolating from z c ¼ 3 to z c ¼ 4 as the ratio of frictional to sliding bonds changes. This method of counting is now known as generalized isostaticity [35,36]. Such bounds are indeed observed in experiments [37]. Note that increasing the friction coefficient μ increases q, and in addition, that one cannot smoothly interpolate between frictional and frictionless packings as one cannot smoothly interpolate between 2 and 3 local d.o.f.
To construct a lattice model for frictional particle packings, we consider a honeycomb lattice with additional NNN bonds. This modified honeycomb lattice has a maximum coordination number of z max ¼ 9. We define p as the probability of bond or contact occupation. The reason we employ the honeycomb lattice with next-nearest neighbors is because we can explore geometry to determine whether or not it is relevant for determining the nature of the phase transition. We do so by constructing and studying two different models for bond occupation; see Fig. 3. For the first model, we fully occupy the honeycomb backbone such that p ¼ 1=3 and then occupy the additional NNN bonds occupied randomly such that p ≥ 1=3. We dub this first strategy of bond occupation HC1. We also implement a second strategy of bond occupation in which the bonds, both nearest-neighbor and next-nearest neighbor, are occupied at random, which we dub HC2. For HC1, since the honeycomb backbone is fully occupied, the central forces on each particle can be balanced, which is required by local mechanical stability in frictionless, but not frictional, packings. Therefore, HC1 will allow us to more readily

HC1 HC2
FIG. 3. Schematic frictional rigidity percolation models HC1 and HC2. These random networks are constructed by either adding next-nearest neighbor (NNN) bonds to an occupied honeycomb lattice (HC1) or adding random first-and secondneighbor bonds (HC2), all with probability p. Double bonds denote frictional or gearlike bonds, which occur with probability q, while single bonds denote sliding bonds.
compare with the geometry of frictionless packings in order to see how frictionless differs from frictional. The frictional bonds are then randomly assigned with probability q, and periodic boundary conditions are implemented. Now we address the frictional versus sliding bonds for this lattice model. Since frictional bonds randomly occur with probability q, using Eq. (2) the critical occupation probability predicted by Maxell constraint counting is Equation (3) tells us how p c depends on q, therefore, denoting a phase transition line between floppy and rigid phases just as in generalized isostaticity.
Once the frictional and sliding bonds have been identified, we construct a constraint network in which frictional bonds below the Coulomb criterion are denoted as double bonds in the constraint network and sliding bonds at the Coulomb threshold are denoted as single bonds in the constraint network. We then play the (3,3) pebble game on this constraint network in which the first number denotes the number of local d.o.f. and the second number denotes the number of trivial global d.o.f., which does depend on boundary conditions. However, Ref. [27] found that changing the number of trivial global d.o.f. from 3 to 2 due to periodic boundary conditions did not significantly affect the rigid cluster analysis for both frictional and frictionless particles, and so we stick with the (3,3) and (2,3) pebble games.
We illustrate the (3,3) pebble game algorithm using several very simple constraint networks in Fig. 4. A more detailed explanation can be found in Appendix A. Examples of the rigid clusters we find below, at, and above the rigidity transition for both HC1 and HC2 are shown in Fig. 5. To compare frictional RP with central-force RP, we complement our analysis with an approach where any double bond is converted to a single bond and a (2,3) pebble game is played since each site contains now only 2 d.o.f. Examples of the corresponding rigid clusters are shown in Fig. 6. Finally, we implement finite-size scaling analysis to quantify the transition for HC1 and HC2 for different q's for the frictional (3,3) pebble game and for the central-force (2,3) pebble game.

Results
Spanning probability.-We first identify the location of the rigidity transition by determining whether or not there exists at least one spanning rigid cluster in the x or y direction as both p and q are varied. We do this for all four variants, HC1 and HC2 for both the frictional (3,3) game and the frictionless (2,3) game. For HC2, we study q ¼ 0.5 and q ¼ 1.0, the two extreme values of q. For q ¼ 1, isostaticity predicts p c ð1Þ ¼ 1 3 . For the HC1 model, this is identical to the initial occupation of the honeycomb lattice backbone, and for this regular lattice, we expect one spanning rigid cluster with a unity probability of spanning, i.e., a first-order transition. Therefore, for HC1, we study q ¼ 0.5 and q ¼ 0.7 and do not explore the limit q → 1 since q ¼ 1 is a special case. Figure 7 plots the probability that the system contains at least one spanning rigid cluster as a function of p for different system lengths L.  The probability of spanning for HC1 at p c ðqÞ for both q ¼ 0.5 and q ¼ 0.7 is approximately 0.6. Since this value is significantly less than unity, our findings suggest a continuous transition for the onset of the spanning rigid cluster. Typically, the value of probability of spanning at the transition is not a universal quantity and depends on details of the model. For the frictionless version of HC1, shown in Fig. 7(c), the crossing point is more difficult to determine, but we estimate it to be near 0.448(1).
For the HC2 version of the model (bottom row of Fig. 7), we again find crossing points near the predicted generalized isostaticity counting since the formula also applies to this variant of the bond occupation. Since there is no ordered honeycomb backbone that is initially occupied, we explore both the lower and upper bounds of q, i.e., q ¼ 0.5 and q ¼ 1.0 [Figs. 7(d) and 7(e)]. Our results can be found in the first column of the table in Fig. 9. We note that there is greater discrepancy of the estimated p c from generalized isostaticity for HC2 than for HC1. We also note that the probability of spanning at the transition (the crossing point) now differs between q ¼ 0.5 and q ¼ 1.0, which does not imply a different universality class because the crossing point depends on details of the lattice. The frictionless version of HC2 is plotted in Fig. 7 Correlation length.-The correlation length ξ quantifies how two distant particles or sites interact. In a continuous transition, the correlation length diverges at the transition, while near the critical point, ξ ∼ jp − p c ðqÞj −ν on either side of the transition, where ν is the correlation length exponent. In a finite-size system and near the transition, ξ is replaced by the system length L. For each realization, after this replacement, the system has a finite-size critical point p L c ðqÞ when the system contains a spanning rigid cluster, with jp L c ðqÞ − p ∞ c ðqÞj ∝ L −1=ν . Since the location of the transition fluctuates for each realization, we therefore obtain a distribution of finite-size critical points as observed in Fig. 7. The standard deviation of this distribution Δ yields a measurement of the correlation length exponent [38]. More precisely, Using error function fits to the data in Fig. 7  shows the width of the transition for the six different variations of the honeycomb lattice model. Spanning rigid cluster.-We now study the properties of the spanning rigid cluster using P ∞ , the fraction of occupied bonds in the spanning rigid cluster. Figure 8(a) shows P ∞ for increasing p for different system lengths for HC1 with q ¼ 0.5. We note that P ∞ at p c ð0.5Þ decreases as the system size increases, which, again, suggests that rigidity transition here is continuous. The behavior of this curve just above the critical point p c ðqÞ is described by the order parameter exponent β with P ∞ ∼ ½p − p c ðqÞ β for p ≥ p c ðqÞ, for an infinite-size system. As long as L ≫ ξ, the equation applies and P ∞ ∼ ξ −β=ν . However, when L ≪ ξ, the length scale will be set by L such that P ∞ ∼ L −β=ν . We therefore introduce a universal scaling function fðL=ξÞ that interpolates between these two regimes, or with fðL=ξÞ ¼ ðL=ξÞ −β=ν for L ≪ ξ and a constant for L ≫ ξ. The universal scaling functionfðL=ξÞ can be obtained by rescaling as is done in Fig. 8(b) for q ¼ 0.5, with ν ¼ 1.56 and β ¼ 0.18 used as fitting parameters to obtain the optimal collapse. This estimate for ν is consistent with our previous measurement for the same q from ΔðLÞ in Fig. 8(e). The collapse supports the notion of a continuous rigidity transition. We implement the same protocol for the remaining cases to look for a continuous rigidity transition.
With the exception of the frictionless version of HC1, shown in Figs. 8(c) and 8(d), the order parameter exponent does not vary too much between the different models, as summarized in Fig. 9(c), though given the smallness of β, it is more difficult to measure as precisely as ν. The small value of β ≈ 0.07 for the frictionless version of HC1 perhaps suggests that this model is similar to the square and kagome lattices with next-nearest neighbors studied in Ref. [16]. There, a hybrid transition was found, where the onset of the spanning cluster was discontinuous, but with a diverging correlation length, though the correlation length exponent appeared to be unity. In addition to the order parameter exponent, one can also measure the fractal dimension of the spanning cluster to determine whether or not it is, indeed, fractal. To test for this possibility, the fractal dimension is determined by measuring the number of bonds in the spanning rigid cluster M as a function of system length such that M ∼ L d f . In Fig. 8 .05, so we observe little change in the fractal dimension with q, at least for these system sizes, provided q is not close to unity. Similar fractal dimensions were found for both the frictional and frictionless versions of the HC2 version of the model and are listed in the table in Fig. 9(c).
Nonspanning rigid clusters.-In connectivity percolation, one typically investigates the nonspanning cluster size distribution, defined as the number of finite clusters of size s per lattice site or bond, or n s [38]. At the transition, n s ∼ s −τ , where τ is the cluster size exponent. For connectivity percolation, we have the inequality τ > 2 strictly. How can we understand this result? We start with P ∞ s¼1 sn s ðpÞ þ P ∞ ðpÞ ¼ p. Since P ∞ ðp c Þ ¼ 0 for connectivity percolation, then P ∞ s¼1 sn s ðp c Þ ¼ p c at the transition. Using the assumption that n s ∼ s −τ and converting the sum to an integral, τ > 2 for convergence to a finite value, i.e., p c .
In rigidity percolation, the situation is more complex because there are nonspanning rigid clusters, spanning rigid cluster(s), and floppy regions. If both the nonspanning rigid cluster size distribution and the floppy cluster size distribution are power laws independently at the transition, each exponent associated with the respective size distribution should be greater than 2. On the other hand, if one of the exponents is less than 2, that would suggest a natural cutoff for that type of cluster and the more tenuous structure could still facilitate a continuous transition. If both exponents are less than 2, then perhaps this aspect of the transition is discontinuous.
As detailed below, we cannot rely on hyperscaling relations for RP. Therefore, we keep track of the nonspanning rigid clusters only and posit that their size distribution also behaves as a power law at the transition with exponent τ. If τ < 2, then there is presumably a characteristic cutoff for the nonspanning rigid cluster sizes at large enough sizes with coupling to the floppy regions perhaps driving the continuity of the transition. Figure 9(a) shows the probability for having a nonspanning rigid cluster of size s as p is increased through the transition point for HC1 with q ¼ 0.5 on a log-log scale. Below the transition point, there are many small rigid clusters in the system. As p increases, they merge into larger ones and the distribution broadens to approach a linear function on a log-log scale; the downward trend of the tail is due to finitesize effects. We obtain τ ¼ 1.90 AE 0.03 < 2 from a linear fit to the relevant part of the curve. Above p c ðqÞ, the spanning rigid cluster "swallows" the nonspanning rigid clusters and ultimately, as p is increased far beyond the transition point, there is only one spanning rigid cluster left. We have measured τ for the six different cases and find a persistent difference between the frictional and frictionless case in that τ > 2 for the frictionless cases, while τ < 2 for all frictional versions [see Figs. 9(a) and 9(b)] indicative of rather different ways the rigid clusters merge and grow in the two cases.

Rigid cluster merging mechanisms
The results of our finite-size scaling analysis are summarized in the table in Fig. 9. We also include the exponents for central-force rigidity percolation using the (2,3) pebble game on the triangular lattice (denoted as TL) and for connectivity percolation (denoted as CP) on the triangular lattice for comparison. For the frictional versions implementing the (3,3) pebble game, we find that HC1 and HC2 appear to be in the same universality class, with the exception of the special case of HC1 at q ¼ 1, in which a discontinuous transition emerges as discussed earlier in Sec. II. We also conclude that exponents associated with HC2 and the (2,3) pebble game are in the same universality class as the exponents for central-force rigidity percolation on the triangular lattice obtained about 20 years ago. On the other hand, we find that the exponents associated with HC1 and the (2,3) pebble game are potentially more related to the square lattice plus braces (i.e., next-nearest neighbors) in which a hybrid transition was found [16], so that this case is special, just as HC1 with q ¼ 1 is special.
So while our (2,3) pebble game results are consistent with prior central-force rigidity percolation results, our new frictional rigidity percolation compels us to ask the question, what mechanism(s) could drive frictional rigidity percolation and central-force rigidity percolation to be in distinct universality classes? To begin to answer this question, we ask the following question: How do two rigid clusters combine to form one larger rigid cluster?
Unlike in connectivity percolation, in rigidity percolation two independently rigid clusters cannot become one rigid cluster by joining via a single bond. For a frictional rigidity percolation example with q ¼ 1, consider two triangles, which are individually rigid. If Even one double bond and one single bond connecting to the two triangles generates one merged rigid cluster. The two spatially distinct bonds leading to rigid clusters merging in frictional rigidity percolation does not hold for central-force rigidity percolation. For central-force rigidity percolation, at least three bonds are needed to merge two independently rigid clusters. One bond fixes the distance between the two rigid clusters, the second bond the relative rotation between them, and the third the shearing, provided the bonds are not all parallel with respect to each other. In the frictional case, the one double bond between the two rigid clusters fixes both the distance and the relative rotation. Rigid hinges are another means by which two rigid clusters can merge at a point and still be rigid. In central-force rigidity percolation, hinges consisting of single bonds between rigid clusters are always floppy, and so rigid hinges cannot exist. However in frictional rigidity percolation, this is not the case, at least for a hinge composed of all double bonds; see Figs. 10(b) and 10(d).
We conclude that two double bonds and the rigid hinge (composed of double bonds) are distinct means of propagating rigidity in frictional rigidity percolation that do not occur in central-force rigidity percolation. Both frictional motifs are more spatially localized than their central-force analogs even in the absence of floppy regions. The presence of floppy regions indeed complicates matters, as they can become rigid as well due to the merging of rigid clusters. This rigidification of floppy regions can then trigger other rigid clusters to merge, until the rigidity cascade is complete; see Figs. 10(c), 10(e), and 10(f). It is these very local motifs of connecting rigid clusters that can participate in a rigidity cascade to give rise to nonlocal, or distant, rigidity due to the addition of one bond. While both types of models contain such an effect, we believe that the zerodimensional rigid cluster connector (the hinge) and the three-versus two-bond rigid cluster connectors could potentially account for the difference in exponents reminiscent of correlated percolation models such as kcore percolation, where the k ¼ 2 behavior is very different from the k ¼ 3 behavior [39]. To more thoroughly understand how the rigidity propagates through the system for frictional RP as compared to central-force RP is combinatorially tricky. In the next section we study a hierarchical lattice where we are easily able to perform such a task.
In contrast, to the ν and τ exponents, the structure of the ultimate spanning rigid cluster appears to not differ as dramatically between the two cases. Specifically, we observe little variation in β for the different models studied. We also observe little variation in d f for both HC1 and HC2 studied with just central forces and with frictional forces. We address this finding after presenting our hierarchical lattice results.

B. Hierarchical lattices
While we have presented predominantly numerical results so far for frictional RP and argued for a distinct universality class from central-force RP, one exactly solvable RP model is RP on hierarchical lattices. We can therefore analytically determine if indeed the central-force RP is in a different universality class than frictional RP. To do so, we first review prior results using the (2,3) pebble game with central forces only and then generalize to the frictional version.

Review: Central forces only
It has been previously shown that central-force RP transitions in such lattices exhibit a continuous rigidity transition [40,41]. To understand this finding, let us start with the generation of a particular hierarchical lattice known as the Berker lattice [41]. Given two points and a bond as in Fig. 11(a), replace the bond with some base structure to generate a first-generation hierarchical structure. This replacement continues ad infinitum to arrive at a network with an infinite number of sites between two initial points. For this particular lattice, the nth generation contains 8 n bonds (with the exception of n ¼ 0). To embed this lattice in two dimensions, the bond length decreases with each generation.
To analyze rigidity in this hierarchical lattice, assume each bond has a probability p < 1 to be occupied. In the n ¼ 0 graph, the probability of having a spanning rigid cluster between the two ends (black dots) is p 0 ¼ p. In the n ¼ 1 graph, the probability of being rigid between two ends can be found by subgraph counting: If all eight bonds are occupied in the n ¼ 1 network, there is a spanning rigid cluster between the two ends. The probability of such a structure existing is p 8 , while the probability for any bond belonging to the spanning rigid cluster is 1. All other subgraphs that contain a spanning rigid cluster, as determined by the (2,3) pebble game, and their respective probabilities are listed in Fig. 11(b). Summing up all ways of having a spanning rigid cluster between the two ends of the n ¼ 1 graph, we obtain Given the hierarchical structure of the lattice, it is trivial to generalize this relation to from which we can solve for a fixed point, p c ¼ 0.9446, as the system approaches the thermodynamic limit, i.e., p nþ1 ¼ p n . In Fig. 11(c), p n as a function of p n1 is plotted for the first four generations. We observe that the curves cross at p ¼ p c and p n will converge to a step function which jumps from 0 to 1 at p c as n goes to infinity. Meanwhile, we use P R ðpÞ to denote the probability for a bond to belong to the spanning rigid cluster. The recurrence relation for P R ðpÞ is and near p c , λ ¼ 0.9554 < 1, demonstrating that the probability of a bond belonging to the spanning rigid network will approach zero as p approaches p c . This trend suggests a continuous transition. Expanding about p c in both Eqs. (7) and (8) leads to ðp nþ1 − p c Þ ¼ λ 1 ðp n − p c Þ and P nþ1 ðpÞ ¼ λ 2 P n ðpÞ, such that λ 1 ¼ b 1=ν , λ 2 ¼ b −β=ν , and λ 3 ¼ b d f , where b is the length rescaling factor from one generation of the hierarchical lattice to the next. For the Berker lattice, λ 3 ¼ 8. We can therefore determine β ¼ − logðλ 2 Þ= logðλ 1 Þ and νd f ¼ logðλ 3 Þ= logðλ 1 Þ, which are both quantities that are independent of b, resulting in β ¼ 0.078 and νd f ¼ 3.533.

Frictional forces
Let us now consider a "frictional" hierarchical lattice with double and single bonds to denote frictional and sliding contacts. Double bonds are introduced at random with probability q. When double bonds are taken into account, they affect subgraph constraint counting in the hierarchical lattice as we now play the (3,3) pebble game to determine whether or not a subgraph has a spanning rigid cluster. Since there is an increased number of possible subnetworks in the frictional case given that the occupied bonds can be either double or single bonds, let us first discuss the q ¼ 1 case. Here, there are several additional types of subgraphs containing a spanning rigid cluster that were not allowed in the central-force case, as shown in Fig. 12(a), which we can easily identify. For instance, one of the subgraphs is not allowed in the central-force case because it contains a hinge structure. In the frictional case, the frustrated loops of odd numbers of vertices, or "gears," prevent rotation. These additional rigid subgraphs contribute an additional 16p 6 ð1pÞ 2 to the probability of having a spanning rigid cluster above the central force case. Therefore, the counting for q ¼ 1 leads to the recursion relation In the limit n → ∞, we find the unstable fixed point p c ðq ¼ 1Þ ¼ 0.8533, in addition to two stable fixed points at p ¼ 0 and p ¼ 1. Moreover, we can compute P R;nþ1 ðp; q ¼ 1Þ to arrive at such that λ 2 ðq ¼ 1Þ ¼ 0.3511. Since λ 2 ðq ¼ 1Þ < 1, the rigidity transition is continuous with d f ν ¼ 3.181. This value is indeed distinct from the central-force value on the same lattice, thereby indicating two different universality classes. Now we consider q < 1. After keeping track of what subgraphs are rigid between the two black circles in Figs. 11(b) and 12(a) in the presence of both double and single bonds, we obtain With q ¼ 1, the unstable fixed point occurs at p ¼ p c ðq ¼ 1Þ ¼ 0.8533 with the two stable fixed points at p ¼ 0 and p ¼ 1. Therefore, p n will converge to a step function as n → ∞. However, for 0.8465 < q < 1, both the unstable and nonzero stable fixed points, p lower and p higher , respectively, are smaller than 1, so that p n will converge to a step function which jumps at p lower from 0 to p higher . The reason that p higher is not unity in these cases is because p denotes a double or a single bond such that when q ¼ 1, p ¼ 1 translates to all double bonds; however, when q < 1 and p ¼ 1, p higher depends on the ratio of double to single bonds. When q ¼ 0.8465, p lower ¼ p higher , which means when q ≤ 0.8465, p n will always converge to zero and rigidity transition will vanish entirely, showing that the existence of a rigidity transition in this hierarchical lattice very much depends on q. We also compute d f ν (for q > 0.8465) and find that its value depends on q, as shown in Fig. 12(b). And while there is one value of q at which the two d f ν values are the same, β may also be different. So, again, we find analytical evidence for two distinct universality classes. In addition, the fact that the correlation exponent depends continuously on q is not necessarily unique, as has been found in Ising models on hierarchical lattices; see, e.g., Ref. [42]. This sensitivity is presumably due to the special nature of the hierarchical lattice, as detailed in Appendix B.

III. MINIMAL RIGIDITY PROLIFERATION
For our numerical studies on the honeycomb lattice, we found that the order parameter exponent β and the fractal dimension of the spanning rigid cluster at the transition d f are not very distinct between frictional and central-force RP. Are they in fact the same, signaling features of superuniversality for the structure of the spanning rigid cluster at the transition? Or is it the case that these exponents are indeed different but the distinction is small, making it hard to detect studying finite-sized systems? If we can find a model where the order parameter exponents are actually the same-a model where the bonds are strategically placed as opposed to randomly placed, for example-this strengthens the potential for superuniversality rather than relying purely on numerical analysis, which has its limitations.
So let us now explore more explicitly connections between frictional RP and central-force RP via a subset of rigid cluster configurations using an algorithmic approach rather different from finite-size scaling. As will become clear below, connectivity percolation also enters the picture, since if we can construct spanning rigid clusters in the same way as geometrically connecting clusters, then we have evidence for superuniversality across all three models.
We first review invasion percolation, which is motivated by the problem of one fluid displacing another from a random, porous medium [29]. More importantly for us, invasion percolation allows one to create a spanning cluster on a lattice that has the same properties as a spanning cluster in connectivity percolation. Next, we review the Henneberg moves [30], which are used to grow a large minimally rigid network (a Laman graph) from a small minimally rigid network in the central force case. We then extend the Henneberg moves to include frictional forces and ultimately unify the two concepts, invasion percolation and Henneberg moves. The final algorithm that we introduce, minimal rigidity proliferation, allows us to grow minimally rigid networks that span a frictional system, and only grow such networks.
Invasion percolation is a modified version of connectivity percolation where the spanning cluster grows along the path of smallest weights, with the following algorithm.
(1) Assign uniformly distributed random numbers ranging from 0 to 1 to bonds on a lattice as their weights.
(2) Occupy an initial bond, and create a list of all its neighbors. This list creates a boundary of bonds.
(3) Occupy the bond from the list that has the smallest weight. (4) Update the list so that it contains all unoccupied nearest neighbors of occupied bond.
(5) Repeat 3 and 4, until the occupied cluster spans the entire lattice. The above algorithm reduces to the Leath algorithm [43], which creates the spanning cluster for bond connectivity percolation for p > p c in the following limit: Instead of occupying the boundary bond with the smallest weight, all boundary bonds whose weight is less than p are accepted into the cluster, and then the boundary list is updated. The algorithm terminates when there are no bonds on the boundary with weights less than p. This modification from invasion percolation to the Leath algorithm does not affect the large-scale structure of the spanning cluster; i.e., they remain part of the same universality class [44].
Let us also review the Henneberg moves for building a minimally rigid network with central forces only. A minimally rigid graph in this case is also known as a Laman graph. Minimal rigidity occurs when the d.o.f. match the constraints and there are no rendundant bonds, as determined through a (2,3) pebble game. Starting from such a network GðN; N B Þ with N B bonds and N sites, one can extend it using two basic Henneberg moves, as illustrated in Fig. 13 (top).
(1) Add one site and two bonds between this site and two points in G, then G 0 ðN þ 1; N B þ 2Þ is the new minimally rigid network (type I move), (2) or add one site and three bonds between this site and three prior sites in G, then delete a prior bond between two of the selected three prior sites (type II move). Both moves simultaneously add 2 d.o.f. and two constraints, which results in a minimally rigid graph by induction. Now we generalize, for the first time, Henneberg moves for the (3,3) pebble game in order to propagate minimal rigidity. We focus on two cases: q ¼ 1=2 and q ¼ 1. For q ¼ 1=2, we consider only a type I move by adding a site and then adding three bonds, one double bond and one single bond (see Fig. 13, middle). This move perpetuates minimal rigidity since no dependent constraints are introduced. For q ¼ 1, i.e., all double bonds, we consider two type II moves in series, if you will, by adding two sites, where the first site connects to two existing sites and the second new site must attach to the initial new site as well as an older site. Then, any one of the double bonds between the first new site and either old site is removed, though not the bond between the two new sites, to preserve minimal rigidity (see Fig. 13, bottom).
Having discussed invasion bond percolation and the "growing" of minimal rigidity via generalized Henneberg moves, we are now ready to introduce minimal rigidity proliferation. We first focus on (2,3) minimal rigidity and then address (3,3) minimal rigidity.

A. Central-force case
To create a spanning minimally rigid cluster as defined by the (2,3) pebble game, we combine the Henneberg move type I and invasion bond percolation in the following algorithm (see Fig. 14 for an illustration).
(1) Assign uniformly distributed random numbers ranging from 0 to 1 to bonds on the honeycomb lattice as their weights. (2) Begin by occupying a random triangle between three closest sites and create a list of all nearest and nextnearest neighbor bonds of these sites. (3) Determine the sum of the weights of any two bonds from the sites on the list that join at one site and the existing sites in the graph, then find the smallest sum and occupy those two bonds. (4) Update the bond list such that it contains any unlisted nearest and next-nearest neighbor bonds of the newly added site. (5) Repeat 3 and 4, until the graph spans the lattice. Though the graph is grown by adding two bonds at a time, as opposed to one, we still expect that this process will fall under the connectivity percolation universality class. Why? Because adding two bonds (with their additive weights) at a time involves a simple rescaling of time in which two bonds are added in one time step as opposed to two bonds in two time steps. With this simple rescaling of time, we do not expect the structure of the spanning cluster to change between standard invasion percolation and minimal rigidity proliferation; i.e., they are the same. Also note that in MRP there are no floppy regions and so there is no distant, or nonlocal, rigidity.
In Fig. 15 we show an example of a spanning minimally rigid cluster on the honeycomb lattice with NNN bonds using minimal rigidity proliferation. We measure the fractal dimension by computing the number of sites in the cluster M as a function of total number of sites N, as shown in   . 16(a), and obtain M ¼ N 0.958 . Since N ¼ L 2 , this leads to a fractal dimension of the spanning rigid cluster at the critical point of d f ¼ 1.916, which is consistent with connectivity percolation. Figure 16(b) shows P ∞ , the fraction of the system in the spanning cluster, converging to zero when system becomes infinitely large, suggesting a continuous transition just as with connectivity percolation.

B. Frictional case
Since we have extended the Henneberg moves to the (3,3) pebble game for q ¼ 1=2 and q ¼ 1, we can generalize minimal rigidity proliferation to the frictional case. For the q ¼ 1 case, two Henneberg type II moves are made in sequence to arrive at one growth step. With bond removal, it is not immediately clear that the minimally rigid cluster growth results in the same cluster structure as the q ¼ 1=2 case, and so we leave this for future study. However, since the type I move for q ¼ 1=2 corresponds precisely to the type I move for the central-force case just with one single and one double bond, we expect the same configurations as above. Thus, both central-force percolation and frictional RP collapse to an identical construction in this case. Since we have already argued in the central-force case that adding two bonds at a time is a simple rescaling of time from adding one bond at a time, we also expect the q ¼ 1=2 frictional process to be in the same universality class as connectivity percolation via transitivity. In other words, within this subset of growing minimally rigid configurations, we expect superuniversality to emerge: All three universality classes collapse into one all with the same exponents.
Should we have expected this result? Indeed, we should have because this strategic bond occupation does not contain floppy regions, so there is no nonlocal rigidity, and the minimal rigidity constraint maps to adding either one or two bonds at a time. In other words, MRP is a simpler model than random bond occupation given that there is no nonlocal rigidity. Interestingly, transfer matrix methods (not focusing on minimally rigid clusters) argued that connectivity and central-force percolation were in the same universality class, but their results were later discounted [5,6]. We now perhaps have some understanding as to why some exponents appear to be quite close in value in that one can add floppy, or redundant, bonds in some perturbative manner and interpolate between the two limits.

IV. DISCUSSION
We have now expanded the notion of rigidity percolation to include friction in two dimensions with the extension of  the (2,3) pebble game to the (3,3) pebble game and the incorporation of double bonds representing contacts below the Coulomb threshold. In doing so, we have uncovered a new universality class in the realm of rigidity percolation, namely that of frictional RP, which is directly compared with central-force RP on the same lattice. Such a direct comparison between two universality classes has not been possible until now. By expanding the scope of RP, the direct comparison presented here should help to formulate a more general framework for rigidity transitions just as there exists a general framework for spin systems, with and without disorder, to understand how a transition in the three-state Potts model is in a different universality class from the one in the Ising model. We make this direct comparison between central-force RP and frictional RP on honeycomb lattices with additional next-nearest bonds. We find different correlation length and nonspanning rigid cluster size distribution exponents ν and τ, respectively, between the two cases, but a statistically similar order parameter exponent and fractal dimension of the spanning rigid cluster at the transition. Given the different ν and τ, we propose that local motifs, such as two double bonds and a rigid hinge composed of double bonds, are ways to connect rigid clusters in frictional RP that are distinct from central-force RP. Neither construct is rigid in central-force RP, and additional supporting bonds are necessary. The less strict rigid cluster connection mechanisms in frictional RP compared to central-force RP potentially drive the distinction between universality classes. For the hierarchical lattice, not only are there two different universality classes, frictional and central-force RP, that can be shown analytically, the exponents also depend continuously on the fraction of double bonds.
Motivated by the small difference in order parameter exponent in central-force and frictional RP, we combined Henneberg moves [here extended to the (3,3) pebble game] and invasion percolation to construct another new model, minimal rigidity proliferation, that can be implemented for both types of forces. The rigid cluster in MRP grows in a simple fashion, unlike in RP where rigid clusters surrounded by floppy regions can lead to a rigidity cascade. With minimal rigidity proliferation, there are no floppy or even redundant bonds-the spanning rigid cluster is built in a "clever" way, which is to be contrasted with the tuning by pruning approaches [45], where springs are removed while conserving the bulk modulus, for example, and the jamming graph approach [46], where minimally rigid clusters follow the geometric constraint of local mechanical stability. With the strategic bond growth in minimal rigidity proliferation, the order parameter exponent is the same across connectivity, central-force RP, and frictional RP. This would be the first time superuniversality is observed in RP in a way that goes beyond transfer matrix methods [5]. Our work also suggests that looking at minimally rigid configurations-a subset of all possible configurations within RP-represents a new way of viewing phase transitions in the sense that nested within two distinct universality classes there could be an underlying superuniversality establishing deeper connections between the classes than previously thought.
Since frictional RP was devised to explore the nature of the jamming transition in frictional particle packings, this work compels us to make a rather strong claim that the rigidity transition in frictionless particle packings with purely repulsive central forces is of a different nature than the rigidity transition in frictional particle packings. In fact, the frictionless case with purely repulsive central forces may indeed be a very special case because even rigid cluster analysis of particle packings with both attractive and repulsive central forces indicates a continuous transition [47]. In frictionless packings, there are no redundant bonds, which makes the constraint counting rather straightforward. However, in frictional packings, redundant bonds emerge such that the constraint counting is more intricate and, therefore, perhaps non-mean-field. It will be interesting to apply the frictional (3,3) pebble game to experimental frictional particle packings to test the applicability of our approach as well as to compare the rigid clusters with dynamical matrix calculations. And very recently, the frictional (3,3) pebble game has been applied to frictional packing derived beam networks to predict fracture locations near the brittle-ductile transition [48].
Finally, we are currently exploring a limitation of the (3,3) pebble game [49]. Specifically, if there are four particles forming a square and all four contacts are below Coulomb threshold, then we have a square with all double bonds (like the middle image in Fig. 4 without the diagonal bond). From the (3,3) pebble game perspective, this configuration is floppy, and there is one floppy mode where the particles are in a gearing motion. However, these four particles are rigid under strain, since the pure spin mode does not couple to translations, and so one can play a (3,4) pebble game if one is interested only in translational rigidity. More generally, odd loops of double bonds do not contain this pure spin mode, while even loops of double bonds do. This complication can be addressed by keeping track of even and odd loops of double bonds. Any odd loop intersecting an even loop destroys the gearing mode, and we relabel the loop as even. If there are no even loops after looking at intersections of even and odd loops, then the original (3,3) pebble game is robust at all length scales. Near the transition where system-spanning length scales dominate, the initial version of the (3,3) pebble game is also robust as long as there is no cluster-spanning set of even loops, which is unlikely due to the intersection with odd loops such as triangles. Note also that the low-energy normal modes of rigid frictional packings show a rough equipartition between rotational and translational d.o.f. and do not contain any purely rotational modes [36].
In closing, our work opens up many new avenues for exploration in rigidity percolation with new constraint counting methods and the discovery of potentially new universality classes. It also invites us to explore not only rigid regions but floppy regions as well, which may be the key in constructing field theories of continuous rigidity percolation transitions. Finally, our new optimal rigid cluster growth algorithms do not waste material and, therefore, perhaps have a chance of being realized in living matter as well as provide mechanical examples for decision-based cluster growth that may draw links with explosive percolation [50].

ACKNOWLEDGMENTS
The authors would like to thank D. Lester for useful discussion, D. M. Sussman for helpful comments on an earlier draft, and comments from two of the referees that motivated us to restructure the presentation of the manuscript to make it more digestible. J. M. S. would like to acknowledge financial support from NSF-DMR-CMMT-1507938, and S. H. would like to acknowledge support from BBSRC Grant No. BB/N009150/1.

APPENDIX A: PEBBLE GAME
Maxwell constraint counting discussed in Sec. II assumes that every bond or constraint is an independent one. However, not every bond is an independent constraint in a random network. There may exist some redundant bonds. In order to more accurately locate the critical point where RP occurs in two-dimensional networks by keeping track of independent and redundant constraints, one can invoke the pebble game. This algorithm was described in Ref. [9] and is rooted in the following Laman condition: A two-dimensional network with N sites is minimally rigid if and only if it has 2N − 3 bonds and no subnetwork of k sites has more than 2k − 3 bonds [34]. To implement the Laman condition numerically requires checking all possible subnetworks, which is computationally expensive. The pebble game is a more computationally efficient method with a running time proportional to the number of sites times the number of bonds.
Here are a few more details of the pebble game. In a network extracted from a frictionless particle packing, since each site has 2 local d.o.f. and there are 3 global d.o.f., one plays the (2,3) pebble game. Initially, there are two pebbles on each site, then these pebbles are assigned or assigned to bonds one by one based on specific rules. The rules stem from an alternate version of the Laman condition, namely, that the bonds in the network are independent from each other if and only if for each bond in a new network formed by quadrupling each bond in the original network has no induced subnetwork of k sites and greater than 2k bonds. With this reformulation, one can check when a new bond is added to the existing set of independent bonds is itself independent via quadrupling the bond in question and invoking the Laman condition. To do this, the pebble game quadruples the new bond and tries to find a pebble covering for the 4 new bonds. If a pebble covering is not found, the new bond is not an independent constraint from the others. More specifically, the pebble game is as follows.
(1) Start with a set of covered bonds and add a new bond. (2) Look at the sites emanating from the new bond. If any of those sites has a free pebble, use it to cover the bond. Give a direction to this bond such that it points away from the site that has given up the pebble. Continue with another copy of the new bond. If the pebbles of the neighboring sites already cover existing bonds, then search for free pebbles in the directed network of existing edges. Once a free pebble is found, swap pebbles and reverse the arrows on the bonds appropriately, so that the new bond is covered. Repeat this three more times. If a free pebble is found for each of the 4 copies (the quadrupled bond), then remove three of the copies and retain one bond (with its pebble and its direction) since it is added to the existing set of independent bonds. If no free pebble is found for any of the four copies, then the new bond is not independent of the current set and it is not added to the independent set of bonds.
(3) Once all the bonds in the network have been tested, if 2N − 3 independent bonds are found, then the network is minimally rigid. If there are less than 2N − 3 independent bonds and no free pebbles, the network is overconstrained, or simply rigid, and if there are less than 2N − 3 independent bonds and free pebbles, the network is underconstrained, or floppy. To identify rigid clusters in the network, one introduces a new cluster label for an unlabeled bond and gathers three pebbles at its two incident sites. Then, three free pebbles are temporarily pinned down and the two incident sites marked as rigid. For each new nearest-neighbor site (to the two incident sites), a pebble search is performed. If a free pebble is found, the nearest-neighbor site is not mutually rigid with respect to the initial bond nor is any other site that was encountered during a pebble rearrangement; all of these sites are floppy with respect to the initial bond. However, if a free pebble is not found, the site is mutually rigid with respect to the initial bond as well as all other sites that make up the failed pebble search and so these sites are marked as rigid. Then the next-nearest neighboring sites are visited until all nearest neighbors to the set of rigid sites have been marked floppy. All bonds between pairs of sites marked as rigid are given the same cluster label. Finally, floppy and rigid marks are removed from all sites (since a site is not unique to a rigid cluster) and the process continues until there are no unlabeled bonds. In mapping out the rigid clusters, there will be two types of bonds: isostatic bonds and redundant bonds. Isostatic bonds are critical for maintaining the rigidity of the cluster, while redundant bonds can be removed without changing the overall rigidity. Only the redundant bonds can carry stress.
For the frictional case, we must incorporate the additional rotational d.o.f. for each particle into the pebble game. In addition, to account for the additional constraints due to tangential forces in the frictional case, we introduce a second bond for each frictional contact into the network. The pebble game then explores the network to see if that additional rotational d.o.f. can be independently constrained. This second bond in the network is only added to frictional contacts below the Coulomb threshold, i.e., where the normal and tangential forces are independent of each other. For contacts at the Coulomb threshold, the tangential and normal forces are no longer independent, so that only one bond in the network is needed. We, therefore, arrive at a (3,3) pebble game where contacts below the Coulomb criterion are denoted as double bonds in the network and contacts at the Coulomb threshold are denoted as single bonds in the network. Two very simple networks were discussed earlier.

APPENDIX B: GENERAL HIERARCHICAL LATTICES
Let us define a basic network motif with N B bonds in a hierarchical lattice as the general first-generation network and denote it by G 0 . Then perform the subnetwork counting, as in Fig. 11(b). Assuming that we find a n rigid subnetworks which have n bonds less than the full network motif G 0 , then, generically, the recurrence relation between two generations is p nþ1 ¼ X N B n¼0 a n p N B −n n ð1 − p n Þ n : Usually a 0 ¼ 1, and in the specific case we discussed in the main text, a 1 ¼ 8, a 2 ¼ 6, a 3 ¼ 2, and the others are zeros. The critical point is determined by the crossing point of plots of Eq. (B1) and p nþ1 ¼ p n . In Fig. 17, we can see that we need at least the first two terms of Eq. (B1) to obtain a crossing point and that they dominate the remaining terms in determining critical point p c . With the same method, we can obtain λ 1 by taking the derivative of Eq. (B1), or and use d f ν ¼ logðN B Þ= logðλ 1 Þ to find d f ν. Table I lists how p c , λ 1 , and d f ν change when we add higher-order terms to first two terms in Eq. (B1). Now let us investigate the first two terms in more detail. We have p nþ1 ¼ p N B þ a 1 p N B −1 ð1 − pÞ. To obtain a critical point, we require that To make this equation solvable in the range [0, 1], we can rewrite it as We know that p c is a number between 0 and 1, so a solution requires that a 1 > N B − 1. Since a 1 is the number of rigid subnetworks when just one bond is taken away from G 0 , we have a 1 ≤ N B . Ultimately, we obtain the equality which gives a rough criterion whether a general hierarchical lattice has a critical point, based on simple subnetwork counting.