# The statistics of frictional families

###### Abstract

We develop a theoretical description for mechanically stable frictional packings in terms of the difference between the total number of contacts required for isostatic packings of frictionless disks and the number of contacts in frictional packings, . The saddle order represents the number of unconstrained degrees of freedom that a static packing would possess if friction were removed. Using a novel numerical method that allows us to enumerate disk packings for each , we show that the probability to obtain a packing with saddle order at a given static friction coefficient , , can be expressed as a power-series in . Using this form for , we quantitatively describe the dependence of the average contact number on friction coefficient for static disk packings obtained from direct simulations of the Cundall-Strack model for all and .

###### pacs:

83.80.Fg,45.70.-n,81.05.RmGranular media are fascinating, complex materials that display gas-, liquid-, and solid-like behavior depending on the boundary and driving conditions. Frictional forces are crucial for determining the structural and mechanical properties of granular media in the solid-like state makse_compress (). For example, friction plays an important role in setting the angle of repose angle (), determining the width of shear bands in response to applied stress shear (); xu (), and enabling arches to form and jam hopper flows pak ().

For static packings of frictionless spherical particles, it is well known that the minimum contact number required for mechanical stability gao () is , where is number of contacts among particles in the force-bearing backbone of the system and is the spatial dimension. However, at nonzero static friction coefficient , fewer contacts are required for mechanical stability with and in the large- and limits.

Several computational studies have measured the contact number as a function of for packings of frictional disks and spheres using ‘fast’ compression algorithms that generate amorphous configurations silbert1 (); makse (); makse1 (). In particular, these studies find and in the and limits, respectively, for bidisperse disks silbert (); papanikolaou (). For intermediate values of , smoothly varies between and . However, it is not currently known what determines the contact number distribution for each and form of for a given packing preparation protocol. The ability to predict the functional form of the contact number with is important because controls the mechanical somfai () and vibrational bertrand () properties of granular packings.

In this Letter, we develop a theoretical description for static frictional packings at jamming onset in terms of their ‘saddle order,’ or the number of contacts that are missing relative to the isostatic value in the zero-friction limit, . In contrast, previous studies used packings as the reference henkes (). Using a novel numerical procedure (the ‘spring network’ method) that allows us to enumerate packings for each and molecular dynamics (MD) simulations of the Cundall-Strack model cundall () for frictional disks, we show that characterizes the dimension of configuration space that static packings occupy. Frictional packings with contacts form one-dimensional lines in configuration space, packings with populate two-dimensional areas in configuration space, and packings with larger form correspondingly higher-dimensional structures in configuration space. We assume that the probability for obtaining a static packing with saddle order at a given , , is proportional to the volume occupied by force- and torque-balanced th order saddle packings in configuration space. We find that can be written as a power-series in , , where are the normalized coefficients of the power series. Using this form, we are able to quantitatively describe the dependence of the average contact number on the friction coefficient for static disk packings obtained from MD simulations of the Cundall-Strack model over a wide range of and in the large system limit.

We generated static packings of bidisperse (- mixtures of particles with equal mass and diameter ratio ) frictional disks in square cells with periodic boundary conditions using two methods. First, we implemented a packing-generation algorithm in which the system is isotropically compressed or decompressed (followed by energy minimization) to jamming onset gao () at packing fraction that depends on the saddle order. Pairs of overlapping disks and interact via repulsive linear spring forces in the direction of the center-to-center separation vector . We implemented the Cundall-Strack model for the frictional interactions. When disks and come into contact, a tangential spring is initiated with a force that is proportional to the tangential (perpendicular to ) displacement between the disks. The tangential displacement is truncated so that the Coulomb threshold, , is always satisfied. When the disk pairs come out of contact, we set to zero.

We characterize each disk packing in configuration space by plotting the second invariant of the distance matrix, versus , where and are the - and -coordinates of particles and . (Note that is invariant to uniform translations and rotations, as well as particle-label permutations, of the system.) The plot of versus in Fig. 1 (a) for packings with and illustrates several important features. First, in the limit, packings occur as distinct points in configuration space (or versus ) gao (). Second, as increases, packings form one-dimensional lines in configuration space that emanate from packings. The packings that are stabilized at low are displaced in configuration space from the packings as highlighted in Fig. 1 (b). In contrast, the packings that occur at large approach the packings. Thus, we find that the lengths of the lines increase with . (Fig. 1 (c)) and higher-order saddle packings populate areas and higher-order volumes in configuration space.

We also developed a numerical technique (‘spring network’ method) to enumerate packings at each . The method is best explained using an example. In Fig. 2 (a), we show an packing of frictionless disks with contacts, which corresponds to the packing circled in the lower right corner of the plane in Fig. 1 (a) and (b). To systematically generate packings with contacts, we break one of the contacts in this packing (e.g. the contact between disks and in Fig. 2 (a)) and constrain its separation to be , while the other contacts are constrained to be . With these constraints and as a function of , we implement the successive compression and decompression packing-generation algorithm gao () to find packings at jamming onset, . This procedure is repeated for each of the other contacts in the packing in Fig. 2 (a) to yield the , branches in Fig. 1 (b), and then for each of the packings. As shown in Fig. 1 (b), we find overlap between the branches from the spring network method and the packings generated from simulations of the Cundall-Strack model. and higher-order saddle packings can be obtained using a similar procedure, except multiple contacts are broken, as shown in Fig. 1 (c). Thus, a family emanates from each packing with branches with dimension in configuration space for each saddle order.

We are interested in determining the probability to obtain an th order saddle packing at a given averaged over families. We will assume that is proportional to the volume in configuration space of the th order saddles that can be force- and torque-balanced with friction coefficient or less:

(1) |

where is a lengthscale required to make Eq. 1 dimensionally correct. To measure , we first employ the spring network method to generate a grid of points for each branch of saddles of order . At each grid point characterized by , we determine the minimum friction coefficient required to achieve mechanical equilibrium for that configuration, using Monte-Carlo moves to search the null-space of the force- and torque-balance matrix shaebani (). The allowed configuration volume is determined by integrating over the contour (as shown in Fig. 2 (c)) such that for a given th order branch. The total admissible volume in configuration space is obtained by summing the volumes over all th order branches. We show in Fig. 3 that scales linearly with for and .

Thus, from the results in Fig. 3, we postulate the following form for the normalized probability for an th order saddle:

(2) |

where and the highest order saddle is in 2D. The coefficients , where is the number of packings for a given and the normalized coefficients . A reasonable estimate for the number of branches stemming from each packing is the permutations of , . We show in Fig. 4 (a) that Eq. 2 with and yields qualitatively correct results for the measured probabilities to obtain a given th order saddle from MD simulations of the Cundall-Strack model for (Fig. 4 (b)). For example, zeroth order saddles are most highly probable for small , and the highest order saddles are most probable for . However, as shown in Fig. 4 (b), we obtain a much better fit to the data from the Cundall-Strack model using

(3) |

which indicates an excess of th order saddles for small that likely originates from rattler particles that join the force-bearing network during compression.

The average contact number can be obtained by calculating . Our strategy is to use the results for small systems (i.e. ) to predict and versus for large . We show in Fig. 4 (c) that the predictions from Eqs. 2 and 3 agree quantitatively with with no additional parameters for and obtained from the Cundall-Strack model. Thus, we have developed a method to calculate for large systems by enumerating frictional families in small systems.

The calculations presented in this Letter provide a framework for addressing several important open questions related to frictional packings. For example, why does the crossover from the low- to high-friction limits in the average contact number and packing fraction occur near for disks papanikolaou (); silbert () compared to for spheres silbert (); makse () for fast compression algorithms. In addition, using the methods described above, we will be able to determine how the crossover from low- to high-friction behavior depends on the compression rate and degree of thermalization. Such calculations are crucial for developing the ability to design granular assemblies with prescribed structural and mechanical properties.

We acknowledge support from NSF Grant No. CBET-0968013 (MS) and DTRA Grant No. 1-10-1-0021 (CO and SP). This work also benefited from the facilities and staff of the Yale University Faculty of Arts and Sciences High Performance Computing Center and NSF Grant No. CNS-0821132 that partially funded acquisition of the computational facilities.

## References

- (1) H. A. Makse, D. L. Johnson, and L. M. Schwartz, Phys. Rev. Lett. 84 (2000) 4160.
- (2) L. Vanel, D. Howell, D. Clark, R. P. Behringer, and E. Clément, Phys. Rev. E 60 (1999) 5040(R).
- (3) D. Fenistein, J. W. van de Meent, and M. van Hecke, Phys. Rev. Lett. 92 (2004) 094301.
- (4) N. Xu, C. S. O’Hern, and L. Kondic, Phys. Rev. E 72 (2005) 041504.
- (5) K. To, P.-Y. Lai, and H. K. Pak, Phys. Rev. Lett. 86 (2001) 71.
- (6) G.-J. Gao, J. Blawzdziewicz, C. S. O’Hern, and M. D. Shattuck, Phys. Rev. E 80 (2009) 061304.
- (7) L. E. Silbert, D. Ertas, G. S. Grest, T. C. Halsey, and D. Levine, Phys. Rev. E 65 (2002) 031304.
- (8) P. Wang, C. Song, C. Briscoe, K. Wang, and H. A. Makse, Physica A 389 (2010) 3972.
- (9) C. Song, P. Wang, and H. A. Makse, Nature 453 (2008) 629.
- (10) L. E. Silbert, Soft Matter 6 (2010) 2918.
- (11) S. Papanikolaou, C. S. O’Hern, and M. D. Shattuck, Phys. Rev. Lett. 110 (2013) 198002.
- (12) E. Somfai, M. van Hecke, W. G. Ellenbroek, K. Shundyak, and W. van Saarloos, Phys. Rev. E 75 (2007) 020301(R).
- (13) T. Bertrand, C. F. Schreck, C. S. O’Hern, and M. D. Shattuck, “Vibrations in jammed solids: Beyond linear response,” (unpublished) (2013); xxx.lanl.gov/pdf/1307.0440.pdf.
- (14) S. Henkes, M. van Hecke, and W. van Saarloos, Europhys. Lett. 90 (2010) 14003.
- (15) P. A. Cundall and O. Strack, Geotechnique 29 (1979) 47.
- (16) M. R. Shaebani, T. Unger, and J. Kertész, Phys. Rev. E 79 (2009) 052302.