Frictional rigidity percolation and minimal rigidity proliferation: From a new universality class to superuniversality

We introduce two new 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. For frictional rigidity percolation, we construct rigid clusters in two different lattice models using a $(3,3)$ pebble game, while taking into account contacts below and at the Coulomb threshold. The first lattice is a honeycomb lattice with next-nearest neighbors, the second, a hierarchical lattice. For both, we generally find a continuous rigidity transition. Our numerical results suggest that, for the honeycomb lattice, the exponents associated with the transition found with the frictional $(3,3)$ pebble game are distinct from those of a central-force $(2,3)$ pebble game. 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. However, the closeness of the order parameter exponent between the two cases hints at potential superuniversality. To explore this possibility, we construct a bespoke cluster generating algorithm invoking generalized Henneberg moves, dubbed minimal rigidity proliferation. The minimally rigid clusters the algorithm generates appear to be in the same universality class as connectivity percolation, suggesting superuniversality between all three types of transitions. Finally, the hierarchical lattice is analytically tractable and we find that the exponents depend both on the type of force and on the fraction of contacts at the Coulomb threshold. These combined results allow us to compare two universality classes on the same lattice via rigid clusters for the first time to highlight unifying and distinguishing concepts within the set of all possible rigidity transitions in disordered systems.


I. INTRODUCTION
At the heart of every rigidity transition is the emergence of a spanning rigid cluster-a cluster in which the bonds of an underlying contact network, be it particles or springs, are rigid with respect to each other. For disordered systems, the starting point of choice has become randomly-diluted spring networks with centralforce interactions [1][2][3][4]. For example, as bonds are randomly diluted from a triangular lattice, be it regular or one in which the lattice points have been slightly randomized, i.e. generic, the system goes from rigid with a non-zero shear modulus to floppy with zero shear modulus [5][6][7][8]. Underlying this mechanical phase change is the transition from a system with a spanning rigid cluster to a system with no spanning rigid cluster, as identified by the combinatorial (2,3) pebble game [9]. The location of this rigidity transition approximately occurs where the number of degrees of freedom are frozen out by the number of constraints-the number of force-balance equations, otherwise known as Maxwell constraint counting [10].
The nature of the rigidity transition in the centralforce, randomly bond-diluted triangular lattice is a continuous one with a correlation length exponent, ν = 1.21±0.06, an order parameter exponent β = 0.18±0.02, and a fractal dimension of the spanning rigid cluster, d f = 1.86±0.02 [6,8]. Note that while the lattice in these studies was a generic one, these exponents were obtained via the combinational pebble game and so are independent of whether or not the underlying triangular lattice is regular or generic, unlike elastic exponents or other prop-erties extracted from, say, the dynamical matrix [5,6]. From now on, we will assume any underlying lattice is regular unless otherwise specified. The exponents are slightly different from connectivity percolation transition with two-dimensional exponents ν = 4/3, β = 5/36, and d f = 91/48 [11] Despite the small difference in exponents, it was eventually argued that the rigidity percolation transition is in a different universality class from connectivity percolation since there are nonlocal interactions and rigidity is a vector problem given that forces are vectors, unlike connectivity percolation [6,8].
On the other hand, lattices with no loops, i.e. Bethe lattices, are amenable to analytical treatment and demonstrate that the spanning rigid cluster at the transition is not fractal [12,13]. 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 twodimensional case [14]. Finally, central-force models with next-neighbor springs can exhibit hybrid rigidity transitions exhibiting both continuous and discontinuous features [15]. With this rather varied set of phase transition outcomes depending just on the type of lattice, a general solution to the rigidity percolation problem is not yet clear, if it is at all possible.
Rigidity percolation with bond-bending forces has also been studied, thereby adding yet another dimension to the problem [16][17][18][19][20]. Numerical simulations of twodimensional systems measuring elastic properties suggest that bond-bending forces drive the transition into a different universality class [21]. However, since the pebble game has not yet been applied to bond-bending forces even in two-dimensions, a direct comparison between bond-bending and central-force rigidity percolation in terms of the order parameter exponent and the correlation length exponent has yet to be even made. In other words, how the nature of the rigidity transition depends on the type of force is even less understood than the type of lattice.
Particle packings also undergo a rigidity transition as a function of the packing fraction [22][23][24] and possess the additional feature of contact network rearrangements that is not allowed in a randomly-diluted spring network. The rigidity transition in such systems has been typically labeled a jamming transition. Nonetheless, the transition is from a zero to non-zero bulk modulus as the packing fraction increases, suggesting the emergence of a spanning rigid cluster [23]. This suggestion was made explicit by extracting the contact network of a two-dimensional frictionless particle packing at the jamming transition, where there are only central forces, and constructing the rigid clusters from this network via the (2,3) pebble game to find that at the onset of rigidity/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 [25].
If one were to turn to frictional systems, how does one compute rigid clusters from networks abstracted from two-dimensional frictional particle packings? A 2016 PRL by two of us outlined a new pebble game algorithm for doing so for both translational and rotational degrees of freedom [26] and applied it to molecular dynamics simulations of frictional particle packings at a fixed packing fraction experiencing slow shear. As the packing goes from jammed to unjammed and back again, rigid clusters were determined at constant strain intervals. The size of the largest rigid cluster indicated a continuous transition, and so did the observation of a roughly powerlaw distribution of rigid cluster sizes near the jamming transition [26]. Within the spanning rigid clusters, regions of floppiness were found, suggesting partial rigidity. These floppy regions are reminiscent of the rigid clusters found in the canonical bond-diluted triangular lattice. The floppy regions also appeared to be physically relevant as the non-affine motion of the particles was larger outside the rigid cluster as compared to inside. We present four such rigid cluster images close to the transition for four different values of the friction coefficient in Fig. 1.
One of the most natural questions arising from this recent study of frictional rigid clusters is whether or not the rigidity transition is actually a continuous one with regards to the onset of the spanning rigid cluster. A continuous rigidity transition would be very different from the frictionless case where there are only repulsive, central forces. Since the molecular dynamics simulation performed earlier does not allow one to tune the system to be arbitrarily close to the rigidity transition, it is difficult to assess the nature of the onset of the spanning rigid cluster using such simulations.
In this article, we construct a lattice model of a frictional packing based on a honeycomb lattice with nextnearest neighbors. We denote some bonds as below the Coulomb threshold (frictional) and some bonds as at the Coulomb threshold (sliding) and implement the frictional (3,3) pebble game to construct rigid clusters and study the nature of the rigidity transition, now in the context of frictional rigidity percolation to broaden the applicability of rigidity percolation. Finite-size-scaling provides strong evidence that the frictional percolation transition is indeed a continuous transition, with the exception of some limiting cases. The exponents that we compute numerically are distinct from those found for central-force rigidity percolation, therefore signifying a change in universality class between the two cases. We focus on rigid hinges and other one-dimensional connectors as mechanisms for propagating rigidity that are different between the two cases. However, both the order parameter exponents and the fractal dimension are statistically similar across all our countinuous transitions. This is a hint of superuniversality, which is usually considered in the context of symmetry classes in ordinary continuous phase transitions driven by symmetry breaking, even though here, there is no apparent symmetry that is explicitly broken. Our numerical results are summarized in Fig. 8.
We support our finite-size scaling analysis near the transition by algorithmically exploring the possibility of superuniversality within a subset of configurations, namely those that are minimally rigid. By (1) using concepts from invasion percolation [27], used to build only spanning clusters in connectivity percolation, and (2) extending the Henneberg moves [28], which are used to grow a minimally rigid network with central-forces only, we contruct a new algorithm to grow a minimally rigid network with frictional forces, which we dub minimal rigidity proliferation. Through this, we demonstrate that there is indeed a way to construct minimally rigid spanning clusters whose structure is the same between frictional rigidity percolation, central-force rigidity percolation, and even connectivity percolation, making the case for superuniversality, at least for such clusters.
While the honeycomb lattice approach is numerical, we also employ analytical techniques to study the frictional rigidity transition in hierarchical lattices, and compare these results with our honeycomb lattice. The hierarchical lattice solution expands the number of analytical results in the field of rigidity percolation. We conclude with a discussion about the onset of rigidity in frictional particle packings as compared to frictionless ones, investigate how our two-dimensional results connect (or do not connect) with three-dimensional results, and finally, how our results impact the field of rigidity percolation more broadly.

HC1 HC2
FIG. 2. Schematic 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 second neighbor bonds (HC2), all with probability p. Double bonds denote frictional/gear-like bonds, which occur with probability q, while single bonds denote sliding bonds.

A. Model
To motivate the model we begin with by reviewing Maxwell constraint counting in two-dimensional frictionless packings with N particles and average coordination number z [10]. The total number of degrees of freedom is 2N , while there are Nsz 2 force-balance constraints. When the number of degrees of freedom 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 [29][30][31].
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 [32]. 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 distinquish 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 degrees of freedom, there are 3(N − 1) total degrees of freedom (subtracting out the trivial global degrees of freedom in which there is no relative motion between the particles). Moreover, the interparticle forces yield z(1+q)N 2 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, then all bonds are frictional and z c = 3. If q = 1/2, then z c = 4. Therefore, Eq. 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 [33,34]. Such bounds are indeed observed in experiments [35]. Note that increasing µ increases q, and in addition, that one cannot smoothly interpolate between frictional and frictionless packings as one cannot smoothly interpolate between 2 and 3 degrees of freedom.
To construct a lattice model for frictional particle packings, we consider a honeycomb lattice with additional next-nearest neighbor (NNN) bonds. This modified honeycomb lattice has a maximum coordination number of z max = 9. We define p as the probability of bond/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 Figure II A. 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 or next-nearest-neighor, 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 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. 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 degrees of freedom and the second number denotes the number of trivial global degrees of freedom, which does depend on boundary conditions. However, Ref. [26] found that changing the number of trivial global degrees of freedom from 3 to 2 due to periodic boundary conditions did not signficantly 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. 3. A more detailed explanation can be found in Appendix A. Examples or the rigid clusters we find below, at, and above the rigidy transition for both HC1 and HC2 are shown in Figure 4. To compare frictional rigidity percolation with central force rigidity percolation, 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 two degrees of freedom. Examples of the corresponding rigid clusters are shown in Figure 5.
We now implement the common tools percolation analysis for the spanning rigid cluster to quantify the transition for HC1 and HC2 for different qs for the (3,3) pebble game and for the (2, 3) pebble game.

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. Figure 6 plots the probability that the system contains at least one spanning rigid cluster as a function of p for different system lengths L. Panel a presents data for HC1 with q = 0.5, while panel b presents data for HC1 with q = 0.7. In both subfigures, different curves with different system sizes cross near a particular value of p, which indicates the location of the transition point denoted hereafter as p c (q). In particular, p c (0.5) ≈ 0.447 (1) and p c (0.7) ≈ 0.396(1). These two critical points are very close to the results from the generalized isostaticity counting in Equation 3 with p c (0.5) = 4 9 ≈ 0.444 and p c (0.7) = 20 51 ≈ 0.392. The probability of spanning for HC1 at p c (q) for both q = 0.5 and for 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 panel c, the crossing point is more difficult to determine but estimate it to be near 0.448(1).
For the HC2 version of the model (bottom row of Fig. 6), 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 (panels d and e). Our results can be found in the first column of the Table in Figure 8. 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 panel f.

Correlation length
The correlation length ξ quantifies how two distant particles/sites interact. In percolation, it can be extracted from a two-point correlation function for the probability of occupied sites some distance from each other participating in the same rigid cluster. In a continuous transition, the correlation length diverges at the transition, while near the critical point, ξ ∼ |p − p c (q)| −ν 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 |p L c (q) − p ∞ c (q)| ∝ 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 Figure 6. The standard deviation of this distribution, ∆, yields a measurement of the correlation length exponent [36]. More precisely, Using error function fits to the data in Figure 6, we numerically differentiate the curves and fit to Gaussians to compute ∆(L) and extract the correlation length exponents ν = 1.58 ± 0.13 for HC1 with q = 0.5 and ν = 1.48 ± 0.20 for q = 0.7. Both values are within one standard deviation of each other. For the frictionless version of HC1, we obtain ν = 1.50 ± 0.07. For HC2, we find ν = 1.48 ± 0.05 for q = 0.5, ν = 1.43 ± 0.04 for q = 1.0, and ν = 1.33 ± 0.05 for the frictionless version. Figure 7e shows 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. Fig. 7a 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. 7b 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. 7e. 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 Fig. 7c-d, the order parameter exponent does not vary too much between the different models, as summarised in Fig. 8c, 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. [15]. 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 Figure 7f, we see that when q = 0.5, d f = 1.81 ± 0.06. For q = 0.7, d f = 1.80 ± 0.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 frictionles versions of the HC2 version of the model and are listed in Table 1.

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/bond, or n s [36]. At the transition, n s ∼ s −τ , where τ is the cluster size exponent. For connectivity percolation, τ > 2 strictly. How can we understand this result? The inequality τ > 2 in connectivity percolation is a consequence of the following. We start with ∞ s=1 sn s (p) + P ∞ (p) = p. Since P ∞ (p c ) = 0 for connectivity percolation, then ∞ 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 non-spanning rigid clusters, spanning rigid cluster(s), and floppy regions. There are no floppy regions in connectivity percolation. 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.
Also, when looking at how nonspanning rigid clusters grow/enlarge, the addition of a bond between two finite rigid clusters does not imply a larger rigid cluster formed by simple addition as it does in connectivity percolation, which makes it difficult to study how two nonspanning rigid clusters merge and become one. We illustrate the complexity of such a merging in Figure 9c-d where adding one double bond merges five rigid clusters and some floppy regions into a single spanning rigid cluster.
Merging through simple addition also what allows one to derive relations between the nonspanning cluster size distribution exponent and the fractal dimension of the spanning rigid cluster, for example, otherwise known as a hyperscaling relation. Given the nontrivial merging for rigidity percolation, we cannot rely on the hyperscaling relations developed for connectivity percolation. We return to the consequences of this phenomenon below.
Here, 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 τ , as above. 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 8a shows the probability for having a nonspannning 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 finite size effects. We obtain τ = 1.90±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 Fig. 8c) indicative of rather different ways the rigid clusters merage and grow in the two cases.

Summary
The results of our finite-size scaling analysis are summarized in the Table in Fig. 8. We also include the exponents for rigidity percolation using the (2,3) pebble game on the triangular lattice (denoted as T L) and for connectivity percolation (denoted as CP ) on the triangular lattice for comparison. For the frictional versions, i.e. 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 rigidity percolation on the triangular lattice obtained about twenty 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 [15], so that this case is special, just as HC1 with q = 1 is special.
Focusing now on frictional rigidity percolation versus central-force rigidity percolation, we conclude that they represent two distinct universality classes based on the rather different values of ν and τ . The order parameter exponents β are not very distinct, which could either be significant, or be due to smallness of the number which makes it difficult to discriminate between models. Setting the order parameter exponent aside for now, let us ask what mechanism 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, for both the (2,3) and (3,3) pebble games, two independently rigid clusters cannot become one rigid cluster by joining via a single bond. For a (3,3) pebble example with q = 1, consider two triangles, which are individually rigid. If they are now joined by a dou-ble bond, 18 degrees of freedom of the now 6 particles, minus 3 global degrees of freedom, are compared with 14 constraints from 7 bonds to ultimately give one floppy mode. However, for the (3,3) pebble game two independent rigid clusters connected by two double bonds makes a new rigid cluster, which is not the case for the (2,3) pebble game. Rigid hinges are another means by which two rigid clusters can merge at a point and still be rigid. In the (2,3) pebble game version, hinges consisting of single bonds between rigid clusters are always floppy, and so rigid hinges cannot exist. However, in the (3,3) pebble game, this is not true, at least for a hinge comprised of all double bonds, see Fig. 9a-b.
We conclude that two rigid double bonds and the rigid hinge (composed of double bonds) are distinct means of propagating rigidity in the (3,3) pebble game (frictional rigidity percolation) that do not occur in the (2,3) pebble game. Both motifs are more localized than possible propagators in the (2,3) pebble game and allow for rigid clusters to merge at a point or along a line, even in the absence of floppy regions. The presence of floppy regions complicates matters, see the example in Fig. 9c-d. Here, the addition of one double bond brings transmits rigidity across five nonspanning rigid clusters with connected via floppy bonds prior to the additional one double bond. Given these local motifs of connecting rigid clusters, one can understand why the sizes of nonspanning rigid clusters are typically larger for the frictional rigidity percolation as compared to central-force rigidity percolation.
And yet, even though the rigid cluster formation differs between the frictional and central-force rigidity percolation, can the differences in the order parameter exponent be that minimal? In addition to observing 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. Typically, for continuous transitions, there exist hyperscaling relations between the different exponents so that if τ is different, then one would expect d f to be different and, therefore, β to be different as well. Such relations are rigorous for connectivity percolation. For rigidity percolation, they have yet to even be derived since the problem is considerably more complex. For example, to derive the hyperscaling relation between τ and d f in connectivity percolation, there are only two types of connected clusters-those that are finite (nonspanning) and those that are spanning. In rigidity percolation, there are however three types of entitites with the floppy regions playing a role in determining how two different rigid clusters join up to form one larger rigid cluster. Moreover, for τ < 2, we expect a characteristic size cut-off. Without such hyperscaling relations, it is difficult to evaluate the significance of our findings for β and d f .

III. MINIMAL RIGIDITY PROLIFERATION
As mentioned above, the order parameter exponents and the fractal dimension are not very distinct between frictional and central force rigidity percolation. Are they in fact the same, signalling features of superuniversality for this aspect of the transition? Or is it the case that the order parameter exponents are distinct but the distinction is small, making it hard to detect? If we can find a limiting case where the order parameter exponents are actually the same, this strengthens the case for superuniversality rather than relying purely on numerical analysis. So let us now explore more explicitly connections between frictional rigidity percolation and central-force rigidity percolation via a subset of 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.
To find out if this approach is feasible, we need to construct our new algorithm step by step. Let us first review invasion percolation, which is motivated by the problem of one fluid displacing another from a random, porous medium [27]. 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 will review the Henneberg moves [28], 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 will then extend the Henneburg moves to include frictional forces and ultimately unify the two concepts, invasion percolation and Henneburg moves. The final algorithm that we introduce, minimal rigidity proliferation (MRP), 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. The above algorithm reduces to the Leath algorithm [37], 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 [38].
Let us also review the Henneberg moves for building a minimally rigid network associated with frictionless particles, i.e. for central forces only. A minimally rigid graph in this case is also known as a Laman graph. Minimal rigidity occurs when the degrees of freedom 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 Figure 10 (top): • add one site and two bonds between this site and two points in G , then G (N + 1, N B + 2) is the new minimally rigid network (Type I move).
• 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 two degrees of freedom 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 Figure 10, 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 Figure 10, 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 will first focus on (2, 3) minimal rigidity and then address (3, 3) minimal rigidity. 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. 11 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 next-nearest 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, 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.

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 rescaling of time in which two bonds are added in one time step as opposed to two bonds in two time steps. Moreover, such a local modification from invasion percolation shouId not alter the long wavelength behavior which captures the universality class. In Fig. 12 we show an example of a spanning minimally rigid cluster on the honeycomb lattice with NNN bonds using minimal rigidity proliferation. We also 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 Figure 13a, 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. Figure 13b shows P ∞ , the fraction of the system in the spanning cluster, converging to zero when system becomes infinitely large, suggesting a continuous transition. Networks constructed in this way correspond to the frictionless 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 gowth 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, we expect the same configurations as above, just with half single and half double bonds. Both central-force percolation and frictional rigidity percolation collapse to an identical construction in this case. We have already argued that adding two bonds at a time involves a rescaling of time from adding one bond at a time for central-force percolation, so that we also expect the q = 1/2 frictional process to be in the same universality class as connectivity percolation. In other words, within this subset of minimally rigid configurations, we expect superuniversality to emerge: All three universality classes collapse into one! 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 have some understanding as to why some exponents appear to be quite close in value.

IV. HIERARCHICAL LATTICES
While we have presented predominantly numerical results so far, one exactly solvable rigidity percolation model is rigidity percolation on hierarchical lattices. However, such lattices have the particular property that the exponents depend on details of the lattice. To determine if this property prevails in the frictional (3, 3) pebble game, we will first review prior results using the (2, 3) pebble game with central forces only and then generalize.

A. Review: central forces only
It has been previously shown that central force rigidity percolation transitions in such lattices are a continuous transition [39,40]. To understand this finding, let us start with the generation of a particular hierarchical lattice known as the Berker lattice [40]. Given two points and a bond as in Figure 14a, 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 n th 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 Every type of subnetwork is a way to obtain a spanning rigid network between two ends(black dots) and its probability is calculated, as well as the probability for an occupied bond to be in the spanning rigid cluster. (c) First four pn as function of p. pn tends to converge to a step function at pc = 0.9446 which jumps from 0 to 1.
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 Figure 14b. 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. 14c, 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.
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. 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 type of subgraphs containing a spanning rigid cluster that were not allowed in the central force case, as shown in Fig. 15a. 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 = λ 2 P R,n , such that λ 2 (q = 1) = 0.3511. Since λ 2 (q = 1) < 1, the rigidity transition is continuous with d f ν = 3.181. Now we consider q < 1. After keeping track of what subgraphs are rigid between the two black circles in Figs. 14b and 15a in the presence of both double and single bonds, we obtain n (84q 7 − 210q 6 + 96q 5 ) + 16p 6 n q 6 (11) + p 5 n (−8q 5 + 10q 4 ). 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. 15b. The fact that the correlation exponent depends continuously on q is not necessarily unique as has been found in Ising models on hierarchical lattices, for example [41]. This sensitivity is presumably due to the special nature of the hierarchical lattice, as detailed in Appendix B.

V. 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, which are representative of contacts that are below the Coulomb threshold. In doing so, we have uncovered a new universality class in the realm of rigidity percolation, namely that of frictional rigidity percolation, which is to be directly compared with central-force rigidity percolation on the same lattice. The apparent sensitivity of the nature of the rigidity transition to the type of lattice and type of force has so far made such a direct comparison between two universality classes not possible. For instance, a direct comparison of central-force rigidity percolation to bond-bending rigidity percolation in terms of the order parameter exponent, etc., has not been done because the pebble game has not yet been applied to bond-bending forces. By expanding the scope of rigidity percolation, the direct comparison presented here should help to formulate a more general framework for rigidity transitions.
More specifically, we make a direct comparison between central force rigidity percolation and frictional rigidity percolation on honeycomb lattices with additional next-nearest bonds. We find different correlation length and rigid cluster size distribution exponents between the two cases but a statistically similar order parameter exponent but could ultimately be demonstrated to be distinct with larger system size studies. Given the different correlation length and rigid cluster size distribution exponents, 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 rigidity percolation that are distinct from central-force rigidity percolation. The rigid hinge is a zero-dimensional contact between two-dimensional rigid clusters, while the double bond is a one-dimensional contact. In centralforce rigidity percolation, there must be additional supporting bonds to connect up rigid clusters. The lowerdimensional rigid cluster connection mechanisms in frictional rigidity percolation (as compared to central-force rigidity percolation) perhaps drive the distiction between universality classes. For the hierarchical lattice, not only are there two different universality classes for each respective model, the exponents depend continuously on the fraction of the double bonds in the frictional case.
It is interesting to note that in three-dimensions that central-force percolation is thought to exhibit a firstorder rigidity transition [14], though Laman's theorem (as implemented via the pebble game) is not rigorous in three-dimensions. Perhaps the mechanisms for connecting up rigid clusters by small numbers of additional supporting bonds are few and far between, thereby demanding a first-order transition, or one needs to simulate much larger systems to observe any fractal nature of the rigid spanning cluster. In addition to higher dimensions, one could also ask about how rigid clusters connect up in the general (k, l) pebble game with and without double bonds. Presumably other mechanisms for connecting rigid clusters emerge to warrant additional universality classes not yet discovered. Or, possibly, all k > 3 collapse to frictional rigidity percolation, at least in two-dimensions.
Since the difference between the order parameter exponent is small between central-force and frictional rigidity percolation, we combined Henneberg moves (extended here to the (3, 3) pebble game) and invasion percolation to construct another new model, minimal rigidity proliferation (MRP), that can be implemented for both central-force and frictional rigidity percolation. The rigid cluster in MRP grows in an additive fashion, unlike rigidity percolation where rigid clusters surrounded by floppy regions may not always connect up in a way that leads to larger rigid clusters. With minimal rigidity proliferation, there are no redundant bonds-the spanning rigid cluster is built in a "clever" way, which is to be contrasted with the tuning by pruning approaches [42], where springs are removed in a way that results in the least change in the bulk modulus, for example, and the jamming graph approach [43] where minimally rigid clusters with the geometric constraint of local mechanical stability are generated. In minimal rigidity proliferation, the order parameter exponent is the same across connectivity, central-force rigidity percolation, and frictional rigidity percolation. This would be the first time superuniversality is observed in rigidity percolation in a way that goes beyond transfer matrix methods. Our work also suggests that looking at minimally rigid configurations-a subset of all possible configurations within rigidity percolation-represents a paradigm shift in the way phase transitions are viewed 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 rigidity percolation 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 indicate a continuous transition [44]. In frictionless packings, there are no redundant bonds, which makes the constraint counting rather straightforward. However in frictional packings, the redundant bonds emerge such that the constraint counting is more intricate and, therefore, perhaps more non-mean-field. It will be interesting to apply our constraint counting algorithm to experimental frictional particle packings to test the applicability of our approach as well as to compare the rigid clusters with dynamical matrix calculations, for example.
Finally, we are currently exploring a limitation of the (3, 3) pebble game [45]. 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 Figure 3 without the diagonal bond). From the (3,3) pebble game this configuration is floppy with one floppy mode that is a pure spin mode and the particles spinning as gears. 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 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 remedied by keeping track of even and odd loops of double bonds. Any odd loop intersecting an even loop drives the even loop to become odd, if you will. 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 lengthscales. Near the transition where systemspanning lengthscales 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 degrees of freedom and do not contain any purely rotational modes [34].
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 [46].
The authors would like to thank D. Lester for useful discussion and D. M. Sussman for helpful comments on an earlier draft. JMS would like to acknowledge finanical support from NSF-DMR-CMMT-1507938 and SH would like to acknowledge support from BBSRC grant BB/N009150/1.

Appendix A: The pebble game
The constraint counting in Section II assumes that every bond/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 rigidity percolation 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 [32]. To implement the Laman condition numerically requires checking all possible subnetworks, which is comptuation-ally 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 degrees of freedom and there are 3 global degrees of freedom, one plays the (2,3) pebble game. Initially, there are two pebbles on each site, then these pebbles are assigned/covered 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, the network formed by quadrupling the bond has no induced subnetwork of k sites and greater than 2k − 3 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 re-arrangement, all 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 changed the overall rigidity. Only the redudant bonds can carry stress.
For the frictional case, we must incorporate the additional rotational degree of freedom 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 degree of freedom 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, like in Fig. 14b, and assume 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 Usually a 0 = 1, and in the specific case we discussed before, 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. 16, we can see that we need at least the first two terms of Eq. B1 to obtain a crossing point, and that they contribute much more than other terms in determining critical point p c . By same method, we can obtain λ 1 by taking derivative Now let's 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 p c that 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.