# Network analysis of particles and grains

Network analysis of particles and grains Abstract The arrangements of particles and forces in granular materials have a complex organization on multiple spatial scales that range from local structures to mesoscale and system-wide ones. This multiscale organization can affect how a material responds or reconfigures when exposed to external perturbations or loading. The theoretical study of particle-level, force-chain, domain and bulk properties requires the development and application of appropriate physical, mathematical, statistical and computational frameworks. Traditionally, granular materials have been investigated using particulate or continuum models, each of which tends to be implicitly agnostic to multiscale organization. Recently, tools from network science have emerged as powerful approaches for probing and characterizing heterogeneous architectures across different scales in complex systems, and a diverse set of methods have yielded fascinating insights into granular materials. In this article, we review work on network-based approaches to studying granular matter and explore the potential of such frameworks to provide a useful description of these systems and to enhance understanding of their underlying physics. We also outline a few open questions and highlight particularly promising future directions in the analysis and design of granular matter and other kinds of material networks. Glossary of Terms Granular materials: Granular materials are collections of discrete, macroscopic particles that interact with each other through contact (rather than long-range) forces. Importantly, these systems are non-equilibrium: the particles are large enough to avoid rearrangement under thermal fluctuations, and they lose energy through frictional and inelastic interactions with neighboring particles. Particulate materials: Like granular materials, particulate materials are collections of discrete, macroscopic elements. However, the elements making up the system may be entities—such as bubbles, foams, colloids, or suspensions—that include multiple phases of matter. The term “particulate material” is a more general one than “granular material”. Packing fraction: The fraction of a granular material that is physically occupied by the particles. One calculates the packing fraction as the ratio of the total volume of all particles to the total volume of the region that contains the particles. The packing fraction is also sometimes called the “packing density” or “volume fraction”. Force chain: Force chains are typically described as the subset of inter-particle contacts in a granular material that carry the largest forces in the system. They often form filamentary networks that align preferentially with the principal stress axes under which a material is loaded. Jamming: As certain system parameters change, disordered, particulate materials typically undergo a transition from an underconstrained, liquid-like state to a rigid, solid-like state characterized by the onset of mechanical stability. The transition to/from a jammed state often arises through increases/decreases in packing fraction (or in contact number), which can occur due to an applied load. A formal theory of jamming exists for idealized situations (with soft, frictionless, and spherical particles). Isostatic: A jammed packing is isostatic when it has exactly the minimum number of contacts that are required for maintaining mechanical stability through force balance and torque balance. One typically examines isostaticity by calculating the mean number of contacts per particle. “Hyperstatic” and “hypostatic” packings have more and fewer contacts, respectively. Mono/bi/poly-disperse: A particulate material is monodisperse if it is composed of particles of a single species (i.e., particles with the same size, shape, and material properties). Bidisperse materials have particles of two species, and polydisperse materials can either have particles from three or more discrete species or have particles with a continuum of properties. Structural rigidity theory: For a specified mechanical structure that is composed of fixed-length rods connected to one another by flexible hinges, structural rigidity theory studies the conditions under which the associated structural graph (in which the edges correspond to the rods and the nodes correspond to the hinges) is able to resist deformations and support applied loads. Stress: A stress is a force applied to an object’s surfaces. (The units measure a force per unit area.) Shear stress arises from the component of the applied force that acts in a direction parallel to the object’s cross-section, and normal stress arises from the perpendicular component. Strain: Strain is the fractional (unitless) deformation of a material that arises due to an applied stress. One calculates strain from the relative displacement of particles in a material, excluding rigid-body motion such as translation or rotation. Like stress, strain has both shear and normal components. Pure shear: One can use the term “pure shear” to describe either stresses or strains in which an object is elongated along one axis and shortened in the perpendicular direction without inducing a net rotation. Axial compression: In axial compression, one applies inward forces to an object in one direction (uniaxial), two directions (biaxial), or all directions (isotropic compression). These forces result in uniaxial strain, biaxial strain, or isotropic strain, respectively. Cyclic shear/compression: These terms refer to procedures in which repeated cycles of shear or compression are applied to the same system. Shear band: A shear band is a narrow region of a particulate material in which most of the strain is localized, whereas other regions remain largely undeformed. A shear band is also sometimes called a region of “strain localization”. Strain softening/hardening: As a material is loaded and undergoes deformation, continuing deformation can become either easier (strain softening) or harder (strain hardening). Eventually, after much deformation, the material can reach a critical state in which there are no further changes in the resistance to deformations. Stress ratio: The stress ratio, which is analogous to Coulomb’s Law, is the ratio of shear to normal stresses. Frictional failure occurs when the shear force exceeds the product of the normal force and the coefficient of friction. Photoelasticity/birefringence: Photoelasticity is an optical technique for quantifying internal stresses in a material. This experimental methodology is based on the transmission of polarized light through “birefringent” materials, which have preferentially fast and slow directions for the propagation of light. DEM or MD simulations: The Discrete (or Distinct) Element Method and Molecular Dynamics simulations are related numerical techniques that compute the motions of all particles in a system (such as in a granular material). In both methods, a computer algorithm treats each particle as an object subject to Newton’s laws of motion, where forces consist of body forces (e.g., gravity) and those that arise from interactions with the object’s neighbors. 1. Introduction Granular materials comprise a subset of the larger set of particulate matter [1–6]. People engage with such materials—which include sands, beans, grains, powders such as cornstarch, and more—often in their daily lives. One can define a granular material as a large collection of discrete, macroscopic particles that interact only when in contact. Granular materials are inherently non-equilibrium in two distinct ways, characterized by (1) the lack of rearrangement under thermal fluctuations and (2) the loss of energy through frictional and inelastic dissipation during contact between grains. Nonetheless, they phenomenologically reproduce equilibrium states of matter, exhibiting characteristics of solids (forming rigid materials), liquids (flowing out of a container) or gases (infrequent contacts between grains), depending on the type and amount of driving. In this review, we focus mainly on granular solids and slow (non-inertial) flows [7]; these are dense materials in which sustained inter-particle contacts provide the dominant contribution to material properties. The functional properties of granular materials are related in a non-trivial way to the complex manner in which particles interact with one another and to the spatial scales (particle, chain, domain and bulk) and time scales over which those interactions occur. For example, pairs of particles can exert force on one another in a local neighbourhood. However, as particles push on adjacent particles, the combined effect can transmit forces over long distances via important mesoscale structures commonly called force chains [8, 9]. The idea of networks has been invoked for many years to help provide a quantitative understanding and explanation of force-chain organization [10–14]. Broadly speaking, force chains form a network of filamentary-like structures that are visually apparent in images from experiments, like the one shown in Fig. 1. In such images, the brighter particles carry larger forces [8, 16]. Furthermore, force chains tend to align preferentially along the principal stress axes [17]. It can be helpful to think of a force-chain network as the backbone of strong forces that span a system, providing support for both static [18] and dynamic [16] loading. However, weaker forces can also play a stabilizing role, much as guy-wires do on an aerial tower [19, 20]. Fig. 1. View largeDownload slide Force chains in an experimental granular system. A photoelastic image of a quasi-two-dimensional (quasi-2D) packing of photoelastic disks that were subjected to pure shear. The photoelastic image allows one to visualize the force pattern in a material. Bright particles carry the strongest forces, and one can observe that a network of force chains tends to align along the principal stress axes. Modern techniques allow one to determine vector contact forces at each inter-particle contact. We adapted this figure, with permission, from [15]. Fig. 1. View largeDownload slide Force chains in an experimental granular system. A photoelastic image of a quasi-two-dimensional (quasi-2D) packing of photoelastic disks that were subjected to pure shear. The photoelastic image allows one to visualize the force pattern in a material. Bright particles carry the strongest forces, and one can observe that a network of force chains tends to align along the principal stress axes. Modern techniques allow one to determine vector contact forces at each inter-particle contact. We adapted this figure, with permission, from [15]. It is also possible for sets of particles to cluster together into larger geographical domains, with potentially distinct properties, that can have weak structural boundaries between them [21]. At the largest scale, granular materials as a whole exhibit bulk properties, such as mechanical stability or instability in response to sheer or compression [22]. All the aforementioned spatial scales are potentially relevant for understanding phenomena such as transmission of acoustic waves [23], thermal conductivity and heat transfer [24], electrical properties [25] and more. The time scales of interactions in granular materials are also important, and they can vary over many orders of magnitude. For example, in systems under compression, statistical fluctuations of grain displacements depend fundamentally on the length of the strain step (i.e., ‘increment’) over which one makes measurements, as fluctuations over short windows are consistent with anomalous diffusion and those over longer windows are consistent with Brownian behaviour [26]. The principled study of such diverse characteristics and organization in a single system can be very challenging, and the development of appropriate physical, mathematical, statistical and computational models is important to attain a mechanistic understanding of granular materials. Traditionally, it has been common to model granular materials using either particulate-based or continuum-based frameworks [20]. However, both of these approaches are often implicitly agnostic to intermediate-scale organization, which is important for understanding both static granular packings [27] as well as granular dynamics [28]. Recently, tools from network science [29, 30] and related mathematical subjects—which include approaches that can account explicitly for mesoscale structures [31–34]—have been used successfully to study properties of granular materials across multiple spatial and temporal scales. The most common representation of a network, an important idea for the study of complex systems of interacting entities [35], is as a graph [30]. A graph consists of a set of nodes (to represent the entities) and a set of edges, each of which represents an interaction between a pair of entities (or between an entity and itself). Increasingly, more complicated network representations (such as multilayer networks [36]) are also employed. Moreover, there is also a growing recognition that it is important to consider the impact of other features, such as spatial embedding and other spatial effects [37], on network structure and dynamics, rather than taking an approach that promises that ‘one size fits all’. Network science offers methods for quantitatively probing and analysing large, interacting systems whose associated networks have heterogeneous patterns that defy explanations attained by considering exclusively all-to-all, regular, or lattice-like interactions [29]. There are several open problems in granular physics that may benefit from network-science approaches. In particular, because granular materials have multiple relevant length and time scales [16, 38–41], it can be challenging to model and quantify their structural organization, material properties and responses to external loads [42–46]. However, although complex, the pairwise inter-particle interactions that underlie and govern the structure and behaviour of granular systems (and other particulate matter) render them amenable to various network representations and network-based analyses. Smart and Ottino [47] were among the first to explicitly suggest and formalize the use of ideas from network science to begin to study some of the difficult questions in granular physics. In their article, they highlighted the ability of a network-based perspective to complement traditional methods for studying granular materials and to open new doors for the analysis of these complex systems. One place in which network analysis may be especially useful is in quantifying how local, pairwise interactions between particles in a granular packing yield organization on larger spatial scales (both mesoscale and system-level). For example, in sheared or compressed granular packings, such organization can manifest as force chains or other intermediate-sized sets of particles that together comprise a collective structure. Network science provides approaches to extract and quantitatively characterize heterogeneous architectures at microscale, mesoscale and macroscale sizes, and one can use these methods to understand important physical phenomena, including how such multiscale organization relates to bulk material properties or to spatial and temporal patterns of force transmission through a material. Network-based approaches should also be able to provide new insights into the mechanisms that govern the dynamics of granular materials. For example, as we will discuss, network analysis can provide quantitative descriptions of how the structure of a dense granular material evolves as a system deforms under external loads (such as those induced by compression, shear, tapping or impact), and it can be helpful for describing certain aspects of complex dynamics (such as granular flows). Furthermore, it is important to consider the relative time scales of dynamics (e.g., sound propagation) on networks and the evolution of network structure. We also expect ideas from temporal networks [48] or adaptive networks [49] to be fruitful for studying faster structural dynamics, and investigation of granular dynamics in general should benefit from the development of both novel network representations and methods of network analysis that are designed specifically to understand temporally evolving systems. Another important problem in the study of granular materials is to predict when and where a granular system will fail. There has been some progress made in this area using network-based approaches, but it is important to continue to develop and apply tools from network analysis and related areas to gain a deeper understanding of which network features regulate or are most indicative of eventual failure. Another exciting direction for future work is to combine network-based approaches with questions about material design. In particular, can one use network-based approaches to help engineer granular systems—or other materials that are amenable to network representations—with desired and specialized properties? It is also important to note that network-based representations and methods of analysis can provide insightful descriptions of granular materials with various additional complexities, such as systems that are composed of differently-shaped particles, three-dimensional (3D) materials and so on. This flexibility makes the application of tools from network science a powerful approach for studying the structural properties and dynamics of granular systems. Such a framework also allows one to compare network architectures in diverse situations, such as between simulations and experiments, across systems that are composed of different types of particles or exposed to different loading conditions, and more. Exploiting these capabilities will yield improved understanding of which properties and behaviours of granular materials are general, versus which are specific to various details of a system. With the continued development of physically-informed network-analysis tools, network-based approaches show considerable promise for further development of both qualitative and quantitative descriptions of the organization and complex behaviour of granular materials. The purpose of our article is to review the nascent application of network theory (and related topics) to the study of granular materials. We begin in Section 2.1 with a mathematical description of networks. In Section 2.2, we briefly review a set of measures that one can calculate on graphs and which have been useful in past investigations of granular materials. In Section 3, we review several ways in which granular materials have been represented as networks, and we discuss investigations of such networks to quantify heterogeneous, multiscale organization in granular materials and to understand how these systems evolve when exposed to external perturbations. We also point out insights into the underlying physics that have resulted from network-based investigations of granular matter. We close in Section 4 with some thoughts on the many remaining open questions, and we describe a few specific future directions that we feel are important to pursue. We hope that our review will be helpful for those interested in using tools from network science to better understand the physics of granular systems, and that it will spur interest in using these techniques to inform material design. 2. Network construction and characterization 2.1 What is a network? It is often useful to model a complex system as a network, the simplest type of which is a graph [29, 50]. For most of our article, we will use the terms network and graph synonymously, but the former concept is more general than the latter.1 A graph $$G$$ consists of nodes (i.e., vertices), where pairs of nodes are adjacent to each other via edges (i.e., links). We denote the set of nodes by $$\mathscr{V}$$ and the set of edges by $$\mathscr{E}$$. A node can also be adjacent to itself via a self-edge (which is also sometimes called a self-loop), and a multi-edge describes the presence of two or more edges that are attached to (i.e., incident to) the same pair of nodes. (Unless we state otherwise, we henceforth assume that our networks have neither self-edges nor multi-edges.) The number of nodes in a graph is the size of the graph, and we also use the word ‘size’ in the same way for other sets of nodes. A subgraph of a graph $$G$$ is a graph constructed using a subset of $$G$$’s nodes and edges. An edge between two nodes represents some sort of relationship between them. For example, edges can represent information flow between different parts of the internet [53], friendship or other social interactions between people [54], trading between banks [55], anatomical or functional connections between large-scale brain regions [56, 57], physical connections between particles in contact [21, 58] and so on. Edges can be either unweighted or weighted, and they can be either undirected or directed [29]. In an unweighted (i.e., binary) network, an edge between two nodes is assigned a binary value (traditionally $$0$$ or $$1$$) to encode the absence or presence of a connection. In a weighted network, edges can take a variety of different values to convey varying strengths of relationships between nodes. In an undirected network, all edges are bidirectional, so one assumes that all relationships are reciprocal. In a directed network, however, edges have a direction that encodes a connection from one node to another. An adjacency matrix is a useful way to represent the information in a graph. For an unweighted and undirected graph, an element of the adjacency matrix $$\mathbf{A}$$ of an $$N$$-node network is $$A_{ij} = \left\{\begin{array}{@{}ll} 1, \text{ if there is an edge between nodes}\,i\, \text{and}\,j, \\[3pt] 0, \text{ otherwise, } \end{array} \right.$$ (2.1) where $$i,j \in \{1, \dots, N\}$$. For a network in which nodes do not have labels, one can apply a permutation to $${\bf A}$$’s rows and columns—one uses the same permutation for each—to obtain another adjacency matrix that represents the same network. In this article, when we refer to ‘the’ adjacency matrix $${\bf A}$$ of a graph with unlabelled nodes, we mean any one of these matrices, which have the same properties (spectra, etc.). For a weighted graph, if nodes $$i$$ and $$j$$ are adjacent via an edge, we denote the corresponding edge weight by $$w_{ij}$$ (which is usually given by a non-negative real number (e.g., the value of the normal or tangential component of the force between two contacting particles)2). An element in the associated weighted adjacency matrix (which is sometimes called a weight matrix) $$\mathbf{W}$$ is $$W_{ij} = \left\{\begin{array}{@{}ll} w_{ij}, \text{ if there is an edge between nodes}\,i\, \text{and}\, j, \\[3pt] 0, \text{ otherwise.} \end{array} \right.\$$ (2.2) For the more general case of a weighted, directed graph, if there is an edge from node $$j$$ to node $$i$$, then we let $$w_{ij}$$ represent the weight of that edge [29]. The associated weighted and directed adjacency matrix $$\mathbf{W}$$ is $$W_{ij} = \left\{\begin{array}{@{}ll} w_{ij}, \text{ if there is an edge from node}\,j\, \text{to node}\, i, \\[3pt] 0, \text{ otherwise.} \end{array} \right.\$$ (2.3) An adjacency matrix associated with an undirected network is symmetric, but an adjacency matrix for a directed network need not be (and is symmetric if and only if all directed edges are reciprocated). In the present review, we primarily consider undirected networks, although we will occasionally make remarks about directed situations. For weighted graphs, it is often also important to consider a binary adjacency matrix $$\mathbf{A}$$ associated with a weight matrix $$\mathbf{W}$$. Note that $$\mathbf{A}$$ captures only the connectivity of nodes (i.e., their adjacencies), irrespective of how strongly they interact with each other. In terms of $$\mathbf{W}$$, the corresponding binary network (which can be either directed or undirected) is $$A_{ij} = \left\{\begin{array}{@{}ll} 1, \text{ if } W_{ij} \neq 0, \\[3pt] 0, \text{ otherwise.} \end{array} \right.\,$$ (2.4) It is common to use terms like network topology when discussing structural properties of $$\mathbf{A}$$, and sometimes one uses terms like network geometry when discussing properties that also depend on edge weights. Because we will also employ ideas from subjects like algebraic topology (see Section 2.2.10), we will need to be very careful with such terminology. The network representations that have been used to study granular matter (and other kinds of materials) employ diverse definitions of edges (both weighted and unweighted and both directed and undirected), and some generalizations of graphs have also been considered. See Fig. 2 for a schematic showing possible choices of nodes and edges for a network representation of a granular packing. A variety of tools and measures from network analysis have been used to study granular networks. We discuss some of these ideas in Section 2.2. Fig. 2. View largeDownload slide From a packing to a network. (a) A sample packing of grains. (b) Representation of particles as network nodes. (c) Representation of contacts as network edges. (d) Graph representation of nodes and edges. The edges highlighted in green illustrate a cycle, which is a loop of physical contacts. The edges highlighted in peach illustrate a set of triangular loops, which are minimally-rigid structures in the context of structural rigidity. Edge weights often represent contact forces, which we illustrate here with different line widths. The degree $$k$$ of a node is equal to the number of edges attached (i.e., incident) to that node, and the strength (i.e., weighted degree) $$s$$ of a node is given by the sum of the weights of the edges attached to that node. One can (e) encode an unweighted (i.e., binary) graph as an unweighted adjacency matrix and (f) encode a weighted graph as a weighted adjacency matrix. Fig. 2. View largeDownload slide From a packing to a network. (a) A sample packing of grains. (b) Representation of particles as network nodes. (c) Representation of contacts as network edges. (d) Graph representation of nodes and edges. The edges highlighted in green illustrate a cycle, which is a loop of physical contacts. The edges highlighted in peach illustrate a set of triangular loops, which are minimally-rigid structures in the context of structural rigidity. Edge weights often represent contact forces, which we illustrate here with different line widths. The degree $$k$$ of a node is equal to the number of edges attached (i.e., incident) to that node, and the strength (i.e., weighted degree) $$s$$ of a node is given by the sum of the weights of the edges attached to that node. One can (e) encode an unweighted (i.e., binary) graph as an unweighted adjacency matrix and (f) encode a weighted graph as a weighted adjacency matrix. 2.2 Some tools for characterizing granular networks Network theory [29] provides myriad ways to characterize and quantify the topological and geometrical organization of complex networks. Thus, different network methods can reveal different important features of the underlying system, and these features, in turn, can help explain how a system behaves in certain situations. In the context of granular matter, for example, it is often desirable to understand the stability of a material, mechanical responses to external stresses, or wave propagation through a system. Recent investigations have demonstrated that network analysis can inform understanding of the mechanisms that underlie these phenomena. In this section, we discuss several network concepts; and in Section 3, we describe how they have been used for the study of granular materials. We are, of course, not presenting an exhaustive list of methods from network science. See [29] and other books and reviews (and references therein) for discussions of other tools. For simplicity, we primarily give definitions for undirected networks, though many of the ideas that we present also have directed counterparts. We start with basic network diagnostics, and also discuss some more complicated methods. 2.2.1 Degree. One local property of a network is node degree. In an undirected network, a node’s degree is equal to the number of edges that are attached to it (see Fig. 2d). We denote the degree of node $$i$$ as $$k_{i}$$, and we recall that $$N$$ denotes the total number of nodes. For an unweighted graph with adjacency matrix $${\bf A}$$, one can calculate $$k_{i}$$ with the formula $$k_{i} = \sum_{j=1}^{N} A_{ij}.$$ (2.5) One can generalize the idea of degree to strength (i.e., weighted degree) using a weight matrix $$\mathbf{W}$$ [59, 60]. The strength $$s_i$$ of node $$i$$ is equal to the sum of the weights of the edges that are attached to that node: $$s_{i} = \sum_{j=1}^{N} W_{ij}.$$ (2.6) Its associated degree $$k_i$$ is still given by Eq. (2.5). A common network representation of granular materials is to treat particles as nodes and physical contacts between particles as either unweighted or weighted edges (see Fig. 2a–c). In this representation, node degree and node strength quantify information at the scale of single particles. One can compute the mean degree (a global property) of a network by calculating $$\langle{k}\rangle = \frac{1}{N}\sum_{i} k_{i},$$ (2.7) and one can similarly compute the mean strength with the formula $$\langle{s}\rangle = \frac{1}{N}\sum_{i} s_{i}.$$ (2.8) In an undirected network, mean degree is related to $$N$$ and the total number $$m$$ of edges through the relation $$\langle{k}\rangle = \frac{2m}{N}$$. It is sometimes useful to characterize a network using its degree distribution $$P(k)$$ [29], which gives the probability that a node has degree $$k$$. When one represents a granular packing as a contact network (see Fig. 2 and Section 3.1), which is a binary network (i.e., unweighted network), the degree $$k_i$$ of node $$i$$ is known more commonly in the physics literature as its contact number or coordination number$$Z_i$$. If every node has the same degree (or even if most nodes have this degree), such as in a regular lattice, one often refers to the mean coordination number $$Z$$ of the lattice or packing. It is well-known that coordination number is related to the stability of granular systems [61] and plays a critical role in the jamming transition [62–64], a change of phase from an underconstrained state to a rigid state that is characterized by the onset of mechanical stability.3 We discuss these ideas further throughout the review. 2.2.2 Walks and paths. In a network, a walk is an alternating sequence of nodes and edges that starts and ends at a node, such that consecutive edges are both incident to a common node. A walk thus describes a traversal from one node to another node (or to itself) along edges of a network. A path is a walk that does not intersect itself or visit the same node (or the same edge) more than once, except for a closed path, which starts and ends at the same node (see Section 2.2.3). One can compute the number of (unweighted) walks of a given length from a binary, undirected adjacency matrix $${\bf A}$$ [29]. The length $$l$$ of an unweighted walk is defined as the number of edges in the associated sequence (counting repeated edges the number of times that they appear). Letting $$\Xi^{l}_{ij}$$ denote the number of walks of length $$l$$ between nodes $$i$$ and $$j$$, one calculates $$\Xi^{l}_{ij} = [\mathbf{A}^{l}]_{ij}.$$ (2.9) Various types of random walks yield short paths between nodes in a network, and such ideas (and their relation to topics such as spectral graph theory) are very insightful for studying networks [65]. In an undirected network, a path from node $$i$$ to node $$j$$ is necessarily also a path from node $$j$$ to node $$i$$. However, this is not typically true in directed networks, which have sometimes been utilized in studies of granular force chains (see e.g., [13]). Depending on a network’s structure, there may be several or no paths between a given pair of nodes. An undirected network is called connected if there exists a path between each pair of nodes, and a directed network is called strongly connected if there is a path between each pair of nodes [29]. (A directed network is called weakly connected if its associated undirected network is connected, and a strongly connected network is necessarily also weakly connected.) In networks of granular packings, both the existence and lengths of paths can impact system behaviour. A connected network consists of a single component. When a network has multiple components, it is common to study one component at a time (e.g., focusing on a largest connected component (LCC), which is one that has the largest number of nodes). The length of an unweighted path is the number of edges in the associated sequence, and it is sometimes also called the hop distance (or occasionally, unfortunately, the topological distance). Paths in a network can also be weighted by defining some (possibly abstract) notion of distance associated with the edges of the network. For example, in a spatially-embedded network [37], distance may refer to actual physical distance along an edge, in which case the length of a weighted path in the network is given by the sum of the physical distances along the sequence of edges in the path. However, one can also consider ‘distance’ more abstractly. For example, in a transportation or flow network, one can define a distance between two adjacent nodes to be some measure of resistance between those nodes, and then the length of a weighted path in such a network is given by the sum of the resistances along the sequence of edges in the path.4 We use the term network distance to indicate a distance between two nodes (which can be either unweighted or weighted) that is computed by summing along edges in a path. A geodesic path—that is, a shortest path (which need not be unique)—between two nodes can be particularly relevant (though other short paths are often also important), and a breadth-first search (BFS) algorithm [67] is commonly employed to find geodesic paths in a network. The diameter of a graph is the maximum geodesic distance between any pair of nodes. Denoting the shortest, unweighted network distance (i.e., shortest-path distance) between nodes $$i$$ and $$j$$ as $$d_{ij}$$, the mean shortest-path distance $$L$$ between pairs of nodes in a graph is [68] $$L = \frac{1}{N(N-1)}\sum_{i,j \,\, (i \neq j)} d_{ij}.$$ (2.10) Note that one must be cautious when computing the mean shortest-path distance on disconnected networks (i.e., on networks that have more than one component), because the usual convention is to set the distance between two nodes in different components to be infinite [29]. Therefore, in a network with multiple components, the mean shortest-path distance $$L$$ from Eq. (2.10) is infinite. One solution to this problem is to compute $$L$$ for each component separately. Another network notion that relies on paths is network efficiency [69, 70] $$E = \frac{1}{N(N-1)}\sum_{i,j\,\, (i \neq j)}\frac{1}{d_{ij}}.$$ (2.11) One can generalize measures based on walks and paths to incorporate edge weights [70, 71]. Letting $$d^{w}_{ij}$$ denote the shortest, weighted network distance between nodes $$i$$ and $$j$$, one can define a weighted mean shortest-path distance $$L^{w}$$ and weighted efficiency $$E^{w}$$ as in Eqs. (2.10) and (2.11), respectively, but now one uses $$d_{ij}^{w}$$ instead of $$d_{ij}$$. The network efficiency $$E$$ is a normalized version of the Harary index of a graph [72]. Additionally, the convention $$d_{ij} = \infty$$ (or $$d^{w}_{ij} = \infty$$) if there is no path from $$i$$ to $$j$$ allows one to use Eq. (2.11) (or its weighted counterpart $$E^{w}$$) on connected graphs or on graphs with more than one component. For both unweighted and weighted scenarios, large values of network efficiency tend to correspond to small values of mean shortest-path length, and vice versa. One can also readily generalize notions of paths, distances, and efficiency to directed networks [69–71]. In later sections, we will describe the use of paths, walks and related ideas for investigating the structure of granular materials and their response to perturbations—including how these quantities change as a granular packing is compressed and goes through the jamming transition [73]—and we will also describe their use in specific applications, such as in understanding heat transfer through a granular material [24]. 2.2.3 Cycles. A cycle (i.e., a closed walk) in a network is a walk that begins and ends at the same node [30]. As with other walks, one way to characterize a cycle is by calculating its length. An $$l$$-cycle is a cycle in which $$l$$ edges are traversed (counting repeated edges the number of times that they appear in the cycle). A simple cycle is a cycle that does not include repeated nodes or edges, aside from one repetition of the origin node at the termination of a closed cycle. Thus, for example, a simple 3-cycle in an undirected network is a triangle. For the remainder of this review, we assume that cycles are simple cycles, unless we explicitly state otherwise. In the context of a granular packing, one can directly map particle contact loops—sets of physically-connected grains arranged in a circuit—to cycles in a corresponding graphical representation (see Fig. 2d). The length $$l$$ is odd for an odd cycle and even for an even cycle. We briefly note a few related concepts that are used to examine cycles in graphs because of their relevance to several network-based studies of granular materials. These are the notions of cycle space, cycle basis, and minimum cycle basis [30, 74]. The cycle space of an undirected graph is the set of all simple cycles in a graph along with all subgraphs that consist of unions of edge-disjoint simple cycles (i.e., they can share nodes but not edges) [75, 76]. A cycle basis is a minimal set of simple cycles such that any element of the cycle space can be written as a symmetric difference of cycles in the cycle basis [75]. Finally, for unweighted networks, a minimum cycle basis is a basis in which the total length of all cycles in the basis is minimal. For weighted networks, it is a basis in which the sum of the weights of all cycles in the basis is minimal. Minimum cycle bases can provide useful information about the structure and organization of cycles in a network, so several algorithms have been developed to extract them (see, for example, [77, 78]). Once one has determined a minimum cycle basis, one can examine the distribution of cycle lengths or define measures to quantify the participation of different nodes in cycles of different lengths. For example, Walker et al. [79, 80] defined the concept of a cycle-participation vector$$X_{i}^{\mathrm{cycle}} = [x_{i}^{0},x_{i}^{3},\dots,x_{i}^{l}]$$ for each node $$i$$. The elements of this vector count the number of cycles of each length in which node $$i$$ participates. In this definition, $$x_{i}^{3}$$ is the number of 3-cycles in which node $$i$$ participates, $$x_{i}^{4}$$ is the number of 4-cycles in which node $$i$$ participates, and so on (up to cycles of length $$l$$). If node $$i$$ is not part of any cycle, then $$x_{i}^{0} = 1$$ and $$x_{i}^{j} = 0$$ for all $$j \geq 3$$; otherwise, $$x_{i}^{0} = 0$$. One reason to examine cycles in granular networks [58, 73, 81–84] is that they can help characterize mesoscale structural features of a network. Cycles that are non-trivial involve more than a single node, but they do not typically embody global structures of a large network. This makes them appealing for studying network representations of granular materials, because mesoscale features seem to play a role in the behaviour of these systems [21]. Perhaps the most important motivation, however, is that cycles appear to be relevant for the stability and rigidity of a system. Specifically, in the context of structural rigidity theory, 3-cycles tend to be stabilizing structures that can maintain rigidity under applied forces [85], whereas 4-cycles can deform under applied forces (see Section 3.1.2). In Section 3, we discuss in more detail how cycles can help characterize granular systems. 2.2.4 Clustering coefficients. Clustering coefficients are commonly-used diagnostics to measure the density of triangles either locally or globally in a network [29]. For an unweighted, undirected network, the local clustering coefficient $$C_{i}$$ is usually defined as the number of triangles involving node $$i$$ divided by the number of triples centred at node $$i$$ [68, 86]. A triple is a set of three nodes that can include either three edges (to form a 3-cycle) or just two of them. In terms of the adjacency matrix and node degree, the local clustering coefficient is $$C_{i} = \frac{\sum_{hj} A_{hj} A_{ih} A_{ij}}{k_{i}(k_{i} - 1)}$$ (2.12) for $$k_i \geq 2$$ (and $$C_{i} = 0$$ if $$k_i \in \{0,1\}$$). One can then calculate a global clustering coefficient of a network as the mean of $$C_{i}$$ over all nodes: $$C = \frac{1}{N} \sum_{i}^{N} C_{i}.$$ (2.13) There is also another (and simpler) common way of defining a global clustering coefficient in a network that is particularly useful when trying to determine analytical approximations of expectations over ensembles of random graphs [29, 87]. The notion of a local clustering coefficient has also been extended to weighted networks in several ways [59, 88–90]. In one formulation [59], a local, weighted clustering coefficient $$C^{w}_{i}$$ is defined as $$C^{w}_{i} = \frac{1}{s_{i}(k_{i} - 1)} \sum_{j,h} \frac{(W_{ij} + W_{ih})} {2} A_{ij} A_{ih} A_{jh}$$ (2.14) for strength $$s_i > 0$$ and degree $$k_i \geq 2$$. The quantity $$C^{w}_{i} = 0$$ if either $$s_i = 0$$ (so that $$k_i = 0$$) or $$k_i = 1$$. Recall that $$\mathbf{W}$$ and $$\mathbf{A}$$ are, respectively, associated weighted and unweighted adjacency matrices. The mean of $$C^{w}_{i}$$ over all nodes gives a weighted clustering coefficient $$C^{w}$$ of a network. As we will discuss later (see Sections 3.1.3 and 3.2.4), clustering coefficients have been employed in several studies of granular materials. For example, they have been used to examine stability in granular networks [58, 73, 82–84]. See Fig. 3a for an example of the spatial distribution of a clustering coefficient in a granular packing. Fig. 3. View largeDownload slide Network diagnostics reveal different types of structure in granular packings. The two panels show the same compressed collection of photoelastic disks. One can measure the forces between particles and use the magnitude of the force between particle $$i$$ and particle $$j$$ to weight the edge between node $$i$$ and node $$j$$ in a weighted, undirected graph representation of the packing. (a) One can compute a weighted clustering coefficient at each network node (i.e., particle) to probe local structure in the packing. (b) One can compute a weighted betweenness centrality at each network node to probe how often its associated particle lies on a weighted shortest path between any two particles in the packing. We adapted this figure, with permission, from [21]. Fig. 3. View largeDownload slide Network diagnostics reveal different types of structure in granular packings. The two panels show the same compressed collection of photoelastic disks. One can measure the forces between particles and use the magnitude of the force between particle $$i$$ and particle $$j$$ to weight the edge between node $$i$$ and node $$j$$ in a weighted, undirected graph representation of the packing. (a) One can compute a weighted clustering coefficient at each network node (i.e., particle) to probe local structure in the packing. (b) One can compute a weighted betweenness centrality at each network node to probe how often its associated particle lies on a weighted shortest path between any two particles in the packing. We adapted this figure, with permission, from [21]. 2.2.5 Centrality measures. In network analysis, one often calculates centrality measures in efforts to quantify the importances of particular nodes, edges, or other structures in a network [29]. Different types of centralities characterize importance in different ways. The degree centrality (i.e., degree) of a node, for example, is simply the number of edges attached to it (see Section 2.2.1). A few other types of centrality that have been used to study granular materials are closeness centrality, node betweenness centrality, edge betweenness centrality, and subgraph centrality. Notions of closeness centrality of a node measure how close that node is to other nodes in a network [91]. For a given node $$i$$, the most common notion of closeness is defined as the inverse of the sum over the shortest-path lengths from node $$i$$ to all other nodes $$j$$ in a network. That is, node $$i$$’s closeness centrality is $$H_{i} = \frac{N-1}{\sum_{j \neq i} d_{ij}}.$$ (2.15) Note that if we use the convention that the distance between two nodes in different components is infinite, then Eq. (2.15) only makes sense for connected networks. For any network with more than one component, Eq. (2.15) yields a closeness centrality of $$0$$. The geodesic node betweenness centrality of node $$i$$ is the fraction of geodesic paths (either unweighted or weighted) between distinct nodes (not including $$i$$) that traverse node $$i$$ [91]. Let $$\psi_{gh}(i)$$ denote the number of geodesic paths from node $$g$$ to node $$h$$ that traverse node $$i$$ (with $$i \not \in \{g,h\}$$), and let $$\psi_{gh}$$ denote the total number of geodesic paths from node $$g$$ to node $$h$$. The geodesic node betweenness centrality of node $$i$$ is then $$B_{i} = \sum_{g,h;\,g \neq h}\frac{\psi_{gh}(i)}{\psi_{gh}}, \quad i \not\in \{g,h\}.$$ (2.16) Geodesic node betweenness can probe the heterogeneity of force patterns in granular networks. See Fig. 3b for an example spatial distribution of a geodesic node betweenness centrality in an experimental granular packing. One can also compute a geodesic edge betweenness centrality by calculating the fraction of shortest paths (either unweighted or weighted) that traverse it [92]. Let $$\psi_{gh}(i,j)$$ denote the number of geodesic paths from node $$g$$ to node $$h$$ that traverse the edge that is attached to node $$i$$ and node $$j$$, and let $$\psi_{gh}$$ denote the total number of geodesic paths from node $$g$$ to node $$h$$. The geodesic edge betweenness centrality of this edge is then $$B^{e}_{ij} = \sum_{g,h;\, g \neq h} \frac{\psi_{gh}(i,j)}{\psi_{gh}}.$$ (2.17) Another measure of node importance is subgraph centrality$$Y$$ [93, 94], which quantifies a node’s participation in closed walks of all lengths. Recall from Section 2.2.2 that one can write the number of length-$$l$$ walks from node $$i$$ to node $$j$$ in terms of powers of the adjacency matrix $$\mathbf{A}$$. To calculate closed walks of length $$l$$ that begin and end at node $$i$$, we take $$i = j$$ in Eq. (2.9). The subgraph centrality of node $$i$$, with a specific choice for how much we downweight longer paths, is then given by $$Y_{i} = \sum_{l=0}^{\infty} \frac{[\mathbf{A}^{l}]_{ii}}{l!}.$$ (2.18) Because shorter walks are weighted more strongly than longer walks in Eq. (2.18), they contribute more to the value of subgraph centrality (In other contexts, centrality measures based on walks have also been used to compare the spatial efficiencies of different networks, and such ideas are worth exploring in granular materials [72].). One can also express subgraph centrality in terms of the eigenvalues and eigenvectors of the adjacency matrix [93]. Let $$v^{i}_{\alpha}$$ denote the $$i{\mathrm{th}}$$ component of the $$\alpha{\mathrm{th}}$$ eigenvector $$\boldsymbol{v}_{\alpha}$$ of $$\mathbf{A}$$, and let $$\lambda_\alpha$$ denote the corresponding $$\alpha{\mathrm{th}}$$ eigenvalue. One can then write $$Y_{i} = \sum_{\alpha=1}^{N} (v^{i}_{\alpha})^{2} e^{\lambda_\alpha}.$$ (2.19) One can then calculate a mean subgraph centrality $$Y$$ by averaging $$Y_i$$ over the nodes in a network. In one study of granular materials [58], a subgraph centrality was examined for weighted networks by considering the eigenvalues and eigenvectors of the weight matrix $$\mathbf{W}$$ in Eq. (2.19). One can also compute network bipartivity$$R$$ [95] to quantify the contribution to mean subgraph centrality $$Y$$ from closed walks of even length. In particular, the network bipartivity $$R_{i}$$ of node $$i$$ is $$R_{i} = \frac{{Y}^{\mathrm{even}}_{i}}{Y_{i}},$$ (2.20) where $$Y^{\mathrm{even}}_{i}$$ is the contribution to the sum in Eq. (2.18) from even values of $$l$$ (i.e., even-length closed walks). As with other node diagnostics, one can average bipartivity over all nodes in a network to obtain a global measure, which we denote by $$R$$. In Section 3, we will discuss calculations of closeness, betweenness, and subgraph centralities in granular packings. Obviously, our discussion above does not give an exhaustive presentation of centrality measures, and other types of centralities have also been used in studies of granular materials (see for example, [21]). 2.2.6 Subgraphs, motifs and superfamilies. One can interpret the local clustering coefficient in Eq. (2.12) as a relationship between two small subgraphs: a triangle and a connected triple. Recall that a subgraph of a graph $$G$$ is a graph constructed using a subset of $$G$$’s nodes and edges. Conceptually, one can interpret small subgraphs as building blocks or subunits that together can be used to construct a network. For example, in a directed network, there exist three possible 2-node subgraphs (i.e., dyads): the dyad in which node $$i$$ is adjacent to node $$j$$ by a directed edge, the dyad in which node $$j$$ is adjacent to node $$i$$ by a directed edge, and the dyad in which both of these adjacencies exist. In a directed, unweighted graph, there are 13 different connected 3-node subgraphs [96] (see Fig. 4). Fig. 4. View largeDownload slide Subgraphs in networks. We show all 13 different 3-node subgraphs that can occur in a directed, unweighted graph. A motif is a subgraph that occurs often (typically relative to some null model) in a particular network or set of networks [97]. We reproduced this figure, with permission, from [96]. Fig. 4. View largeDownload slide Subgraphs in networks. We show all 13 different 3-node subgraphs that can occur in a directed, unweighted graph. A motif is a subgraph that occurs often (typically relative to some null model) in a particular network or set of networks [97]. We reproduced this figure, with permission, from [96]. The term motif is sometimes used for a small subgraph that occurs often in a particular network or set of networks (typically relative to some null model, such as a randomly rewired network that preserves the original degree distribution) [96–99]. Borrowing terminology from genetics, these motifs appear to be overexpressed in a network (or set of networks). Unsurprisingly, the number of $$n$$-node subgraphs increases very steeply with $$n$$, so identifying subgraphs in large networks is computationally expensive, and many algorithms have been developed to estimate the number of subgraphs in an efficient (though approximate) way; see, for example, [96, 100–104]. In applying algorithms for motif counting to data, one seeks to identify subgraphs that are present more often than expected in some appropriate random-network null model. The over-representation of a motif in a network is often interpreted as indicative of its playing a role in the function of that network (though one has to be cautious about drawing such conclusions). For example, 3-node motifs can form feedforward loops in which there are directed edges from node $$i_1$$ to node $$i_2$$, from node $$i_2$$ to node $$i_3$$ and from node $$i_3$$ to node $$i_1$$. The identification and characterization of motifs has yielded insights into the structure and function of a variety of systems, including food webs [105], gene-regulation networks of yeast [96], neuronal networks of the macaque monkey [106], and others. For different types of networks, one can also identify so-called superfamilies, which are sets of networks that have similar motif-frequency distributions [98]. There also exists a less-stringent definition of a superfamily; in this definition, one disregards whether a subgraph is a motif in the sense of it being more abundant than expected from some random-graph null model and instead considers a superfamily to be a set of networks that have the same rank-ordering of the number of $$n$$-node subgraphs for some fixed value of $$n$$ [107]. In either case, one can examine different superfamilies to help understand the role that specific motifs (or subgraphs) or sets of motifs (or subgraphs) may have in potentially similar functions of networks in a given superfamily. Subgraphs, motifs and superfamilies have been examined in several studies that applied network analysis to granular materials [80, 108–110]. They have revealed interesting insights into the deformation and reconfiguration that occurs in granular systems for different types of loading conditions and external perturbations. We discuss these ideas further in Sections 3.1.4 and 3.3.3. 2.2.7 Community structure. Many real-world networks have structure on intermediate scales (mesoscales) that can arise from particular organizations of nodes and edges [31–34, 111]. The most commonly-studied mesoscale network property is community structure [31, 32], which describes sets of nodes, called communities, that are densely (or strongly) interconnected to each other but only weakly connected to other dense sets of nodes. In other words, a community has many edges (or large total edge weight, in the case of weighted networks) between its own nodes, but the number and/or weight of edges between nodes in different communities is supposed to be small. Once one has detected communities in a network, one way to quantify their organization is to compute and/or average various network quantities over the nodes (or edges) within each community separately, rather than over the entire network. For example, one can compute the size (i.e., number of nodes) of a community, mean path lengths between nodes in a given community, or some other quantity to help characterize the architecture of different communities in a network. Studying community structure can reveal useful insights about granular systems, whose behaviour appears to be influenced by mesoscale network features [21, 27, 112–117]. Community structure and methods for detecting communities have been studied very extensively [31]. We will briefly discuss the method of modularity maximization [29, 118, 119], in which one optimizes an (occasionally infamous) objective function known as modularity, as this approach has been employed previously in several studies of granular materials (see, e.g., Section 3.2.2). However, myriad other approaches exist for studying community structure in networks. These include stochastic block models (SBMs) and other methods for statistical inference (which are increasingly favored by many scholars) [111], approaches based on random walks (e.g., InfoMap [120]), various methods for detecting local community structure (see, e.g., [121, 122]), edge-based communities [123], and many others. The goal of modularity maximization is to identify communities of nodes that are more densely (or more strongly) interconnected with other nodes in the same community than expected with respect to some null model. To do this, one maximizes a modularity objective function $$Q = \sum_{i,j} [W_{ij} - \gamma P_{ij}] \delta(g_{i},g_{j}),$$ (2.21) where $$g_{i}$$ is the community assignment of node $$i$$ and $$g_{j}$$ is the community assignment of node $$j$$, and where the Kronecker delta $$\delta(g_{i},g_{j})=1$$ if $$g_{i} = g_{j}$$ and $$\delta(g_{i},g_{j})=0$$ otherwise. The quantity $$\gamma$$ is a resolution parameter that adjusts the relative average sizes of communities [124, 125], where smaller values of $$\gamma$$ favor larger communities and larger values of $$\gamma$$ favor smaller communities [126]. The element $$P_{ij}$$ is the expected weight of the edge between node $$i$$ and node $$j$$ under a specified null model. In many contexts, the most common choice is to determine the null-model matrix elements $$P_{ij}$$ from the Newman–Girvan (NG) null model [118, 127, 128], for which $$P^{\mathrm{NG}}_{ij} = \frac{s_{i} s_{j}}{2m},$$ (2.22) where $$s_{i}=\sum_j W_{ij}$$ is the strength (and, for unweighted networks, the degree $$k_i$$) of node $$i$$ and $$m=\frac{1}{2}\sum_{i,j} W_{ij}$$ is the total edge weight (and, for unweighted networks, the total number of edges) in the network. There are several other null models, which are usually based on a random-graph model, and they can incorporate system features (such as spatial information) in various ways [129]. In the next part of this subsubsection, we discuss a physically-motivated null model that is particularly useful for studying granular force networks. Maximizing $$Q$$ is NP-hard [130], so it is necessary to use computational heuristics to identify near-optimal partitions of a network into communities of nodes [124]. Two well-known choices are the Louvain [131] and the Louvain-like [132] locally greedy algorithms, which begin by placing all nodes in their own community, and they then iteratively agglomerate nodes when the resulting partition increases modularity. Because of the extreme near degeneracy of the modularity landscape (a very large number of different partitions can have rather similar values of the scalar $$Q$$), it is often useful to run such an algorithm many times to construct an ensemble of partitions, over which one can average various properties to yield a ‘consensus’ description of community structure [126, 129, 133, 134]. Physical considerations. Community-detection tools, such as modularity maximization, have often been applied to social, biological, and other networks [31, 32]. In applying these techniques to granular materials, however, it is important to keep in mind that the organization of particulate systems (such as the arrangements of particles and forces in a material) is subject to significant spatial and physical constraints, which can severely impact the types of organization that can arise in a corresponding network representation of the material. When studying networks that are embedded in real space or constructed via some kind of physical relation between elements, it is often crucial to consider the spatial constraints—and, more generally, a system’s underlying physics—and their effects on network architecture [37]. Such considerations also impact how one should interpret network diagnostics such as path lengths and centrality measures, the null models that one uses in procedures such as modularity maximization, and so on. The NG null model was constructed to be appropriate for networks in which a connection between any pair of nodes is possible. Clearly, in granular materials—as in other spatially-embedded systems [37]—this assumption is unphysical and therefore problematic. Bassett et al. [126] defined a null model that accounts explicitly for geographical (and hence spatial) constraints in granular materials, in which each particle can contact only its nearest neighbours [27]. In the context of granular networks with nodes representing particles and edges representing forces between those particles, the geographical null model$${\bf P}$$ in [126] has matrix elements $$P_{ij} = \rho A_{ij} ,$$ (2.23) where $$\rho$$ is the mean edge weight in the network and $${\bf A}$$ is the binary adjacency matrix of the network. In this particular application, $$\rho = \overline{f} := \langle\, f_{ij} \rangle$$ is the mean inter-particle force. As we illustrate in Fig. 5, modularity maximization with the geographical null model [Eq. (2.23)] produces different communities than modularity maximization with the NG null model [Eq. (2.22)] [27, 112, 113]. Fig. 5. View largeDownload slide Modularity maximization with different null models can reveal distinct types of community structures in granular networks. In a network representation in which particles are represented by nodes and contact forces are represented by weighted edges, (a) using the NG null model helps uncover contiguous domains in a granular system, and (b) using a geographical null model on the same network helps detect chain-like structures that are reminiscent of force chains. In both panels, nodes (i.e., particles) of the same colour are assigned to the same community. Fig. 5. View largeDownload slide Modularity maximization with different null models can reveal distinct types of community structures in granular networks. In a network representation in which particles are represented by nodes and contact forces are represented by weighted edges, (a) using the NG null model helps uncover contiguous domains in a granular system, and (b) using a geographical null model on the same network helps detect chain-like structures that are reminiscent of force chains. In both panels, nodes (i.e., particles) of the same colour are assigned to the same community. Generalization of modularity maximization to multilayer networks. Although studying community structure in a given granular packing can provide important insights, one is also typically interested in how such mesoscale structures reconfigure as a material experiences external perturbations, such as those from applied compression or sheer. To examine these types of questions, one can optimize a multilayer generalization of modularity to study multilayer granular force networks in which each layer represents a network at a different step in the evolution of the system (for example, at different time steps or at different packing fractions) [113]. In Fig. 6, we show a schematic of a multilayer construction that has been employed in such investigations. See [36, 135] for reviews of multilayer networks (including generalizations of this construction). Fig. 6. View largeDownload slide A schematic of a multilayer network with layer-dependent community structure. In this example, each layer represents a static granular force network in which nodes (i.e., particles) are adjacent to one another via intralayer weighted edges (e.g., representing contact forces). Additionally, the same particle in consecutive layers is adjacent to itself via an interlayer edge of uniform weight $$\omega$$. For clarity, we only show two such couplings, but these interlayer edges exist for all particles (and between all consecutive layers). One can extract communities that change across layers—for example, if layers represent time, these are time-evolving communities—to study mesoscale organization in a granular system and to help understand how it reconfigures due to external loading (such as compression). In this schematic, we use different colours to label the particles that belong to different communities. Note that the same community can persist across several (or all) layers and reconfigure in terms of its constituent particles and the mean strength of its nodes. We reproduced this figure, with permission, from [113]. Fig. 6. View largeDownload slide A schematic of a multilayer network with layer-dependent community structure. In this example, each layer represents a static granular force network in which nodes (i.e., particles) are adjacent to one another via intralayer weighted edges (e.g., representing contact forces). Additionally, the same particle in consecutive layers is adjacent to itself via an interlayer edge of uniform weight $$\omega$$. For clarity, we only show two such couplings, but these interlayer edges exist for all particles (and between all consecutive layers). One can extract communities that change across layers—for example, if layers represent time, these are time-evolving communities—to study mesoscale organization in a granular system and to help understand how it reconfigures due to external loading (such as compression). In this schematic, we use different colours to label the particles that belong to different communities. Note that the same community can persist across several (or all) layers and reconfigure in terms of its constituent particles and the mean strength of its nodes. We reproduced this figure, with permission, from [113]. One way to detect multilayer communities in a network is to use a generalization of modularity maximization [136], which was derived for multilayer networks with interlayer edges between nodes in different layers that represent the same entity. For simplicity, we suppose that all edges are bidirectional. One maximizes $$Q_{\mathrm{multi}} = \frac{1}{2 \eta} \sum_{ijqr} [(\mathscr{W}_{ijq} - \gamma_{q}\mathscr{P}_{ijq}) \delta_{qr} + \omega_{jqr}\delta_{ij}] \delta(g_{iq},g_{jr}),$$ (2.24) where $$\mathscr{W}_{ijq}$$ is the $$(i,j){\mathrm{th}}$$ component of the $$q{\mathrm{th}}$$ layer of the adjacency tensor $$\mathscr{W}$$ [137] associated with the multilayer network, $$\mathscr{P}_{ijq}$$ is the $$(i,j){\mathrm{th}}$$ component of the $$q{\mathrm{th}}$$ layer of the null-model tensor, $$\gamma_{q}$$ is a resolution parameter (sometimes called a structural resolution parameter) for layer $$q$$, and $$\omega_{jqr}$$ is the interlayer coupling between layers $$q$$ and $$r$$. (In the context of multilayer representations of temporal networks, if $$\mathcal{\omega}_{jqr} = \omega$$ for all $$j$$, $$q$$ and $$r$$, one can interpret $$\omega$$ as a temporal resolution parameter.) More specifically, $$\mathcal{\omega}_{jqr}$$ is the strength of the coupling that links node $$j$$ in layer $$q$$ to itself in layer $$r$$. (This type of interlayer edge, which occurs between nodes in different layers that represent the same entity, is called a diagonal edge [36].) The quantities $$g_{iq}$$ and $$g_{jr}$$, respectively, are the community assignments of node $$i$$ in layer $$q$$ and node $$j$$ in layer $$r$$. The intralayer strength of node $$j$$ in layer $$q$$ is $$s_{jq} = \sum_{i} \mathscr{W}_{ijq}$$, and the interlayer strength of node $$j$$ in layer $$q$$ is $$\zeta_{jq} = \sum_{r}{\mathcal{\omega}_{jqr}}$$, so the multilayer strength of node $$j$$ in layer $$q$$ is given by $$\kappa_{jq} = s_{jq} + \zeta_{jq}$$. Finally, the normalization factor $$\eta = \frac{1}{2} \sum_{jq} \kappa_{jq}$$ is the total strength of the adjacency tensor.5 Maximizing multilayer modularity [Eq. (2.24)] allows one to examine phenomena such as evolving communities in time-dependent networks or communities that evolve with respect to some other parameter, and communities in networks with multiple types of edges. Capturing such behaviour has been useful in many applications, including financial markets [128], voting patterns [136], international relations [139], international migration [140], disease spreading [129], human brain dynamics [141–143] and more. In the context of granular matter, multilayer community detection allows one to examine changes in community structure of a force network, in which communities can either persist or reconfigure, with respect to both particle content and the mean strength of nodes inside a community, due to applied loads on a system. 2.2.8 Flow networks. One can examine many natural and engineered systems—such as animal and plant vasculature, fungi, and urban transportation networks [144–150]—from the perspective of flow networks (which are often directed) that transport a load (of fluids, vehicles, and so on) along their edges. It is of considerable interest to examine how to optimize flow through a network [29, 151]. A well-known result from optimization theory is the maximum-flow–minimum-cut theorem [29, 151, 152]: for a suitable notion of flow and under suitable assumptions, the maximum flow that can pass from a source node to a sink node is given by the total weight of the edges in the minimum cut, which is the set of edges with smallest total weight that, when removed, disconnect the source and the sink. A related notion, which applies to networks in which there is some cost associated with transport along network edges, is that of maximum-flow–minimum-cost. In this context, one attempts to find a route through a network that maximizes flow transmission from source to sink, while minimizing the cost of flow along network edges [151, 152]. The maximum-flow–minimum-cut and maximum-flow–minimum-cost problems are usually examined under certain constraints, such as flow conservation at each node and an upper bound (e.g., limited by a capacitance) on flow through any edge. One can examine granular networks from such a perspective by considering a flow of force along a network formed by contacting grains. We discuss relevant studies in Section 3.3.1. 2.2.9 Connected components and percolation. Sometimes it is possible to break a network into connected subgraphs called components (which we introduced briefly in Section 2.2.2). A component, which is sometimes called a cluster, is a subgraph $$G_C$$ of a graph $$G$$ such that at least one path exists between each pair of nodes in $$G_C$$ [29]. Components are maximal subsets in the sense that the addition of another node of $$G$$ to a component destroys the property of connectedness. An undirected graph is connected when it consists of a single component. Networks with more than one component often have one component that has many more nodes than the others, so there can be one large component and many small components. One can find the components of a graph using a breadth-first search (BFS) algorithm [67], and one can determine the number of components by counting the number of $$0$$ eigenvalues of a graph’s combinatorial Laplacian matrix [29]. To study graph components, one can also use methods from computational algebraic topology. Specifically, the zeroth Betti number $$\beta_{0}$$ indicates the number of connected components in a graph [153] (see Section 2.2.10). Percolation theory [29, 154–156], which builds on ideas from subjects such as statistical physics and probability theory, is often used to understand the emergence and behaviour of connected components in a graph [52]. For example, in the traditional version of what is known as bond percolation (which is also traditionally studied on a lattice rather than on a more general network) [157], edges are occupied with probability $$p$$, and one examines quantities such as the size distributions of connected components as one varies the parameter $$p$$, called the bond occupation probability. It is especially interesting to determine a critical value $$p_c$$, called the percolation threshold, at which there is a phase transition: below $$p_{c}$$, there is no percolating component (or cluster), which spans the system and connects opposite sides; above $$p_{c}$$, there is such a cluster [157, 158]. In the latter scenario, it is common to say that there is a ‘percolating network’. In percolation on more general networks, one can study how the size of the largest component, as a fraction of the size of a network, changes with $$p$$. Related ideas also arise in the study of components in Erdős–Rényi random graphs $$G(N,p)$$, in which one considers an ensemble of $$N$$-node graphs and $$p$$ is the independent probability that an edge exists between any two nodes [29, 30, 159, 160]. In the limit $$N \rightarrow \infty$$, the size of the LCC undergoes a phase transition at a critical probability $$p_{c} = 1/N$$. When $$p < p_{c}$$, the ER graph in expectation does not have a giant connected component (GCC); at $$p = p_{c}$$, a GCC emerges whose size scales linearly with $$N$$ for $$p > p_{c}$$. Similarly, for bond percolation on networks, a transition occurs at a critical threshold $$p_c$$, such that for $$p > p_c$$, there is a GCC (sometimes also called a ‘giant cluster’ or ‘percolating cluster’) whose size is a finite fraction of the total number $$N$$ of nodes as $$N \rightarrow \infty$$ [29, 52]. When studying percolation on networks, quantities of interest include the fraction of nodes in the LCC, the mean component size, the component-size distribution and critical exponents that govern how these quantities behave just above the percolation threshold [29, 158, 161]. We will see in Section 3 that it can be informative to use ideas from percolation theory to study the organization of granular networks. For example, it is particularly interesting to examine how quantities such as the number and size of connected components evolve with respect to packing density (or another experimental parameter). [28, 73, 162–169]. Some studies have considered connectivity percolation transitions, which are characterized by the appearance of a connected component that spans a system (i.e., a percolating cluster, as reflected by an associated GCC in the infinite-size limit of a network); or rigidity percolation transitions, which can be used to examine the transition to jamming [170–177]. Rigidity percolation is similar to ordinary bond percolation (which is sometimes used to study connectivity percolation), except that edges represent the presence of rigid bonds between network nodes [178, 179] and one examines the emergence of rigid clusters in the system with respect to the fraction of occupied bonds. One can also study percolation in force networks by investigating the formation of connected components and the emergence of a percolating cluster of contacts as one varies a force threshold, which is a threshold applied to a force-weighted adjacency matrix (representing contact forces between particles) to convert it to a binary adjacency matrix [73, 163–169, 176, 180–182]. (See Section 3.2.4 for additional discussion.) However, it is important to note that when studying networks of finite size, one needs to be careful with claims about GCCs and percolation phase transitions, which are well-defined mathematically only in the limit of infinite system size. 2.2.10 Methods from algebraic topology and computational topology. The tools that we have described thus far rely on the notion of a dyad (i.e., a 2-node subgraph) as the fundamental unit of interest (see Fig. 7a). However, recent work in algebraic topology and computational topology [153, 183–185] offers a complementary view, in which the fundamental building blocks that encode relationships between elements of a system are $$k$$-simplices (each composed of $$k+1$$ nodes), rather than simply nodes and dyadic relations between them (see Fig. 7b). These structures can encode ‘higher-order’ interactions and can be very useful for understanding the architecture and function of real-world networks (e.g., they yield a complementary way to examine mesoscale network features), and they have been insightful in the study of sensor networks [186], contagion spreading [187], protein interactions [188], neuronal networks [51, 189], and many other problems. See [190, 191] for further discussion and pointers to additional applications. The discussion in [192] is also useful. Fig. 7. View largeDownload slide Algebraic topology and clique complexes. (a) An interaction in a graph is a dyad (i.e., a 2-node subgraph) or a node and a self-edge. (b) An alternative fundamental unit is a $$k$$-simplex. A 0-simplex is a node, a 1-simplex is an edge, a 2-simplex is a filled triangle, and so on. (c) A collection of $$k$$-simplices is called a simplicial complex, and one type of simplicial complex that can be used to encode the information in a graph is a clique complex (sometimes also called a flag complex). One constructs a clique complex by taking every $$k$$-clique (a complete subgraph of $$k$$ nodes) in a graph $$G$$ to be a simplex of the same number of nodes. (d) An interesting feature that can occur in a simplicial complex is a cycle, which is a closed arrangement of a collection of $$k$$-simplices. The purple edges in the upper object indicate a one-dimensional cycle that encloses a region filled in by simplices, whereas the purple edges in the lower object indicate a one-dimensional cycle that encloses a hole. (e) One can use a filtration to decompose a weighted graph into a sequence of binary graphs. For example, if one uses edge weight as a filtration parameter, one can represent a weighted graph as a sequence of unweighted graphs, which in turn yields a sequence of unweighted clique complexes. We adapted this figure, with permission, from [57]. Fig. 7. View largeDownload slide Algebraic topology and clique complexes. (a) An interaction in a graph is a dyad (i.e., a 2-node subgraph) or a node and a self-edge. (b) An alternative fundamental unit is a $$k$$-simplex. A 0-simplex is a node, a 1-simplex is an edge, a 2-simplex is a filled triangle, and so on. (c) A collection of $$k$$-simplices is called a simplicial complex, and one type of simplicial complex that can be used to encode the information in a graph is a clique complex (sometimes also called a flag complex). One constructs a clique complex by taking every $$k$$-clique (a complete subgraph of $$k$$ nodes) in a graph $$G$$ to be a simplex of the same number of nodes. (d) An interesting feature that can occur in a simplicial complex is a cycle, which is a closed arrangement of a collection of $$k$$-simplices. The purple edges in the upper object indicate a one-dimensional cycle that encloses a region filled in by simplices, whereas the purple edges in the lower object indicate a one-dimensional cycle that encloses a hole. (e) One can use a filtration to decompose a weighted graph into a sequence of binary graphs. For example, if one uses edge weight as a filtration parameter, one can represent a weighted graph as a sequence of unweighted graphs, which in turn yields a sequence of unweighted clique complexes. We adapted this figure, with permission, from [57]. A collection of simplices that are joined in a compatible way is called a simplicial complex, which is a generalization of a graph that can encode non-dyadic relations [185]. More precisely, and following [51], we define an (abstract) simplicial complex$$\mathscr{X}$$ as a pair of sets: $$V_{\mathscr{X}}$$, called the vertices (or nodes); and $$S_{\mathscr{X}}$$, called the simplices, each of which is a finite subset of $$V_{\mathscr{X}}$$, subject to the requirement that if $$\sigma \in S_{\mathscr{X}}$$, then every subset $$\tau$$ of $$\sigma$$ is also an element of $$S_{\mathscr{X}}$$. A simplex with $$k$$ elements is called a $$(k-1)$$-simplex, and subsets $$\tau \subset \sigma$$ are called faces of $$\sigma$$. Using this notation, a $$0$$-simplex is a node, a $$1$$-simplex is an edge and its two incident nodes (i.e., a dyad), a 2-simplex is a filled triangle, and so on (see Fig. 7b). One type of simplicial complex that can be used to encode the information in a graph is a clique complex (sometimes also called a flag complex); we show an example in Fig. 8. To construct the clique complex of a graph $$G$$, one associates every $$k$$-clique (a complete—i.e., fully connected—subgraph of $$k$$ nodes) in $$G$$ with a $$(k-1)$$-simplex. One can thus think of building the clique complex of a graph $$G$$ as ‘filling in’ all of the $$k$$-cliques in $$G$$ (see Fig. 7c). Note that we use the terms $$k$$-simplex and $$k$$-clique because they are standard, but it is important not to confuse the use of $$k$$ in this context with the use of $$k$$ as the (also standard) notation for node degree. Fig. 8. View largeDownload slide An example force network, an associated filtration over the flag complex (i.e., clique complex), and PDs. (a) In the force network, coloured edges represent the magnitude of the force between contacting particles, which are represented as nodes in the network. In order from smallest to largest, the four values of the force are $$\theta_{1}$$ (dark blue), $$\theta_{2}$$, (cyan), $$\theta_{3}$$ (green), and $$\theta_{4}$$ (red). (b) The flag complex is formed by filling in all 3-particle loops (i.e., triangular loops) with the smallest value of the force along any of its edges. Defining a filtration over the flag complex avoids counting these 3-particle loops. (c–f) The sequence of complexes corresponding to the filtration over the flag complex; one obtains the sequence by descending the four levels of the force threshold $$\theta$$. (g) The $$\beta_{0}$$ persistence diagram $$\mathrm{PD}_{0}$$. (h) The $$\beta_{1}$$ persistence diagram $$\mathrm{PD}_{1}$$. We adapted this figure, with permission, from [165]. Fig. 8. View largeDownload slide An example force network, an associated filtration over the flag complex (i.e., clique complex), and PDs. (a) In the force network, coloured edges represent the magnitude of the force between contacting particles, which are represented as nodes in the network. In order from smallest to largest, the four values of the force are $$\theta_{1}$$ (dark blue), $$\theta_{2}$$, (cyan), $$\theta_{3}$$ (green), and $$\theta_{4}$$ (red). (b) The flag complex is formed by filling in all 3-particle loops (i.e., triangular loops) with the smallest value of the force along any of its edges. Defining a filtration over the flag complex avoids counting these 3-particle loops. (c–f) The sequence of complexes corresponding to the filtration over the flag complex; one obtains the sequence by descending the four levels of the force threshold $$\theta$$. (g) The $$\beta_{0}$$ persistence diagram $$\mathrm{PD}_{0}$$. (h) The $$\beta_{1}$$ persistence diagram $$\mathrm{PD}_{1}$$. We adapted this figure, with permission, from [165]. One important feature of a simplicial complex is the potential presence of cycles.6 A cycle can consist of any number of nodes, and a $$k$$-dimensional cycle is defined as a closed arrangement of $$k$$-simplices, such that a cycle has an empty boundary7. For example, Fig. 7d illustrates a closed arrangement of $$1$$-simplices (i.e., edges) that forms a one-dimensional cycle. It important to distinguish between cycles that encircle a region that is filled by simplices with cycles that enclose a void (which is often called a ‘hole’ for the case of one-dimensional cycles). For example, the set of purple edges in the object in the upper portion of Fig. 7d constitutes a one-dimensional cycle that surrounds a region filled by 2-simplices (i.e., filled triangles), whereas the set of purple edges in the object in the lower portion of Fig. 7d constitutes a one-dimensional cycle that encloses a hole. Characterizing the location and prevalence of void-enclosing cycles in the clique complex of a network representation of a granular packing can offer fascinating insights into the packing’s structure [163]. One way to do this is by computing topological invariants such as Betti numbers [153, 165, 183]. The $$k{\mathrm{th}}$$ Betti number $$\beta_k$$ counts the number of inequivalent $$k$$-dimensional cycles that enclose a void, where two $$k$$-dimensional cycles are equivalent if they differ by a boundary of a collection of $$(k+1)$$-simplices. In other words, the $$k{\mathrm{th}}$$ Betti number $$\beta_k$$ counts the number of non-trivial equivalence classes of $$k$$-dimensional cycles and can thus also be interpreted as counting the number of voids (i.e., ‘holes’ of dimension $$k$$).8 The zeroth Betti number $$\beta_{0}$$ gives the number of connected components in a network, the first Betti number $$\beta_{1}$$ gives the number of inequivalent one-dimensional cycles that enclose a void (i.e., it indicates loops), the second Betti number $$\beta_{2}$$ gives the number of inequivalent 2D cycles that enclose a void (i.e., it indicates cavities), and so on. Another useful way to examine the topological features that are determined by equivalence classes of $$k$$-dimensional cycles (i.e., components, loops, cavities, and so on) is to compute persistent homology (PH) of a network. For example, to compute PH for a weighted graph, one can first decompose it into a sequence of binary graphs. One way to do this is to begin with the empty graph and add one edge at a time in order of decreasing edge weights (see Fig. 7e). More formally and following [189], this process can translate information about edge weights into a sequence of binary graphs as an example of what is called a filtration [185, 190]. The sequence $$G_0 \subset G_1 \subset \dots \subset G_{|\mathscr{E}|}$$ of unweighted graphs begins with the empty graph $$G_0$$, and one adds a single edge at a time (or multiple edges, if some edges have the same weight) in order from largest edge weight to smallest edge weight. (One can also construct filtrations in other ways.) Constructing a sequence of unweighted graphs in turn yields a sequence of clique complexes [195], allowing one to examine equivalence classes of cycles as one varies the edge weight $$\theta$$ (or another filtration parameter). Important values of $$\theta$$ include the weight $$\theta_{\mathrm{birth}}$$ associated with the first graph in which an equivalence class (i.e., a topological feature) occurs (i.e., its birth coordinate) and the edge weight $$\theta_{\mathrm{death}}$$ associated with the first graph in which the feature disappears (i.e., its death coordinate), such as by being filled in with higher-dimensional simplices or by merging with an older feature. One potential marker of the relative importance of a particular feature (a component, a loop, and so on) in the clique complex is how long it persists, as quantified by its lifetime$$\theta_{\mathrm{birth}} - \theta_{\mathrm{death}}$$ (although short-lived features can also be meaningful [190, 192]). A large lifetime indicates robust features that persist over many values of a filtration parameter. Persistence diagrams (PDs) are one useful way to visualize the evolution of $$k$$-dimensional cycles with respect to a filtration parameter. PDs encode birth and death coordinates of features as a collection of persistence points$$(\theta_{\mathrm{birth}},\theta_{\mathrm{death}})$$ in a planar region. One can construct a PD for each Betti number: a $$\beta_{0}$$ PD (denoted by $$\mathrm{PD}_{0}$$) encodes the birth and death of components in a network, a $$\beta_{1}$$ PD (denoted by $$\mathrm{PD}_{1}$$) encodes the birth and death of loops and so on. To demonstrate some key aspects of a filtration, the birth and death of topological features, and PDs, we borrow and adapt an example from [165]. Consider the small granular force network in Fig. 8a; the nodes represent particles in a 2D granular packing, and the coloured edges represent the magnitude of the inter-particle forces (of which there are four distinct values) between contacting particles. In a 2D system like this one, the only relevant Betti numbers are $$\beta_{0}$$ and $$\beta_{1}$$, as all others are $$0$$. In Fig. 8b, we show the flag complex (which is essentially the same as a clique complex [190]) of the granular network, where the colour of a triangle indicates the value corresponding to the minimum force along any of its edges. Computing PH on a flag complex (which has been done in several studies of PH in granular force networks [164–169]) only counts loops that include four or more particles. That is, it does not count 3-particle loops (which are sometimes called ‘triangular loops’). Loops with four or more particles are associated with defects, as they would not exist in a collection of monosized disks that are packed perfectly and as densely as possible into a ‘crystalline’ structure (which has only triangular loops) [165]. In Fig. 8c–f, we show the sequence of complexes that correspond to the filtration over the flag complex. One descends the four threshold levels (i.e., edge weights), beginning with the largest ($$\theta_{4}$$) and ending with the smallest ($$\theta_{1}$$). In Fig. 8g and h, we show the corresponding PDs for $$\beta_{0}$$ and $$\beta_{1}$$. It is helpful to discuss a few features of these diagrams. In $$\mathrm{PD}_{0}$$, we observe four points that are born at $$\theta_{4}$$; these points correspond to the four connected components that emerge at the first level of the filtration in Fig. 8c. Two of the components merge into one component at $$\theta_{3}$$ (see Fig. 8d); this corresponds to the point at $$(\theta_{4},\theta_{3})$$. A new component forms at $$\theta_{3}$$ and dies at $$\theta_{2}$$; this is represented by the point at $$(\theta_{3},\theta_{2})$$ (see Fig. 8d). Additionally, two components born at $$\theta_4$$ die at $$\theta_2$$, corresponding to the two points at $$(\theta_{4},\theta_{2})$$. One can continue this process until the end of the filtration, where there is just a single connected component (see Fig. 8f). This component is born at $$\theta_{4}$$; it persists for all thresholds, and we use [165]’s convention to give it a death coordinate of $$-1$$; this yields the persistence point at $$(\theta_{4},-1)$$. In $$\mathrm{PD}_{1}$$, we observe that a loop emerges at $$\theta_{3}$$ (see Fig. 8d), and it is then filled by triangles at $$\theta_{2}$$ (see Fig. 8e), leading to the point at $$(\theta_{3},\theta_{2})$$. Three more loops are born at $$\theta_{2}$$ and never die (see Fig. 8e); again using the convention in [165], we assign these features a death coordinate of $$0$$, so there are three persistence points at $$(\theta_{2},0)$$. Finally, one more loop appears at $$\theta_{1}$$ and does not die (see Fig. 8e); this is represented by a point at $$(\theta_{1},0)$$. Kramár et al. [165] gave an in-depth exposition of how to apply PH to granular networks, and we refer interested readers to this article for more information. Because studying PH is a general mathematical approach, it can be applied to different variations of force networks and can also be used on networks constructed from different types of experimental data (e.g., digital image data, particle-position data or particle-interaction data). [165] also discussed a set of measures that can be used to compare and contrast the homology of force networks both within a single system (e.g., at two different packing fractions) and across different systems (e.g., if one uses particles of different sizes or shapes), and they explored the robustness of PH computations to noise and numerical errors. In Section 3.2.5, we further discuss applications of methods from algebraic and computational topology to granular materials. 2.3 Some considerations when using network-based methods Because there are many methods that one can use to analyse granular networks and many quantities that one can compute to measure properties of these networks, it is useful to discuss some relationships, similarities and distinctions between them. Naturally, the meaning of any given network feature depends on how the network itself is defined, so we focus the present discussion on the most common representation of a granular system as a network. (See Section 3.3 for discussions of other representations.) In this representation (see Fig. 2), nodes correspond to particles and edges correspond to contacts between particles. Additionally, edge weights can represent quantities such as normal or tangential forces between particles. In this type of granular network, it is important to be aware of which network quantities explicitly take into account spatial information or physical constraints in a system, which consider only network topology, and which consider only network geometry (i.e., both topology and edge weights, but not other information). Granular materials have a physical nature and are embedded in real space, so such considerations are extremely important. For discussions of how these issues manifest in spatial networks more generally, see [37]. One way to explicitly include spatial or physical information into network analysis is to calculate quantities that are defined from some kind of distance (e.g., a Euclidean distance between nodes), whether directly or through a latent metric space, rather than a hop distance. For example, as discussed in Section 2.2.2, one can define the edge length between two adjacent nodes from the physical distance between them, which allows quantities such as mean shortest path length, efficiency and some centrality measures to directly incorporate spatial information. However, traditional network features such as degree and clustering coefficient depend only on network connectivity, although their values are influenced by spatial effects. In Section 2.2.7, we also saw that one can incorporate physical constraints from granular networks into community-detection methods by using a geographical null model, rather than the traditional NG null model, in modularity maximization. Different computations in network analysis can also probe different spatial, topological or geometrical scales. For example, measures such as degree, strength and clustering coefficients are local measures that quantify information about the immediate neighbourhood of a node. However, measures such as the mean shortest path length and global efficiency are global in nature, as they probe large-scale network organization. In between these extremes are mesoscale structures. A network-based framework can be very helpful for probing various types of intermediate-scale structures, ranging from very small ones (e.g., motifs, such as small cycles) to larger ones (e.g., communities) and tools such as PH were designed to reveal robust structural features across multiple scales. Crucially, although there are some clear qualitative similarities and differences between various network-analysis tools (and there are some known quantitative relationships between some of them [29]), it is a major open issue to achieve a precise understanding of the relationships between different network computations. Moreover, in spatially-embedded systems (as in any situation where there are additional constraints) one can also expect some ordinarily distinct quantities to become more closely related to each other [37]. Furthermore, the fact that a granular particle occupies a volume in space (volume exclusion) gives constraints beyond what arises from embeddedness in a low-dimensional space. 3. Granular materials as networks We now review network-based models and approaches for studying granular materials. Over the past decade, network analysis has provided a novel view of the structure and dynamics of granular systems, insightfully complementing and extending traditional perspectives. See [1–5] for reviews of non-network approaches. Perhaps the greatest advantages of using network representations and associated tools are their natural ability to (1) capture and quantify the complex and intrinsic heterogeneity that manifests in granular materials (e.g., in the form of force chains) and to (2) systematically and quantitatively investigate how the structure and organization of a granular system changes when subjected to external loads or perturbations (such as compression, shear, or tapping). In particular, network science and related subjects provide a set of tools that help quantify structure (and changes in structure) over a range of scales—including local, direct interactions between neighbouring particles; larger, mesoscale collections of particles that can interact and reconfigure via more complicated patterns; and system-wide measurements of material (re)organization. It is thought that local, intermediate, and system-wide scales are all important for regulating emergent, bulk properties of granular systems, but it can be difficult to obtain a holistic, multiscale understanding of these materials. For example, particle-level approaches may not take into account collective organization that occurs on slightly larger scales, and continuum models and approaches that rely on averaging techniques may be insensitive to interesting and important material inhomogeneities [42–44, 46]. Network representations also provide a flexible medium for modeling different types of granular materials (and other particulate matter). For example, network analysis is useful for both simulation and experimental data of granular materials, and methods from complex systems and network science can help improve understanding of both dense, quasistatically-deforming materials as well as granular flows. In any of these cases, one often seeks to understand how a system evolves over the course of an experiment or simulation. To study such dynamics, one can examine a network representation of a system as one varies a relevant physical quantity that parameterizes the system evolution. For example, for a granular system in which the packing fraction increases in small steps as the material is compressed quasistatically, one can extract a network representation of the system at each packing fraction during the compression process and then study how various features of that network representation change as the packing fraction increases. Even a particular type of granular system is amenable to multiple types of network representations, which can probe different aspects of the material and how it evolves under externally applied loads. For instance, one can build networks based only on knowledge of the locations of particles (which, in some cases, may be the only information available) or by considering the presence or absence of physical contacts between particles. If one knows additional information about the elements in a system or about their interactions, one can construct more complicated network representations of it. For example, it has long been known that granular materials exhibit highly heterogeneous patterns of force transmission, with a small subset of the particles carrying a majority of the load along force chains [196, 197]. Recall from Section 1 that, broadly speaking, a force chain (which is also sometimes called a force network) is a set of contacts that carry a load that is larger than the mean load [8, 16], and the mean orientation of a force chain often encodes the direction of the applied stress [17]. We illustrated an example of force-chain structure in Fig. 1, and we further discuss force-chain organization in Section 3.2. Because of the nature of the distribution of force values and the interesting way in which forces are spatially distributed in a material, it is often very useful to consider network representations of granular materials that take into account information about inter-particle forces (see Section 3.2) and to use network-based methods that allow one to quantitatively investigate how the structure of a force network changes when one includes only contacts that carry at least some threshold force (see Sections 3.2.4 and 3.2.5). In our ensuing discussion, we describe several network constructions that have been used to study granular materials, indicate how they have been investigated using many of the concepts and diagnostics introduced in Section 2.2, and review how these studies have improved scientific understanding of the underlying, complex physics of granular systems. 3.1 Contact networks A contact network is perhaps the simplest way to represent a granular system. Such networks (as well as the term ‘contact network’) were used to describe granular packings long before explicitly network science-based approaches were employed to study granular materials; see, for example, [198, 199]. The structure of a contact network encodes important information about a material’s mechanical properties. As its name suggests, a contact network embodies the physical connectivity and contact structure of the particles in a packing (see Fig. 2). In graph-theoretic terms, each particle in the packing is represented as a node, and an edge exists between any two particles that are in physical contact with one another. Note that it may not always be possible to experimentally determine which particles are in physical contact, and one may need to approximate contacts between particles using information about particle positions, radii, and inter-particle distances. (See Section 3.5 for details.) By definition (and however it is constructed), a contact network is unweighted and undirected, and it can thus be described with an unweighted and undirected adjacency matrix (see Section 2.1): $$A_{ij} = \left\{\begin{array}{@{}ll} 1, \text{ if particles}\,i\, \text{and}\,j\,\text{are in contact}, \\[3pt] 0, \text{ otherwise.} \end{array} \right.$$ (3.1) Because the organization of a contact network depends on and is constrained by the radii of the particles and their locations in Euclidean space, a contact network is a spatially-embedded graph [37]. In Section 3.2.2, we will see that this embedding in physical space has important consequences for the extraction of force-chain structures via community-detection techniques (see Section 2.2.7). In Fig. 9, we show an example of a contact network generated from a discrete element method (DEM) simulation (see Section 3.5) of biaxial compression [200]. The granular system in this figure is polydisperse, as it has more than two types of particles. (In this case, the particles have different sizes.) If all particles are identical in a granular system, it is called monodisperse; if there are two types of particles in a system, it is called bidisperse. In practice, although the presence or absence of a contact is definitive only in computer simulations, one can set reasonable thresholds and perform similar measurements in experiments [201] (see Section 3.5). It is also important to note that packing geometry and the resulting contact network do not completely define a granular system on their own. In particular, one can associate a given geometrical arrangement of particles with several configurations of inter-particle forces that satisfy force and torque balance constraints and the boundary conditions of a system [202–205]. This is a crucial concept to keep in mind when conducting investigations based only on contact networks, and it also motivates the inclusion of contact forces to construct more complete network representations of granular systems (see Section 3.2). Fig. 9. View largeDownload slide An example of a contact network from a DEM simulation of a densely-packed 2D system of polydisperse, disk-shaped particles. In this case, the granular material was subjected to biaxial compression and constrained to move along a plane. This snapshot corresponds to a network at an axial strain before full shear-band formation. We adapted this figure, with permission, from [58]. Fig. 9. View largeDownload slide An example of a contact network from a DEM simulation of a densely-packed 2D system of polydisperse, disk-shaped particles. In this case, the granular material was subjected to biaxial compression and constrained to move along a plane. This snapshot corresponds to a network at an axial strain before full shear-band formation. We adapted this figure, with permission, from [58]. In the remainder of this subsection, we review some of the network-based approaches for characterizing contact networks of granular materials and how these approaches have been used to help understand the physical behaviour of granular matter. We primarily label the following subsubsections according to the type of employed methodology. However, we also include some subsubsections about specific applications to certain systems. 3.1.1 Coordination number and node degree. One can study a contact network in several ways to investigate different features of a granular system. We begin our discussion by associating the mean node degree of a contact network with the familiar and well-studied coordination number (i.e., contact number) $$Z$$. Although early investigations of granular materials did not consciously make this connection, the mean degree and coordination number are synonymous quantities. The contact degree $$k_i$$ of particle $$i$$ is the number of particles with which $$i$$ is directly in contact, and one can calculate it easily from an adjacency matrix [see Eq. (2.5)]. A contact network is undirected, so its adjacency matrix $${\bf A}$$ is symmetric, and its row sum and column sum each yield a vector of node degrees. The mean degree $$\langle k \rangle$$ of a contact network [see Eq. (2.7)] is then the mean coordination number (i.e., contact number) $$Z$$, and it gives the mean number of contacts per particle. As we noted previously (see Section 2.2.1), $$Z$$ is an important quantity in granular systems because of its connection with mechanical stability and rigidity—which, loosely speaking, is the ability of a system to withstand deformations—and its characterization of the jamming transition [63, 206] and other mechanical properties. In particular, the condition for mechanical stability—that is, the condition to have all translational and rotational degrees of freedom constrained so that there is force and torque balance—in a packing of frictionless spheres in $$d$$ dimensions [61, 63, 64, 206, 207] is $$Z \geq 2d \equiv Z_{\mathrm{iso}}.$$ (3.2) The isostatic number $$Z_{\mathrm{iso}}$$ indicates the condition for isostaticity, which is defined as the minimum contact number that is needed for mechanical stability. One can use the coordination number (which is often tuned by changing the packing fraction $$\phi$$) as an order parameter to describe the jamming transition for frictionless spheres in 2D and 3D [63, 64, 207]. Specifically, there is a critical packing fraction $$\phi_{c}$$ such that below $$\phi_{c}$$, the contact number for these systems is $$Z = 0$$ (i.e. there are no load-bearing contacts), and at the critical packing fraction $$\phi_{c}$$, the contact number jumps to the critical value $$Z_c = Z_{{\mathrm iso}} = 2d$$. One can also generalize the use of the coordination number in order to examine mechanical stability and jamming in granular systems of frictional spheres. In these systems, the condition for stability is \begin{align} Z &\geq Z^{m}_{{\mathrm{iso}}}, \notag \\ Z^{m}_{{\mathrm{iso}}} &\equiv (d+1) + \frac{2N_{m}}{d}, \end{align} (3.3) where $$N_{m}$$ is the mean number of contacts per particle that have tangential forces $$f_{t}$$ equal to the so-called Coulomb threshold—that is, $$N_{m}$$ is the mean number of contacts per particle with $$f_{t} = \mu f_{n}$$, where $$\mu$$ is the coefficient of friction and $$f_{n}$$ is the normal force [64, 207, 208]—and $$Z^{m}_{{\mathrm{iso}}}$$ again designates the condition for isostaticity. Results from experimental systems have demonstrated that contact number also characterizes the jamming transition in frictional, photoelastic disks [201]. The coordination number has been studied for several years in the context of granular materials and jamming, and it is fruitful to connect it directly with ideas from network science. Several recent studies have formalized the notion of a contact network, and they deliberately modeled granular systems as such networks to take advantage of tools like those described in Section 2.2. Such investigations of contact networks allow one to go beyond the coordination number and further probe the rich behaviour and properties of granular materials—including stability and the jamming transition [63], force chains [8, 9, 16], and acoustic propagation [21, 23, 42, 45]. Perhaps the simplest expansion of investigations into the role of the coordination number is the study of the degree distribution $$P(k)$$ of the contact network of a packing. Calculating degree distributions can provide potential insights into possible generative mechanisms of a graph [29, 158], although one has to be very careful to avoid over-interpreting the results of such calculations [209]. In granular physics, it has been observed that the degree distribution of a contact network can track changes in network topology past the jamming transition in isotropically compressed simulations of a 2D granular system [73]. Specifically, the peak of $$P(k)$$ shifts from a lower value of $$k$$ to a higher value of $$k$$ near the transition. Moreover, changes in the mean degree $$\langle k \rangle$$ and its standard deviation can anticipate the onset of different stages of deformation in DEM simulations (i.e., molecular-dynamics simulations) of granular systems under various biaxial compression tests [58, 210]. 3.1.2 Investigating rigidity of a granular system using a contact network. An important area of research in granular materials revolves around attempts to (1) understand how different types of systems respond when perturbed, and (2) determine what features of a system improve its integrity under such perturbations. As we noted in Section 3.1.1, it is well-known that coordination number (and hence node degree) is a key quantity for determining mechanical stability and understanding jamming in granular materials. However, contact networks obviously have many other structural features, and examining them can be very helpful for providing a more complete picture of these systems. To the best of our knowledge, the stability of granular materials was first studied from a graph-theoretic standpoint in the context of structural rigidity theory [172, 179, 211], which has also been applied to amorphous solids more generally [62]. In structural rigidity theory, thought to have been studied originally by Maxwell [212], rods of fixed length are connected to one another by hinges, and one considers the conditions under which the associated structural graphs are able to resist deformations and support applied loads (see Fig. 10). A network is said to be minimally rigid (or isostatic) when it has exactly the number of bars needed for rigidity. This occurs when the number of constraints is equal to the number of degrees of freedom in the system (i.e., when Laman’s minimal-rigidity criterion is satisfied). The network is flexible if there are too few rods, and it is overconstrained (i.e., self-stressed) if there are more rods than needed for minimal rigidity. Triangles are the smallest isostatic structures in two dimensions [213–215]; there are no allowed motions of these structures that preserve the lengths and connectivity of the bars, so triangles (i.e., 3-cycles) do not continuously deform due to an applied force. In comparison, 4-cycles are structurally flexible and can continuously deform from one configuration to another while preserving the lengths and connectivity of the rods (see Fig. 10). Fig. 10. View largeDownload slide Some ideas from structural rigidity theory. (a) An example of a 3-cycle and a 4-cycle that one can examine using concepts from structural rigidity theory by considering the edges to be rods of fixed length that are connected to one another by rotating hinges. (b) Triangular structures are rigid under a variety of applied forces (represented as red arrows), whereas (c) squares can deform under such perturbations. By Laman’s theorem, a 2D network with $$N$$ nodes is minimally rigid if it has exactly $$2N-3$$ edges, and each of its subgraphs satisfies the analogous constraint (so an $$\tilde{N}$$-node subgraph has no more than $$2\tilde{N}-3$$ edges). The network in panel (b) satisfies this criterion, but the network in panel (c) does not. Fig. 10. View largeDownload slide Some ideas from structural rigidity theory. (a) An example of a 3-cycle and a 4-cycle that one can examine using concepts from structural rigidity theory by considering the edges to be rods of fixed length that are connected to one another by rotating hinges. (b) Triangular structures are rigid under a variety of applied forces (represented as red arrows), whereas (c) squares can deform under such perturbations. By Laman’s theorem, a 2D network with $$N$$ nodes is minimally rigid if it has exactly $$2N-3$$ edges, and each of its subgraphs satisfies the analogous constraint (so an $$\tilde{N}$$-node subgraph has no more than $$2\tilde{N}-3$$ edges). The network in panel (b) satisfies this criterion, but the network in panel (c) does not. Extending a traditional network of rods and hinges, concepts from structural rigidity theory yield interesting insights into contact networks of particulate matter. See [216] for a discussion of some of the earliest applications of such ideas to disordered systems and granular materials. Moukarzel [217] used structural rigidity theory to derive conditions for the isostaticity of a granular packing, and he tied the fact that random packings of particles tend to be isostatic to the origin of instabilities in granular piles. Later, similar concepts were used to show that in granular networks, cycles with an even number of edges allow contacting grains to roll without slipping when subject to shear; however, these relative rotations are ‘frustrated’ in cycles with an odd number of edges, so such cycles can act as stabilizing structures in a network [85]. Several later studies (such as [58, 73, 79–81, 83, 84, 210, 218–220]) have confirmed that contact loops are often stabilizing mesoscale features in a contact network of a granular material. We specifically consider the role of cycles in granular contact networks in Section 3.1.3. Another type of network approach for understanding rigidity in granular systems is rigidity percolation [178, 179, 211] (see Section 2.2.9). Feng [170] conducted an early investigation of an idealized version of bond percolation in a granular context. It is now known that hallmarks of this bond-percolation transition occur below isostaticity: Shen et al. [175] identified that a percolating (i.e., system-spanning) cluster of non-load-bearing contacts forms at a packing density below the jamming point. One can also use a rigidity-percolation approach to examine if a network is both percolating and rigid (see Section 2.2.9). Note that a rigid granular network is also percolating, but a percolating network need not be rigid. Rigidity percolation relies on tabulating local constraints via a pebble game [172], which reveals connected, rigid regions (sometimes called ‘rigid clusters’) in a network. In a series of articles [221–224] on simulated packings, Schwarz and co-workers went beyond Laman’s minimal-rigidity criterion to investigate local versus global rigidity in a network, the size distribution of rigid clusters, the important role of spatial correlations, and the necessity of force balance. Building on the above work, [177] recently utilized a rigidity-percolation approach to identify floppy versus rigid regions in slowly-sheared granular materials and to characterize the nature of the phase transition from an underconstrained system to a rigid network. See also the recent article [225]. 3.1.3 Exploring the role of cycles. We now consider the role of circuits (i.e., the conventional network notion of cycles, which we discussed in Section 2.2.3) in granular contact networks. Cycles in a contact network can play crucial stabilizing roles in several situations. Specifically, as we will discuss in detail in this section, simulations (and some experiments) suggest that (1) odd cycles (especially 3-cycles) can provide stability to granular materials by frustrating rotation among grains and by providing lateral support to surrounding particles, and that (2) a contact network loses these stabilizing structures as the corresponding granular system approaches failure. Noting that 3-cycles are the smallest arrangement of particles that can support (via force balance) a variety of 2D perturbations to a compressive load without deforming the contact structure, Ottino [81] studied the effects of friction and tilting on the evolution of contact-loop organization in a granular bed. In their simulations, they implemented tilting by incrementally increasing the angle of a gravity vector with respect to the vertical direction, while preserving the orientation of the granular bed and maintaining quasistatic conditions. In untilted granular packings, they observed that lowering inter-particle friction yields networks with a higher density of 3-cycles and 4-cycles, where they defined the ‘density’ of an $$l$$-cycle to be the number of $$l$$-cycles divided by the total number of particles. By examining the contact network as they varied tilting angle, Ottino [81] also observed that the density of 4-cycles increases prior to failure—likely due to the fracture of stabilizing 3-cycles—and that this trend was distinguishable from changes in the coordination number alone. Cycles have also been studied in the context of DEM simulations of dense, 2D granular assemblies subject to quasistatic, biaxial compression [58, 84, 210, 218]. In many of these studies, the setup consists of a collection of disks in 2D that are compressed slowly at a constant strain rate in the vertical direction, while allowed to expand under constant confining pressure in the horizontal direction [58, 84, 218]. In another variation of boundary-driven biaxial compression, a sample can be compressed with a constant volume but a varying confining pressure [210]. Before describing specifics of the network analysis for these systems, it is important to note that for the previously described conditions, the axial strain increases in small increments (i.e., ‘steps’) as compression proceeds, and one can extract the inter-particle contacts and forces at each strain value during loading to examine the evolution of a system with respect to strain. Additionally, these systems undergo a change in behaviour from a solid-like state to a liquid-like state, and they are characterized by different regimes of deformation as one increases axial strain [200]. In particular, the granular material first undergoes a period of strain hardening, followed by strain softening, after which it enters a critical state. In the strain-hardening regime, the system is stable and the shear stress increases monotonically with axial strain up to a peak value. After the peak shear stress, strain softening sets in; this state is marked by a series of steep drops in the shear stress that indicate reduced load-carrying capacity. Finally, in the critical state, a persistent shear band has fully formed, and the shear stress fluctuates around a steady-state value. The shear band is a region of localized deformation and provides one signature of material failure [226]. Inside the shear band, force chains can both form and buckle [227]. One can also associate increases in the energy dissipation rate of the system with particle rearrangements (such as those that occur during force-chain buckling) and loss of stability in the material. Examining the temporal evolution of cycles in such an evolving granular contact network can reveal important information about changes that occur in a material during deformation. Using DEM simulations (see the previous paragraph), Tordesillas et al. [84] computed the total number of cycles of different lengths in a minimal cycle basis (see Section 2.2.3) of a contact network at each strain state during loading, and they observed that there are many more 3-cycles and 4-cycles than longer cycles in the initial, solid-like state of the system. However, as axial strain increases and one approaches the maximum shear stress, the total number of 3-cycles falls off steeply. (The same is true for 4-cycles, though the effect is less dramatic.) Additionally, during axial-strain steps (i.e., axial-strain ‘intervals’) corresponding to drops in shear stress, Tordesillas et al. [84] observed large increases in the number of 3-cycles and 4-cycles that open up to become longer cycles. In Fig. 11a, we show an example of the evolution of cycle organization with increasing axial strain for a subset of particles from a DEM simulation of a granular material under biaxial compression carried out by Walker and Tordesillas [58]. The authors observed that in this system, both the global clustering coefficient $$C$$ [see Eq. (2.13)] and the mean subgraph centrality $$Y$$ [see Eq. (2.19)] decrease with increasing axial strain, drop sharply at peak shear stress, and then level out (see Fig. 11b and c). Recalling that $$C$$ is a measure of triangle density in a graph and that subgraph centrality measures participation of nodes in closed walks (with more weight given to shorter paths), these results also imply that the loss of small cycles co-occurs with the deformation and failure of a system due to increasing load. Walker and Tordesillas [58] also computed the network bipartivity $$R$$ [95] of the contact network to quantify the contribution to mean subgraph centrality $$Y$$ from closed walks of even length [see Eq. 2.20]. They observed that $$R$$ increases with increasing axial strain, revealing that closed walks of even length become more prevalent during loading (see Fig. 11c). The authors suggested that this trend may be due to a decrease in the prevalence of 3-cycles (which are stabilizing, as discussed in Section 2.2.3 and elsewhere). Tordesillas et al. [218] also examined the stability of cycles of various lengths in both DEM simulations and experimental data, and they observed that, during loading, 3-cycles tend to be more stable (as quantified by a measure of stability based on a structural-mechanics framework [228]) than cycles of other lengths in a minimal cycle basis of the network. Fig. 11. View largeDownload slide Evolution of cycles in a deforming granular material. (a) One can track a subset of particles and their corresponding contact network from a DEM simulation for increasing axial strain values $$|\varepsilon_{22}|$$. The system consists of 5098 spherical, polydisperse particles that were subjected to a quasistatic, biaxial compression test. At the smallest displayed axial strain, the set of particles in this figure yields a network that is composed of 3-cycles and 4-cycles. During loading, contacts are lost and longer cycles arise until only a single 9-cycle remains. (b) One way to quantify these structural changes is by calculating the global clustering coefficient $$C$$ (solid curve), which undergoes a sharp drop at peak stress, signifying the onset of material failure. (The dashed curve shows the standard deviation of the distribution of local clustering coefficients $$C_{i}$$). (c) A decrease in mean subgraph centrality $$Y$$ (solid curve) also illustrates the loss of short cycles during deformation. More specifically, the mean network bipartivity $$R$$ (dashed curve) increases with axial strain, highlighting that, during loading, closed walks of even (respectively, odd) length contribute more (respectively, less) to the mean subgraph centrality of the contact network. We adapted this figure, with permission, from [58]. Fig. 11. View largeDownload slide Evolution of cycles in a deforming granular material. (a) One can track a subset of particles and their corresponding contact network from a DEM simulation for increasing axial strain values $$|\varepsilon_{22}|$$. The system consists of 5098 spherical, polydisperse particles that were subjected to a quasistatic, biaxial compression test. At the smallest displayed axial strain, the set of particles in this figure yields a network that is composed of 3-cycles and 4-cycles. During loading, contacts are lost and longer cycles arise until only a single 9-cycle remains. (b) One way to quantify these structural changes is by calculating the global clustering coefficient $$C$$ (solid curve), which undergoes a sharp drop at peak stress, signifying the onset of material failure. (The dashed curve shows the standard deviation of the distribution of local clustering coefficients $$C_{i}$$). (c) A decrease in mean subgraph centrality $$Y$$ (solid curve) also illustrates the loss of short cycles during deformation. More specifically, the mean network bipartivity $$R$$ (dashed curve) increases with axial strain, highlighting that, during loading, closed walks of even (respectively, odd) length contribute more (respectively, less) to the mean subgraph centrality of the contact network. We adapted this figure, with permission, from [58]. Minimal cycle bases and the easier-to-compute subgraph centrality have also been used to examine fluctuations in kinetic energy in simulations of deforming sand. Walker et al. [79] computed a minimal cycle basis and then constructed cycle-participation vectors (see Section 2.2.3) from a contact network after each strain step (i.e., at each strain state) during loading. They observed that temporal changes in the cycle-participation vectors of the particles between consecutive strain steps are correlated positively with temporal changes in kinetic energy over those steps. They also observed that large values in the temporal changes of particle cycle-participation vectors and particle subgraph centrality occur in the shear-band region. Walker et al. [80] also studied a minimal cycle basis and corresponding cycle-participation vectors to examine structural transitions in a 3D experimental granular system of hydrogel spheres under uniaxial compression. As pointed out in [79], developing quantitative predictors that are based on topological information alone is extremely important for understanding how failure and rearrangements occur in systems in which energy or force measurements are not possible. Examining cycles in contact networks can also shed light on the behaviour of force chains. The stability, load-bearing capacity, and buckling of force chains depend on neighbouring particles (so-called spectator grains) to provide lateral support [19, 20, 229]. Because 3-cycles appear to be stabilizing features, it is interesting to consider the co-evolution of force chains and 3-cycles in a contact network. Such an investigation requires a precise definition of what constitutes a force chain, so that it is possible to (1) extract these structures from a given packing of particles and (2) characterize and quantify force-chain properties. Several definitions of force chains have been proposed; see, for example, [14, 16, 27, 230]. The studies that we describe in the next three paragraphs used a notion of ‘force chains’ from [14, 231], in which force chain particles are identified based on their particle-load vectors (where each particle is associated with a single particle-load vector that indicates the main direction of force transmission). More specifically, a single chain is a set of three or more particles for which the magnitude of each of their particle-load vectors is larger than the mean particle-load vector magnitude over all particles, and for which the directions of the particle load vectors are, within some tolerance, aligned with one another (i.e., they are ‘quasilinear’). We note that it is important in future work to conduct network-based studies of force-chain structure for different definitions of force chains and investigate if there are differences in their associated network properties. Using DEM simulations of a densely packed system of polydisperse disks under biaxial loading—i.e., compressed quasistatically at a constant strain rate in the vertical direction, while allowed to expand under constant confining pressure in the horizontal direction—Tordesillas et al. [84] quantified the co-evolution of force chains and 3-cycles in several ways. For example, they computed a minimal cycle basis (see Section 2.2.3) of a contact network and then examined (1) the ratio of 3-cycles to the total number of cycles in which particles from a force-chain participate and (2) the force chain’s 3-cycle concentration, which they defined as the ratio of 3-cycles involving force-chain particles to the total number of particles in the force chain. When averaged over all force chains, the above two measures decrease rapidly with increased loading. Additionally, Tordesillas et al. [84] observed that force chains that do not fail by buckling (see [200] for how ‘buckling’ was defined) have a larger ratio of 3-cycle participation to total cycle participation than force chains that do buckle. Tordesillas et al. [218] observed, in both DEM simulations of biaxial loading (see above) and 2D photoelastic disk experiments under pure shear, that a particular measure (developed by Bagi [228]) of structural stability of force chains is correlated positively with the mean of the local clustering coefficient [see Eq. (2.12)] over force-chain particles. They also pointed out that 3-cycles are more stable structures than cycles of longer length during loading and that force chains with larger 3-cycle participation tend to be more structurally stable. These observations suggest that cycles—and especially 3-cycles—in contact networks are stabilizing structures that can provide lateral support to force chains. It would be interesting to study these ideas further, and to relate them to structural rigidity theory (see Fig. 10 and Section 3.1.2), especially in light of the difference between 3-cycles (which are rigid) and deformable cycles (e.g., 4-cycles). Cycles have also been examined in granular contact networks that were extracted from DEM simulations of 3D, ellipsoidal-particle systems subject to triaxial compression [219]. Similar to the aforementioned results from simulations of 2D systems with disk-shaped particles, the number of 3-cycles in a minimal cycle basis of the contact networks (and the global clustering coefficient [see Eq. (2.13)]) initially decrease and then saturate with increasing load, and particles in force chains have a larger number of 3-cycles per particle than particles that are not in force chains. Tordesillas et al. [219] also observed that the set of 3-cycles that survive throughout loading tend to lie outside the strain-localization region (where force-chains buckle). The dearth of 3-cycles in certain regions in a material may thus be a signature of strain-localization zones. Another article to note is [210], which examined and compared the temporal evolution of cycles (and several other contact-network diagnostics) in a set of DEM simulations using a variety of different material properties and boundary conditions. In another interesting study, Walker et al. [220] examined the phenomenon of aging [232, 233]—a process in which the shear strength and stiffness of a granular material both increase with time—in collections of photoelastic disks subject to multiple cycles of pure shear under constant volume. Because aging is a slow process, it can be difficult both to uncover meaningful temporal changes in dynamics and to characterize important features in packing structure that accompany aging. To overcome these challenges, Walker et al. [220] first analysed the time series of the stress ratio (using techniques from dynamical-systems theory) to uncover distinct temporal changes in the dynamics of the system. (See [220] for details.) After each small, quasistatic strain step, they also extracted the contact network of the packing at that time to relate aging to changes in topological features of the network structure. As the shear-jammed regime is approached during prolonged cyclic shear, they observed, on average, that force chains are associated with more 3-cycles and 4-cycles from the minimal cycle basis. We have just discussed many articles that concern transitions in granular matter from a solid-like regime to a liquid-like regime. One can also use changes in the loop structure of a contact network to describe the opposite transition, in which a granular material transitions from an underconstrained, flowing state to a solid-like state during the process known as jamming (see Section 3.1.1). Studying 2D frictional simulations of isotropically compressed granular packings, Arévalo et al. [73] examined a granular contact network as they varied packing fraction. They observ