Topological and geometric measurements of force chain structure

Developing quantitative methods for characterizing structural properties of force chains in densely packed granular media is an important step toward understanding or predicting large-scale physical properties of a packing. A promising framework in which to develop such methods is network science, which can be used to translate particle locations and force contacts to a graph in which particles are represented by nodes and forces between particles are represented by weighted edges. Applying network-based community-detection techniques to extract force chains opens the door to developing statistics of force chain structure, with the goal of identifying shape differences across packings, and providing a foundation on which to build predictions of bulk material properties from mesoscale network features. Here, we discuss a trio of related but fundamentally distinct measurements of mesoscale structure of force chains in arbitrary 2D packings, including a novel statistic derived using tools from algebraic topology, which together provide a tool set for the analysis of force chain architecture. We demonstrate the utility of this tool set by detecting variations in force chain architecture with pressure. Collectively, these techniques can be generalized to 3D packings, and to the assessment of continuous deformations of packings under stress or strain.


I. INTRODUCTION
Densely packed granular materials exhibit a rich internal network of physical interactions [2][3][4][5][6][7], which have come to be referred to as force chains (FIG. 1). The structure of these chains is thought to play a critical role in the response of the media to perturbations including acoustic propagation [8][9][10][11] and shear [6]. However, the exact physical mechanisms linking mesoscale network organization to global scale bulk properties are not well understood. Developing measurements that effectively isolate features of force chains responsible for the bulk properties of packings that induce desirable physical behaviors is a vital first step toward designing systems that exhibit those behaviors.
However, there is not currently a widely-accepted definition of what constitutes a force chain. Here, we rely on * cgiusti@seas.upenn.edu † dsb@seas.upenn.edu the data-driven perspective introduced in [1], which offers a fundamental mapping between measurements of granular media and mathematical objects known as weighted graphs [12]. Specifically, constituent particles are represented as network nodes and the normal forces between pairs of particles in contact with one another are represented as network edges whose weight is proportional to the value of the normal force (FIG. 1c). Several studies of granular media have examined properties of similarly constructed networks. For example, contact loops and other topological measures have been examined in a number of systems [13][14][15][16][17][18]. Granular force networks have also been analyzed using tools from persistent homology [19][20][21][22][23][24] as well as with network-science based measures [1,[25][26][27][28][29]. In this study, we express the network architecture as an adjacency matrix A whose elements A ij encode the normal force between node i and node j. Using community detection techniques [30,31] informed by a geographically-constrained null model [32], we extract subgraphs with branch-like characteristics reminiscent of force chains. The ability to identify force chains in a data-driven manner unearths the more fundamental problem of identifying characteristics of force chains that drive or predict bulk properties. One simple approach lies in describing the statistics of individual interactions (force-weighted contacts) between particles. Yet, these statistics do not accurately predict signal transmission through the material, likely due to their naivety with regards to mesoscale architecture or collective dynamics [11]. Alternative options include statistics that directly describe the geometry or topology of the force chains, therefore capturing mesoscale architectural properties of the material [1,27].
In this paper, we discuss three such measures and assess their ability to (i) identify force chain structure in packings of granular particles, and (ii) detect variations in force chain architecture as a function of the applied pressure. This analysis is carried out on both laboratory packings and on simulated packings in two dimensions. We begin with the previously defined "topophysical" statistic known as the gap factor [1], which utilizes a blend of topological and physical information to characterize force chains. Teasing apart the complex relationship between geometry and topology in the definition of the gap factor motivates us to consider a pair of statistics which are respectively purely geometric and topological in flavor: the hull ratio, first described in [33], which measures the density of the packing around each chain; and a novel topological statistic of a network, the topological compactness factor (TCF), which measures physically relevant structure in the topology of the contact network underlying the force chains. The TCF is sensitive to mesoscale features like branching in chains or compact, highly-interconnected portions of the force network, where the system's local stability can be respectively weakened or reinforced [34].
We find that the gap factor and hull ratio are better able to extract force chain structure at each pressure than the purely topological statistic. On the other hand, the TCF is more sensitive to pressure than either of the other two measures. Further, its underlying mathematics can be extended naturally to quantify continuous deformations of packings in variable environmental conditions, thus forming an important tool for examining heterogeneous microstructures in particulate matter. The three statistics together provide a broad picture of the mesoscale structure of the packing, and our findings suggest that both physical and topological information is useful to consider when studying granular systems.

II. METHODS
To illustrate the effectiveness of our statistics, we examine two case studies: (i) granular experiments where friction and gravity play a role, and (ii) frictionless and gravity-less simulations with periodic boundary conditions. In both cases, approximately 600 bidisperse disks with a ratio of 50:50 are confined in two dimensions and interact with each other in a Hertzian-like manner. In Brightness indicates the strength of forces between particles. (c) The corresponding weighted graph overlaid on the photoelastic image, where the thickness of the line segments (edges) are proportional to the normal force between the two particles they connect (nodes). (d) Community structure extracted from the force-weighted network as a function of the resolution parameter γ in the modularity quality function. Colors indicate distinct communities, and warmth corresponds to network force, the amount by which the total inter-particle forces exceed γ times the mean inter-particle force.
each packing, which is jammed under constant pressure, we measure force-contact networks for approximately the same 7 different values of confining pressure.

A. Granular experiments
We perform experiments on a vertical 2D granular system of bidisperse disks that are confined between two sheets of Plexiglas. The particles are 6.35 mm thick, their diameters are d 1 = 9 mm and d 2 = 11 mm (which yields a diameter ratio of approximately 1.22), and they are cut from Vishay PSM-4 photoelastic material to provide measurements of the internal forces. These particles have an elastic modulus of E = 4 MPa. We produce new configurations by rearranging the particles by hand, and we increase the pressure by placing additional brass weights on the top surface of the packing. The values of pressure, which we report in units of the elastic modulus E (recall that the configuration is twodimensional), are 2.7 × 10 −4 E, 4.1 × 10 −4 E, 6.7 × 10 −4 E, 1.1 × 10 −3 E, 2.2 × 10 −3 E, 3.8 × 10 −3 E, and 5.9 × 10 −3 E. See Refs. [11,27,35] for additional details about the ex-periments.
For each of 21 particle configurations and the 7 values of pressure, we compute particle positions and forces using two high-resolution pictures of the system (FIG.  1a,b). We use one image, which we take without polarizers, to determine the particle positions and contacts. (See [36] for a description of the technique.) We take particles to be in contact if they have a measurable-magnitude force based on our photoelastic calculations. Using a second image that we take with polarizers, we then determine the particle contact forces by solving the inverse photoelastic problem [36].

B. Frictionless simulations
We perform numerical simulations of bidisperse frictionless disks with a diameter ratio of 1.22. Particles interact via a Hertzian potential in a box with periodic boundary conditions in both directions and zero gravity [37][38][39]. This model has been well-studied and it is significantly different from our experimental system with friction, gravity, and fixed boundaries. We generate mechanically-stable packings via a standard conjugategradient method [37]. We then perform simulations for a fixed packing fraction and volume, and we analyze 20 mechanically-stable packings at each packing fraction φ. We choose the seven values of the packing fraction so that the mean pressure p at that packing fraction [40]  where the modulus E is defined as the energy scale for the Hertzian interaction divided by the mean Voronoï area of a particle in the packing. The lowest value of φ provides a data point for a jammed packing that is less dense than what is accessible in our experiments.

III. FORCE CHAIN EXTRACTION VIA COMMUNITY DETECTION
Following [1], we represent the forces between particles using a weighted graph and we apply community detection techniques to extract putative force chains from the resulting network. First, we construct a simple, unweighted contact graph B which has as nodes the particles in the system and which has an edge between nodes precisely when the corresponding particles are in contact. We extend this to a weighted force graph W by including edge weights W ij equal to the normal force between the particle i and j.
To extract force chains from the force network W, we need to find sets of particles with strong interparticle forces among themselves. We accomplish this by using community detection techniques [30,31] to partition the vertices into groups (or communities) in a fashion that maximizes a modularity quality function [41]. The particular function we utilize is given by where node i is assigned to community g i , node j is assigned to community g j , δ is the Kronecker delta function, and P ij is a null model. The parameter γ is the structural resolution parameter that can be used to tune the number of communities detected: small values of γ will produce fewer larger communities, while larger values of γ will produce many smaller communities. We optimize modularity for several values of the resolution parameter between γ = 0.1 and γ = 2.1 in steps of 0.2. As in [1,32], we employ a geographical null model given by where ρ is the mean edge weight in the network, to reflect the constraints on network structure induced by the twodimensional packing. Maximizing the modularity quality function partitions the particles in the system into network communities, and by tuning the resolution parameter appropriately, we can identify groups of particles that are visually reminiscent of force chains (FIG.  1d). As the modularity function is non-convex, we use a locally greedy ("Louvain"-like [42]) algorithm [43] to obtain approximations of the optimal partition of nodes into communities. We report results across twenty such optimizations of the modularity quality function per packing for the experimental system, and across one-hundred optimizations per packing for the simulations.

IV. CHARACTERIZING FORCE CHAINS
As we vary the resolution parameter γ, the character of the force chains we extract from a packing change dramatically (FIG. 1d). In order to provide consistent characterizations of structure, we wish to select an "optimal" resolution at which the results of the community detection strongly resemble our heuristic notion of what a force chain should be: heterogeneous and branch-like. Defining statistics to quantify such a branching architecture is an open area of investigation [1]. In what follows, we define three different statistics that we use to characterize the sets of particles found from community detection. We then discuss the ability of these statistics to identify an optimal resolution parameter for force chain extraction, and consider their sensitivity to the pressure applied on the system.
A. The gap factor: a previously defined "topophysical" statistic Before discussing new statistics, we recall the gapfactor, which measures the correlation between physical (a) The gap factor of a force chain measures the discrepancy between the topology of the contact network of a force chain and its physical layout. Force chains for which the two are poorly correlated display branching or curving behavior. (b) The hull ratio of a force chain is the ratio of the area of the constituent particles to those of its convex hull. Curved or branched chains (top) have low hull ratio, while straight or clustered chains (bottom) have a high hull ratio. (c) The topological compactness factor (TCF) measures mesoscale connectivity properties of the contact graph that contain physically relevant information. (left) Binary contact graph for a force chain. (middle) A compact region (grey triangle) is a cluster of particles in the chain which are connected to all of their neighbors in a triangular grid, while branch points (red) are particles with more than two neighbors, none of which share an edge. (right) After collapsing the leaves (green) in the contact graph to a single vertex, we can enumerate fundamental cycles (purple, yellow, blue, orange) in the graph to measure the prevalence of compact regions versus branch points of the chain. distance (measured using the standard Euclidean metric) and "hop distance" (which counts distance measured only along network edges). In this way, the gap factor combines the physical and topological information about the system into a "topophysical" estimate of the branch-like structure of the force chains. Force chains with gaps, branches, and rings have a larger hop distance than physical distance, and the presence of such complicated shapes decreases the correlation between the physical and hop distance (FIG. 2a).
The physical distance between nodes in a force chain is the usual Euclidean distance matrix, which we denote L p .
The hop distance matrix L t of a community is obtained by restricting the binary contact network B to nodes assigned to the community c to obtain a submatrix B c . The (i, j)-entry of L t is then given by the minimum number of edges traversed in any path through B c from node i to node j.
To obtain a global statistic, we weight each community by its size so that larger communities are weighted more heavily than small communities to reflect their relative influence on the system's behavior. We define the gap factor g c of a community c to be where r c is the value of the Pearson correlation coefficient between the strictly upper triangular entries of L t and those of L p , and s max is the size of the largest community. We define the systemic gap factor as where the quantity n >1 is the total number of communities, excluding singletons.

B. The hull ratio
To complement the topophysical gap factor, it is useful to consider a purely geometric quantity that depends only on the relative locations of particles in physical space, and therefore does not use notions of network topology. A simple measure of how tightly packed a group of particles is obtained by taking the ratio of the total area of the particles to the area of the convex hull of the particles (FIG. 2b). Following [33], we define the hull ratio of a community c to be With this definition, R(c) will be close to unity for more compact and linear structures, whereas groups of particles with heterogeneous force chain geometry will have smaller values of R(c). To obtain a summary statistic for the structure of a packing as a whole, we define the weighted hull ratio given by the average of the hull ratios for all the communities, weighted by the number of particles in the community (again excluding singletons).

C. The topological compactness factor
To complement the previously defined gap factor (a "topophysical" statistic) and hull ratio (a physical statistic), we turn now to defining a purely topological statistic that is sensitive to the tendency of the force chain to contain compact regions or to branch (FIG. 2c).
To motivate the development of this new statistic, we first note that the contact structure of chains is strongly reflective of the organization of forces within the system, and that the mesoscale properties of the packing near these compact regions or branch points is of particular structural significance. A compact cluster of particles is one in which the constituent particles exert pairwise forces on one another, forming a roughly triangluar tiling in the 2D system which necessarily constrains its motion. A branch point occurs when the chain diverges into two or more locally distinct paths across a single particle (FIG.  2c). There is theoretical work suggesting that compact regions in the force network increase rotational stability of the system [44], and we posit that branch points increase the propensity for the chain structure to slip under perturbation of the system. Here, we propose a statistic to explicitly assess these two tendencies: the topological compactness factor (TCF).
The TCF is conceptually related to the graph-theoretic clustering coefficient, which measures the percentage of neighbors of each vertex that share an edge, however the TCF focuses attention on the mesoscale structure of each chain rather than its behavior at particular vertices. The TCF derives from techniques drawn from algebraic topology, and is in fact a summary statistic of a much more powerful measurement of the structure of the chain in terms of the homology of its clique complex. However, as it is possible to describe the numerical TCF in purely graph-theoretic terms in the case of a static 2D packing, such as those considered in this paper, we instead provide such a description here and leave the more complete definition to the Appendix A.
While compact regions can be relatively easily ex-tracted from the contact graph by the triangular connectivity patterns they induce, counting branches is an inherently non-local process. Some branches are most easily identified by counting the number of leaves of the force chain graph (FIG 2c (middle)), but other branches may combine to form loops which are entirely internal to the chain. In order to put branch points into the same framework as compact regions, we begin by modifying the graph so that both can be expressed in the language of "cycles". Specifically, let c be a community in the packing, thought of as its underlying contact graph B c .
Recall that a leaf in B c is any vertex that shares an edge with precisely one other vertex (FIG. 2c (middle, green)). Construct a new graph B c from B c in two steps: first, introduce a new vertex π and add an edge from vertex i to π precisely when there is some leaf of B c whose unique edge is to vertex i; second, delete all leaves from B c (FIG. 2c (right)).
Recall that a cycle in a graph G is a path that begins and ends at the same vertex (FIG. 2c (right)). A common difficulty when enumerating cycles in graphs lies in determining which cycles to count, as illustrated by the common "How many rectangles do you see?" puzzle. Here, we wish only to consider a basis for the set of cycles under the operation of composition, and so use the matroid-theoretic notion of fundamental cycles. Recall that a forest is a graph without cycles and that a spanning forest for G is a forest obtained by removing edges from G until no cycles remain, without changing the number of connected components. The number of fundamental cycles in G, written F (G), is defined to be the number of edges which must be removed from G to create a spanning forest [45].
Observe that cycles of length 3, also known as 3cliques, are present precisely where three particles in the packing are all in pairwise contact (FIG. 2c (right, blue)). Such a structure will be our indication that the cycle is induced by membership in a compact region, rather than a branch point; denote by T (G) the number of such cycles in G.
Finally, we define the topological compactness factor of a force chain c to be and set the TCF to be zero when F ( B c ) = 0. Chains for which the TCF is large will tend to consist of compact, highly interconnected regions of particles, while those for which the TCF is small will tend to be long, branching chains. In order to compare packings, we define a summary statistic (following [1]), which we refer to as the systemic topological compactness factor where n is the number of communities {c}, s c is the size of community c and s max is the size of the largest community.

D. Comparison of community measures
It is interesting to compare the previously defined gap factor to the geometric hull ratio and topological compactness factor. In FIG. 3, we show scatter plots of (a) the hull ratio vs. the gap factor and (b) the compactness factor vs. the gap factor. In both plots, points correspond to the mean of each measure across optimizations of the modularity quality function for each packing, pressure, and choice of the resolution parameter between 0.1 and 2.1. Though the hull ratio and gap factor do not appear linearly correlated, it is evident from panel (a) that there is a relationship between their values and we find that they are indeed strongly Spearman correlated (ρ = −0.72). On the other hand, the TCF is only weakly correlated with the gap factor (ρ = 0.37), and this topological statistic thus likely captures different information about community structure than either the gap factor or the hull ratio measures.
The scatter plot of the hull ratio against the gap factor contains two additional features that stand out by eye: the vertical bands at g s ≈ 0 and g s ≈ 0.5 and the diagonal striations near the top of the figure. Investigation of the vertical stripes shows that they occur at small values of the resolution parameter. At low γ, the hull ratio is sensitive to the packing pressure, which can be seen by the vertical stratification of colors in both bands. Moreover, the spread of the vertical stripes along the horizontal axis can be understood from the definition of the weighted gap factor in EQN. 4. At low resolutions, there will either be a single large, compact community or one large, compact community with a few very small communities dispersed throughout the packing. The exact number of communities detected at a given pressure and resolution parameter will in part be dependent on the structure of the particle configuration. If there is a single large community that contains the majority of the particles, then r c ≈ 1 and the gap factor will be close to zero (the first vertical band). However, the gap factor is highly sensitive to the number and size of communities; if a given configuration contains one large community and a just a couple of small communities of size two or three, the weighted gap factor will be noticeably affected. The correlations r c will still be close to one for all communities, but the division by the number of communities (which in this situation is n > 1) will shift the gap factor values to the right (this is the second vertical band). For example if n = 2 and r c ≈ 1 for all communities, then g s will be close to 0.5 when one community is large (size s max ) and the other community satisfies s c s max . A similar issue causes the diagonal striations seen near the top of FIG. 3a. Variations in packing structure will again cause the number and size of communities to differ slightly between different experimental configurations at the same pressure. When the number and size of communities is few (i.e. at high resolutions), the gap factor is sensitive to such small differences and its values can take on a large spread (this causes the set of diagonal bands across the gap factor axis). The hull ratio is more robust against these small differences, and so for fixed resolution parameter and pressure, its values are contained to a smaller range. Both of these features disappear at intermediate values of γ when there are more communities and many communities are intermediately sized. An alternative definition of community averaging in the weighted gap factor could alleviate the large spread in values that can occur for similar packings at low and high resolutions.

A. Force chain identification
We now apply each of these statistics to collections of laboratory and simulated packings under a variety of pressure conditions ranging from 2.7×10 4 E to 5.9×10 3 E. Our first goal is to examine how well these measures can identify force chain architecture from the results of the community detection methods. We compute the value of each statistic over the eleven different resolution parameters between γ = 0.1 and γ = 2.1, and report averages over all optimizations and packings. Beginning with the gap factor, we follow the method described in [1] and extract force chains at each pressure by identifying an optimal value of γ that maximizes the gap factor (FIG. 4a). For the experimental data set, we find that the optimal region occurs between γ = 0.7 and γ = 0.9, and these resolutions also give rise to communities that strongly resemble the force chain architecture seen in the photo elastic disk experiments (FIG. 5a). One can identify optimal resolutions for the simulations as well (FIG.  4b), and we find that these values occur predominantly between γ = 1.1 and γ = 1.3.
A similar process can be carried out with the hull ratio, but this time we take the optimal value of γ at each pressure to be that which minimizes R s . The resolutions that minimize the hull ratio should give rise to more heterogenous structure and less compact and linear structure. The optimal values lie between γ = 0.7 and γ = 0.9 for the experimental data (FIG. 4c), and between γ = 0.9 and γ = 1.1 for the simulations (FIG. 4d). The extrema of the hull ratio are very well defined for the laboratory packings (even more so than those of the gap factor), and it is thus clear where the optimal resolution parameters lie. In addition, the force chain region identified from the hull ratio agrees well with that identified from the gap factor. These findings suggest that physical distance information alone -as instantiated in the hull ratio -may be sufficient for extracting force chain structure once we have constructed communities using the inter-particle contact forces.
The TCF depends only on the network topology and does not take into consideration any physical distances. It is somewhat less straightforward to use this quantity to identify the force chain region of the structural resolution parameter γ than it is with either the gap factor or hull ratio. In the experimental data, the resolution parameters that correspond to maxima in the TCF occur at values similar to those found to maximize the gap factor and minimize the hull ratio (FIG. 4e). In contrast, in the simulated data the resolution parameters that correspond to maxima in the TCF occur at very small resolutions, but local maxima correspond to values found to maximize the gap factor and minimize the hull ratio (FIG. 4f).
Our analysis shows that from the three statistics we examined, those that include notions of physical shape are most adept at quantifying the presence of branching and other heterogeneities known to occur in force chain structure. Both the gap factor (which uses a combination of physical and topological distance information) and the hull ratio (a purely geometric measure) have well defined extrema at resolution parameters that result in groups of particles reminiscent of the force chains we see by eye. On the other hand, the purely topological TCF is less able to identify these optimal resolutions in the experimental data. Based on the outcome of our examination, we consider the force chain region for the laboratory data to occur between resolution parameters of γ = 0.7 and γ = 0.9. For the simulated data, we take a larger range for the force chain region, between γ = 0.9 and γ = 1.3.

B. Sensitivity to packing pressure
In the previous section, we explored the ability of three statistics to identify force chain-like structure in experimental and simulated granular packings over a range of pressures. Another important question is if those same measures can detect differences in packing architecture as a function of the confining pressure, or if they are insensitive to pressure changes. As observed in [1], under increasing pressure, the communities in the force network have (on average) smaller geometric radius for a given resolution parameter γ. In the resolution range where the resulting communities qualitatively resemble force chains (FIG. 5a), it is natural to expect features of the packings that rely on force chain structure to be visible.
We observe that R s and T s both stratify real and simulated packings across pressures, though they do so in different ranges of the resolution parameter. In contrast, the gap factor does not distinguish pressures particularly well at any resolution (FIG. 4a, b), as observed in [1]. For the simulated packings especially, the curves corresponding to each pressure collapse onto one another, making it difficult to identify the relative pressure on the system from the measured statistics. For the laboratory packings, the hull ratio distinguishes pressure best when γ > 1.3 (FIG. 4c). In particular, we observe that in the stratifying regions, the hull ratio increases as pressure increases. This behavior agrees with the notion that increasing pressure generally leads to more compact structure. However, it is important to point out that the stratifying range lies outside the resolution range where communities are most representative of the force chains we observe visually. In FIG. 5b, we show an example of the community structure in the region where the hull ratio can distinguish pressure. The modules are smaller in size and number compared to those extracted at the optimal resolution (FIG. 5a). On the other hand, the purely topological statistic best distinguishes pressure near γ ≈ 1 (FIG. 4e), which is much closer to the force chain region identified in the previous section. In this γ range, the community structure in the experimental packings resembles the types of force chains one might pick out by eye (FIG. 5c). For the simulated packings, the hull ratio can distinguish between most pressures near the minima (i.e. the force chain region), but the relative separation of the curves is quite small (FIG. 4d). Of the three statistics, the TCF is the most sensitive to pressure differences (especially to higher pressures) near the force chain region in the simulations (FIG. 4f).
To understand these results, observe that in both the laboratory and simulated packings, the number of distinct force chains detected drops, while their size increases with pressure ([1], FIG. 6). In the resolution range where communities are most representative of force chains, at higher pressures the resulting chains are long, but branch irregularly. The gap factor, defined as the correlation between the hop distance and the geometric The systemic gap factor as a function of resolution parameter γ across varying pressure on laboratory (a) and simulated (b) packings. Curves are means across trials, and error bars indicate one standard deviation. The maxima of the gap factor correspond to the optimal resolution parameters for identifying community structure reminscent of force chains; but in both the laboratory and simulated packings, the gap factor is relatively homogeneous across pressure. (c, d) The weighted hull ratio as a function of resolution parameter γ and pressure for the laboratory (c) and simulated (d) packings.
The optimal γ is that which minimizes R s , and agrees relatively well with the optimal γ found from the gap factor in both experiments and simulations. However, the laboratory weighted hull ratio can distinguish pressure across packings for large values of γ > 1.3. The simulated hull ratio also partially distinguishes pressures for high values of γ, but is less well differentiated.
(e, f ) The systemic topological compactness factor T s as a function of the resolution parameter γ for the laboratory (e) and simulated (f ) packings. T s provides a strong discrimination of packing pressures at the optimal resolution obtained from the hull ratio. The high values of the TCF for γ ≈ 0.1 is the result of the small number of compact communities which rapidly fragment as γ increases.
distance between nodes in each force chain, will have difficulty detecting the difference between long chains with several short branches and those that are essentially linear. However, the topological measurement treats all branches equally and, as a result, can detect the increase in such motifs. Once γ > 1.3, under high pressure the tendency of strong communities will be to form tightly packed clumps, as opposed to curved structures with large convex hulls. In the experimental setting, the friction between particles can effectively protect branching chains that would otherwise collapse in frictionless numerical settings.

VI. DISCUSSION
A necessary first step toward understanding the physical properties of packed granular materials is the development of statistics that provide insight into the structure of models. Here, we have extended the data driven, network theoretic approach to the study of force chains initiated in [1], extracting putative chains using techniques from community detection. In order to understand their structure, we consider three statistics: the previously defined gap factor, along with a purely geometric measure called the hull ratio which captures much of the same information without relying on network topology, and a novel algebraic-topological measurement called the topo- logical compactness factor which distinguishes pressure in both experimental and simulated packings, showing that the purely topological properties of the force chains vary in a predictable manner with changes in pressure. a. Network-Based Tools for Assessment of Material Architecture Network science provides a natural framework in which to both represent and characterize granular materials. In this paradigm, particles are nodes in a graph and inter-particle forces are edges. Once such a graph is constructed, a number of network-based measures can be used to assess the material architecture at different spatial scales. Network science approaches seem to be especially advantageous in their ability to probe mesoscale features of granular media, of which force chain structure is a prime example. This heterogenous architecture constrains the mechanical properties of granular media under compression and shear [6], and it may also constrain the heterogeneous and nonlinear nature of acoustic signal transmission that particulate and continuum models often fail to describe [9,[46][47][48][49][50][51][52][53][54].
A primary aim of the present work was to thus suggest three distinct measures for the characterization of crucial mesoscale organization in granular networks, and explore their sensitivity to pressure differences in experimental and laboratory systems. A number of previous studies have investigated intermediate scale features of granular systems as well. A large focus has been on the role of contact loops of three or more particles under various conditions. For example, using a graph representation of a granular packing, Rivier [44] demonstrated that the stability of the packing was related to the dynamical frustration of odd cycles in the network. In a study on the co-evolution of cycles and force chains [17], 3-cycles were found to be stabilizing mesoscale structures in the lead-up to force chain buckling [55]. Cycles may also characterize the jamming transition; in particular the number of triangles grows suddenly at the critical packing fraction in simulations of isotropically com-pressed granular systems, suggesting that the presence of 3-cycles is characteristic of a stable, rigid state [15]. The evolving structure of contact loops has been investigated in simulations of tilted granular packings as well [13]. In a study on simulations of tapped granular packings, Arévalo et al. [18] found that polygons in the network could distinguish between equilibrium states of the same packing fraction. Spatially embedded communities probe another mesoscale feature of granular packings. Community structure can predict certain features of acoustic transmission in experimental systems [27], and when determined with a physically informed null model, network communities result in force chain-like structure which can be subsequently analyzed and compared across varying pressures and different types of systems [1]. Examples of other recent studies on granular media that have utilized a network-based analysis to probe local, intermediate and global scales include [14,25,26,28,29,[56][57][58].
The findings mentioned above demonstrate that network science is a useful tool in understanding the complex material properties of granular media. In particular, they emphasize the general need to define new measures that can quantify heterogeneous architectures between the particle-level and bulk scales. An exciting direction for future work would be to examine how the physically informed method of community detection introduced in [1] and used in this study, as well as the statistics introduced here and in prior work, can explain or predict interesting features of acoustic propagation in granular experiments.
b. New Insights from Algebraic Topology Algebraictopological methods have previously been applied to the study of force networks in granular media [21]. The previous approach involved application of persistent homology to the complete, weighted force network of the system (as well as related networks, such as weighting by tangent forces) to extract quantitative measures of the global network structure. Studying these global statis-tics, the authors see significant differences in the evolution of components and loops in the network as packing density changes, and observe a severe decrease in the rate of change of these measures as the system passes the jamming transition. Here, instead, we describe a highly local measure of rigidity in structure of individual force chains, and then aggregate these to obtain an estimate of rigidity and stability throughout the entire packing.
This algebraic-topological approach provides a powerful tool set that that can be extended in a variety of ways. Using the definition of the statistic in terms of the general framework of cliques and chain complexes in Appendix A, it is straightforward to extend this approach to 3D packings by including cliques on four vertices and computing H X, 2 (B) and its corresponding normalization. Even more strikingly, however, because we have lifted the computation to the level of vector spaces with bases defined in terms of elements of the system, we can study how these vector spaces evolve as the system does by translating these changes into linear maps. Techniques from topology such as zig-zag persistence [59], which describes how to "glue together" homology from related networks, can then be used to study the evolution of homological features, thus providing deeper insights into the properties of the dynamic system than simple numerics.
c. Implications for Network Design In future studies, we expect network-based approaches to directly impact the field of material design, which aims to engineer novel materials with desired and optimized physical properties. One notable advantage of network science is that the mathematical tools are agnostic to the physical makeup of the system under inspection, thus offering a framework in which to study a diversity of materials ranging from granular media to biological tissues [60]. This versatility may be especially promising for an emerging branch of material design concerning so-called metamaterials [61][62][63], which are developed via precise control of shape, geometry, orientation, and arrangement, rather than by choice of specific material composition. Metamaterials come in many forms; for example, some are designed with tailored mechanical and acoustic properties [63][64][65][66], others with desired electromagnetic properties [61,67,68], and more. Graph theoretic approaches have potential to aid in the development of numerous materials that have underlying network topology.
From a material design perspective, a second advantage of network science is that it provides a means to quantify structure on all scales, including mesoscale architecture in heterogeneous and disordered materials. As noted previously, this type of structure is often crucial for determining material properties. Furthermore, it is not limited to the granular regime; for example, heterogeneities and structure on the intermediate scale are also known to be important in biological materials (for some examples, see [60,[69][70][71][72][73][74]). One can imagine using a network-based framework to purposefully design heterogeneous materials with desired mesoscale topology that gives rise to unique physical properties. To our knowl-edge, this would be a new perspective, but recent and on-going work in the development of physically informed network tools and models [1,27,32,75] are bringing these objectives closer to fruition.
One approach for material design in a network-based framework may be to utilize multi-objective functions and Pareto optimality [76]. For instance, starting with a physically realizable material network structure and a known correlation between a network statistic and a material property of interest (mechanical stability, for example), one can rewire the network with a cost function to maximize the network statistic that is correlated with the property of interest and anti-correlated with any properties of non-interest (instabilities, for example). One way to explicitly do the rewiring is along Pareto optimal fronts in a network morphospace [77], with rewiring rules that preserve network features and topologies required by physical constraints and laws [78][79][80], and references therein. This method would theoretically select for physically feasible material network designs while simultaneously optimizing topologies to result in purposefully selected material properties. the string v 0 v 1 . . . v k . Cliques with (k +1) vertices can be realized as the convex hull of k points, which generically span a k-dimensional region, which leads to a choice of index; it is standard to index by dimension rather than number of vertices, and we will follow that convention here. Now, the set of cliques decomposes as a disjoint union of cliques of fixed size, and we define sets For our purposes, the most important of these will be X 0 (G), the set of vertices of G; X 1 (G), its edges; and X 2 (G), its triangles, as no larger cliques can appear in a planar network. However, the underlying principle can be generalized to any network, and these may possess larger cliques. Now, construct a sequence of vector spaces and linear maps [82] called the clique chain complex of G, Define the vector spaces C X k (G) to be a vector space with basis elements {e v0v1...vi s.t. v 0 v 1 . . . v i ∈ X i (G)} corresponding to cliques in G. Observe that removing any vertex (and all attached edges) from a clique results in a smaller clique, and that geometrically the collection of all sub-cliques on one fewer vertex forms the boundary of the original clique. This geometry provides us with the definition of the boundary maps δ k : C X k (G) → C X k−1 (G), defined on basis vectors as an alternating sum [83] of all sub-cliques with one vertex fewer, where the hat denotes omission of a vertex, and extended linearly to all of C X i (G). By convention, define δ 0 to be the zero map.
There are two distinguished classes of elements in the vector spaces in the chain complex: cycles, which are elements in the kernel (or null space) of some δ map, and boundaries, which are elements in the image of some δ map. It turns out that all boundaries are cycles: this is perhaps clear geometrically, but one can simply compute that, δ k−1 • δ k is the zero map for each k: Thus, the image of δ k+1 is a vector subspace of the kernel of δ k , and we define the kth homology group of the clique chain complex to be the quotient vector space Elements of the kth homology group are equivalence classes of k-cycles, and the members of each equivalence class differ by boundaries of (k + 1)-cliques. In the twodimensional setting, there is a canonically defined minimal representative of each non-zero equivalence class, which we take to correspond to a fundamental circuit in the force chain (FIG. 2c (right, yellow)). Compact regions will induce cycles that consist of cliques, and these will be equivalent to the trivial cycle (FIG. 2c (right,  blue)). While this machinery in general provides a useful framework in which to discuss force chain topology, we note that the clique homology groups as defined do not fully capture the interesting branching structure of the graph. In particular, branches that do not close to form circuits are ignored (FIG. 2c (right, red)). To repair this defect, we consider instead the quotient of the clique complex obtained by identifying all of the leaves { 1 , 2 , . . . k } ⊆ [N ] (FIG. 2c (middle, green)) to a single point, which creates new cycles that correspond to such branches (FIG. 2c (right, purple). However, the resulting object may not be a clique complex (if this would result in paths between leaves of length two or less, there would be multiple edges between adjacent vertices); to avoid this complexity, we can instead take the quotient on the level of the chain complex, defining a new chain complex where C X, 0 (G) = C X 0 (G)/ e 1 −e 2 , e 1 −e 2 , . . . e 1 −e k is the quotient vector space in which all basis elements corresponding to leaves are identified to a single basis vector and δ 1 is the map induced on the quotient from δ 1 . Now, define the leaf-reduced clique homology groups, H X, k (G) to be the homology groups of this chain complex.
Finally, we consider a numerical measurement of these loops, the leaf-reduced clique Betti numbers of G, given by β X, k (G) def = rank(H X, k (G)).
Because the potential size of a homology group scales non-linearly with the size of the underlying graph, it is expedient to look at an appropriate normalization. The presence of cliques corresponds to the presence of compact regions, and thus makes force chains more stable, while the branch points corresponding to non-trivial cycles reduce stability, so we define the topological compactness factor to be T CF (G) def = 1 − β X, 1 (G) rank(ker(δ 1 )) which corresponds to the graph-theoretic description given in Section IV C.
There are several advantages to this linear-algebraic approach. First, it is straightforward to define an anal-ogous statistic for the second (and higher) Betti numbers for use in the 3D case or for use in non-spatial networks. More interesting, however, is the fact that graph homomorphisms induce linear maps on homology. These maps, in turn, can provide information about the evolution of the topological structure of the system under, for example, shear stress or chagnge in packing density.