Evolution of network architecture in a granular material under compression

As a granular material is compressed, the particles and forces within the system arrange to form complex heterogeneous structures. Force chains are a prime example and are thought to constrain bulk properties such as mechanical stability and acoustic transmission. However, characterizing the dynamic nature of mesoscale architectures in granular systems can be challenging. A growing body of work has shown that graph theoretic approaches may provide a useful foundation for tackling these problems. Here, we extend current approaches by utilizing multilayer networks as a framework for directly quantifying the evolution of mesoscale architecture in a compressed granular system. We examine a quasi-two-dimensional aggregate of photoelastic disks, subject to biaxial compression through a series of small, quasistatic steps. Treating particles as network nodes and inter-particle forces as network edges, we construct a multilayer network by linking together the series of static force networks that exist at each strain step. We then extract the inherent mesoscale structure from the system by using a generalization of community detection methods, and we define quantitative measures to characterize the reconfiguration and evolution of this structure throughout compression. By separately considering the network of normal and tangential forces, we find that they display different structural evolution. To test the sensitivity of the network model to particle properties, we examine whether the method can distinguish a subsystem of low-friction particles within a bath of higher-friction particles. We find that this can be done by considering the network of tangential forces. The results discussed throughout this study suggest that these novel network science techniques may provide a direct way to compare and classify data from systems under different external conditions or with different physical makeup.


I. INTRODUCTION
Granular materials [1] come in many forms, from soils, sands, and grains to powders and pharmaceuticals.However, despite their prevalence, there are still open questions about how the seemingly simple interactions of contact forces lead to the observed emergent behavior of these systems.An active area of research lies in understanding the mechanisms that govern deformation in granular materials subjected to compression and shear.Under both of these external perturbations, the force network of granular systems exhibits complex and inhomogeneous structure in the form of strongly interacting collections of particles known as force chains (see Fig. 1C ) and natural paradigm in which to study granular media.Many of these analyses have focused on the characterization of discrete sets of static granular force networks throughout compression [21][22][23][24][25][26], tapping [27], or tilting [28], using traditional graph metrics such as degree, clustering coefficients, and cycles of different lengths.Other work has probed the dynamical nature of sheared systems by considering time-evolving networks of brokenlinks [29,30], and grain property networks have been used to understand rearrangements in discrete element simulations of compressed systems [31].Methods from algebraic topology, and in particular, persistent homology [32,33] have also been used to quantify the evolution of force networks, providing important insights into the nature of compressed [34][35][36] and tapped [37][38][39] granular materials.
One reason for the promise of network-based measures is that they can assess material architecture at a range of length scales, including the important mesoscale regime.Recent work by Bassett et al. [40] showed that a network clustering technique known as community detection could be used to extract the underlying mesoscale force chain structure from static granular networks.In this study, we extend the static model, and suggest that multilayer networks may be a particularly promising framework in which to simultaneously examine the mesoscale architecture of granular systems, as well as to probe network evolution and reconfiguration in a straightfoward manner.This novel approach has thus far been unexplored.
Multilayer networks encompass several different types of complex graph constructions, and the word can take on a number of meanings depending on the context (see [41,42] for comprehensive reviews).For example, a multilayer network may capture different types of connections between nodes, may quantify interactions between different systems, or may be used to study dynamical processes that occur across time.Here, we focus on a specific subset of these possibilities.In particular, we restrict ourselves to temporal networks with diagonal and ordinal inter-layer couplings.A temporal network consists of a sequential series of static graphs (the layers) ordered such that time-dependence is accounted for.Diagonal couplings mean that a node in one layer is only connected to itself in other layers.Finally, ordinal inter-layer couplings only allow connections between layers that are adjacent to each other in time.In this study, we are interested in the evolution of a granular material as it undergoes biaxial compression.We thus represent a snapshot of the system at a particular point in its evolution as a spatially embedded graph where particles are nodes and inter-particle forces are weighted edges.Repeating this process at several discrete strain steps yields an ordered set of static networks, which can then be combined into a single multilayer graph with the temporal dependence set by the order of the strain steps.
We develop and apply this multilayer network formalism to experimental granular data [43] to (i) extract evolving mesoscale structure from the force network, (ii) understand how this architecture reconfigures throughout compressive (strain) steps, (iii) uncover physical properties of evolving mesostructures and relate them to measures of network rearrangement, and (iv) examine the impact of inter-particle friction on multilayer mesostructures.To achieve this, we represent an ensemble of granular packings as multilayer graphs using the set of force networks (both normal and tangential) obtained from each step above the jamming point.We then use multilayer community detection to extract groups of particles that evolve together throughout the compression procedure.This method allows us to directly characterize the progression of inherent mesoscale organization as a function of strain step.
The outline of this paper is as follows.In Sec.II we describe the granular experiments.Sec.III is dedicated to an explanation of the multilayer network model and community detection, which lay the theoretical foundations for the rest of the paper.A series of results describing the granular network evolution are presented in Sec.IV, and in Sec.V we discuss broader implications of our method and findings, and directions for future work.

II. EXPERIMENTAL METHODS
We study the biaxial compression of a granular monolayer on a nearly frictionless surface provided by an air table (Fig. 1).The system is composed of an inner subsystem and an outer bath which differ only in the interparticle friction coefficient.At the beginning of each compression cycle, the system is in a dilute, unjammed state and is then compressed quasistatically until it is highly stressed.By repeating this protocol many times, we generate an ensemble of configurations for which we record particle positions (Fig. 1B ), and use photoelastic measurements (Fig. 1C ) to calculate contact forces [44].The presence of the subsystem and bath allows us to characterize both the system as a whole, and to investigate the physical differences that exist between the high-and low-friction regions.
The inner subsystem and bath are composed of 100 and 904 particles, respectively.The two systems are both composed of a bidisperse mixture of disks (diameters d s = 11 mm and d l = 15.4 mm, in equal concentrations), which differ only in the inter-particle friction coefficient.The friction coefficient for the inner subsystem is µ < 0.1 and for the high friction bath is µ ≈ 0.8.
A single CCD camera located above the apparatus records three separate images of the system: unpolarized white light for the particle positions, polarized light for recording the photoelastic response, and ultraviolet light for identifying the subsystem particles, as shown in Fig. 1.Using the photoelastic image, we measure the normal and tangential forces at each contact on each disk in the assembly.To calculate the forces, we use a nonlinear least squares optimization algorithm [45,46] to minimize the error between the observed and fitted image of the Experimental setup.(A) A schematic of the apparatus showing two walls bi-axially compressing an array of disk-shaped particles composed of an outer subsystem (black, high µ) and an inner subsystem (red, low µ).(B) The first image is taken with unpolarized white light and is used to locate particle positions.(C) The network of the photoelastic stress pattern, visualized with a second image taken with polarized light.This image allows for the calculation of the normal and tangential forces at each contact.A third image (inset) taken with ultraviolet light is used to identify the subsystem particles, which are tagged with fluorescent ink.
particle.Details and source code are available for download at [44].The third image, taken using a black-light illumination, identifies which particles are low-friction via their fluorescent marking.
Each compression cycle begins with an initial dilute state with global packing fraction Φ ∼ 0.6.In the dilute state, the granular system is unjammed and unable to support stress.Two walls then biaxially compress the system in a series of small steps of constant size (∆Φ = 0.009, equivalently ∆x = 0.3 mm or 0.02r l ).As the system is compressed via these discrete, quasistatic steps, we observe the percolation of force chains throughout both the bath and the subsystem at a value Φ J .This is the onset of rigidity, and as the system is further compressed beyond this point, the contact forces grow in strength and the average number of contacts per particle increases.For each of the 97 configurations we locate Φ J as the step at which photoelastic signals are first present, and consider the evolution of the system due to the subsequent applied strain steps.Further details of the experimental apparatus and measurements were published previously [43].

A. Particle tracking
It is important to note that the construction of a multilayer network requires knowledge about which node (particle) is which from one layer (strain step) to the next [41,42].Under the dynamics described above, we are assured that no particles are removed from or added to the system, and we similarly require that the multilayer network also has this constraint.(In general, multilayer graphs can indeed be constructed for growing systems, where the number of nodes is constantly changing.)In order to correctly identify the particles in each layer, we use the Blair-Dufresne particle tracking algorithm [47].This algorithm, implemented in MATLAB, requires the choice of a "displacement" parameter, which is an estimate of the maximum distance that a particle moves between consecutive frames.As we do not expect large particle displacements in the jammed packings, we initialize the tracking with a displacement parameter value that is much less than the minimum particle diameter d s , and increase the allowed displacement in small steps until all particles are tracked consistently across all steps.

A. Multilayer network representation
In general, temporal multilayer networks are mathematical objects that describe and quantify the evolution of networked systems as a function of time or some other variable of interest [41,42].The construction of such a multilayer network first requires a set of individual networks that describe the system at discrete steps in its evolution.These static "layers" can each be represented as an adjacency matrix, A, describing the connectivity and weight connecting the nodes in the given layer.The set of adjacency matrices may then be combined into a rank-3 adjacency tensor, A [48], to form a multilayer network representation of the system.The elements of A are defined such that where a ijl is the weight between nodes i and j in layer l, and L is the total number of layers (steps).
For the present situation, we are interested in how the force network of a granular packing reconfigures through a sequence of compressive steps above Φ J .We thus wish to represent the system as an ordered, temporal multilayer network, which captures the evolution of the system as a function of compression.To form such a network, we thus let nodes be particles, weighted edges be the forces between contacting particles, and the number of strain steps be the third dimension across which we observe changes in network structure.Using the force information obtained from the photoelastic disk experiments, we construct two multilayer force networks for a given experimental run, one using the normal forces, F n , between particles and another using the magnitude of the tangential forces |F t |.Specifically, we denote the normal force adjacency tensor as A n ijl , defined as in Eq. 1 with a ijl = f n ijl , where f n ijl is the normal force between particles i and j at step l.Similarly, we write the tangential force adjacency tensor as A t ijl , defined as in Eq. 1 with a ijl = |f t ijl |, where f t ijl is the tangential force between particles i and j at step l.This process is repeated for all of the experimental configurations, resulting in a large ensemble of multilayer networks.In Fig. 2, we show an example set of photoelastic images with the corresponding network of normal forces overlaid.
In addition to the intra-layer weights, multilayer networks have a second set of inter -layer edges that connect nodes in different layers of the network.In this study, we consider only diagonal and ordinal inter-layer couplings, meaning that a given particle is connected only to itself (i.e.diagonal coupling) inter-layer edges exist only between adjacent layers (i.e ordinal coupling).The inter-layer weights are crucial to the structure of the network; in Sec.III B 4, we discuss how these couplings are chosen for the granular system at hand.This graphical construction is a powerful approach.Importantly, each layer of the adjacency tensor encodes both the topology (connectivity) as well as the strength of interactions between particles in the system at a given packing fraction.We will see that the extension to a multilayer framework not only allows one to study the static behavior of granular packings, but also promotes a direct investigation and characterization of the evolution of the system throughout the compressive steps.

B. Community detection
A compelling reason to use graphical representations of spatially embedded systems [49] is that network theoretic approaches provide a means to extract and characterize organization that is present at a range of length scales.Granular materials are a prime example of a complex system in which multiple length scales are relevant to a full understanding of the system and to a prediction of bulk properties.For example, the fundamental interactions in these systems are local in the sense that they occur between nearest neighbor particles.But under compression, those local interactions lead to important heterogeneous structure at the mesoscale in the form of force chains [2][3][4][5][6][7][8]50].
In the present study, we aim to both extract mesoscale structures from the force network, and then characterize the reconfiguration that occurs in these structures due to the applied strain steps.In order to accomplish this, we build upon recent approaches that utilize community detection techniques to extract groups of strongly interacting particles from the granular network [26,40], generalizing the methods to the multilayer regime.
A community or module in a network is a set of nodes that are densely interconnected amongst themselves, and relatively weakly connected to other nodes [51].The extraction of community structure is of general interest in network science, as it is thought that these mesoscale units are important to the function of many real systems [51].Several community detection methods exist [51,52]; here, we apply the popular method of modularity maximization, whereby nodes are partitioned into communities via maximization of a quality function known as modularity [53,54].

Single layer modularity maximization
In a single layer network with adjacency matrix A, modularity is given by where c i is the community of node i, c j is the community of node j, P ij is the expected edge weight between nodes i and j under a specified null model, and γ is the structural resolution parameter.In the overall normalization, m = 1 2 ij A ij , which is the weighted degree or strength Φ = 0.7842 Φ = 0.7854 Φ = 0.7831 Φ = 0.7836 Φ = 0.7860 Φ = 0.7866 Φ = 0.7872 Φ = 0.7884 FIG.
3. An example of single layer community structure.When single layer community detection is performed on the series of normal force networks above ΦJ , the groups of particles at each step are highly reminiscent of force chains.However, there is no notion of linking these mesostructures from one compressive step to the next.
of the network.The structural resolution parameter γ allows for the control of size and number of communities: smaller γ leads to larger communities, and larger γ leads to smaller communities.Maximization of Q sing with respect to the assignment of nodes to communities yields a partition in which intra-community connections are as strong as possible relative to the null model.It is important to note, however, that modularity maximization is NP-hard [55] and should therefore be repeated several times for the same network and set of parameters in order to obtain an ensemble of optimizations [56].Throughout this work, we use a Louvain-like locally greedy algorithm for modularity maximization [57,58].

A physically-informed null model
A proper choice of null model is vital to community detection techniques, as it affects both the interpretation and utility of the community structures obtained [56,59].The most commonly used null model in the literature is the Newman-Girvan (NG) model [53,54], which in the case of a static network is given by P ij = kikj 2m , where k i = j A ij is the strength of node i, and m is the total strength of the network, as before.In the Newman-Girvan model (which is also sometimes referred to as the configuration model ), P ij gives the expected edge weight between nodes i and j in a randomized network that has the same weighted degree distribution as the real network.As a randomized version of the real graph, the NG model is most appropriate to use as a null model in situations where all connections between nodes are at least possible.In many physical or spatially embedded systems, however, this is not the case, and there are constraints that prevent the existence of several edges.For example, in the granular networks considered here, edges can only exist between nearest neighbor particles, and it is imperative to consider this fact when designing a null model for these types of systems.A better choice in this instance is the physically informed geographical null model [40,56], defined to be where f is the average inter-particle force (either normal or tangential) in the network and B ij is the contact A schematic of a multilayer network with layer-dependent community structure.Each layer represents a static granular force network in which nodes (particles) are connected to one another via intra-layer weighted edges.These weights can be either the normal contact forces or the absolute value of the tangential contact forces.Additionally, the same particle in consecutive layers is linked to itself with an inter-layer coupling, ω.For clarity, we only show two such couplings, but these inter-layer edges exist between all particles (and across all layers).Evolving communities are extracted from the multilayer network to study the mesoscale organization in the system, and to understand how it reconfigures due to the compressive cycle.The community structure can be determined from the network by maximizing the multilayer modularity, Q multi (Eq.5).In this schematic, the particles belonging to different communities are labeled by different colors.Note that the same community can persist across all layers and reconfigure in terms of particle content and strength throughout compression.

matrix, with elements
Because this null model maintains the contact structure of the real network, it importantly takes into account the physical constraints on the possible patterns of connectivity between particles.Additionally, it selects for strongly connected sets of particles carrying forces larger than γ f .

Multilayer modularity maximization
When community detection is performed on a series of single layer granular networks, such as those obtained here after each applied strain step, the result is a set of independent partitions of particles into communities at each step (Fig. 3).As demonstrated in [40], the communities at a given step correspond to sets of strongly interacting groups of particles, and are clearly reminiscent of physical force chains observed in the photoelastic disk experiments.However, in this scheme, the community structure is not in any way linked from one layer to the next, barring any notion of continuation.Instead, the communities are treated as independent from one another, which is an inaccurate representation of the physics and further challenges our ability to directly capture the evolution of network structure and reconfiguration at the mesoscale.
To form a more dynamic picture of structural network evolution, we investigate the community structure of multilayer granular force networks by applying the recent generalization of modularity maximization to temporal networks [60,61].In this formulation, the multilayer modularity is defined to be (5) where A ijl is the (i, j) component of the l th layer of the adjacency tensor, P ijl is the (i, j) component of the l th layer of the null model tensor, and γ l is the structural resolution parameter for layer l.In addition to γ, the multilayer modularity requires another free parameter, ω (often referred to as an inter-layer coupling or temporal resolution parameter ) which sets the the strength of connections between layers.Namely, ω jlm is the strength of the coupling that links node j in layer l to itself in an adjacent layer m (i.e. the diagonal and ordinal coupling).The quantities c il and c jm are the community assignments of node i in layer l and node j in layer m, respectively.Defining the intraslice strength of node j in layer l as k jl = i A ijl and the strength of node j across layers as w jl = m ω jlm , then the multilayer strength of node j in layer l is given by κ jl = k jl + w jl .Finally, in the overall normalization, µ is the total strength of the adjacency tensor A, given by µ = 1 2 jl κ jl .In Fig. 4, we show a schematic of a multilayer granular force network with evolving community structure.Importantly, the communities can persist across all layers and we can track their reconfiguration in terms of particle content and strength throughout the series of strain steps.
As in the formulation of Q sing , the choice of null model in Q multi is an important one, particularly when considering systems with strict constraints on the allowed connectivity between nodes.For granular networks, we generalize the geographical null model to the multilayer regime.For community detection on the normal force network we use where f n l is the average of the normal component of the inter-particle forces at compressive step l and B ijl is the contact matrix at compressive step l.For the tangential network, we take the null model to be where |f t | l is the average of the absolute value of the tangential component of the inter-particle forces at compressive step l, and B ijl is the contact matrix at that step.
In each layer, we normalize the force network (A n ijl or (|A t ijl |) by the mean inter-particle force in the corresponding layer.Thus, after normalization we have f n l = 1 in Eq. 6 and |f t | l = 1 in Eq. 7, for all l.This normalization ensures that the community structure is not purely driven by the final layer, which will have the largest total edge weight due to it being the most compressed.
The two parameters in the multilayer modularity quality function, ω and γ, must be chosen by the investigator.Recall that the structural resolution parameter, γ, regulates the size and number of detected communities; in this study, we use γ l ≡ γ = 1 for all l.As can be seen from Eq. 5 and Eqs.6 and 7, the physical meaning of this value is that it selects for communities within each layer that have stronger than average force.The choice of ω is an active area of investigation; currently, there is no consensus in the literature on a single broadly-applicable method to determine the inter-layer coupling.In the present work, we make a physically informed choice.As described in the previous section, ω is in general a tensor that can take on different values between each layer or for different nodes.Since to our knowledge this is the first study on multilayer granular networks, we begin by taking the simplest case of a scalar inter-layer coupling, choosing ω to be the same for all particles and all pairs of layers, such that ω jlm ≡ ω for all j, l, m.However, it is important to point out that the inter-layer coupling could be different for different particles.For example, one could alternatively tune the relative value of ω for a given particle based on a particle property, such as total strength.Investigation of more complicated methods for choosing the inter-layer couplings may be an interesting direction for future work.
To proceed, it is necessary to understand the effect of ω on the community structure.There are two limiting cases which are relatively simple to grasp: when ω = 0, there are no connections between layers of the adjacency tensor, and we therefore expect to recover the results of static community detection (Fig. 3).At the other extreme, ω can be made large enough such that the strength of inter-layer connections entirely overwhelms the strength of intra-layer connections, resulting in completely consistent community structure across all compressive steps (that is, no observable changes; see Fig. 17 in the Appendix for an example partition at large ω).To understand what occurs in between these limiting cases, we consider a simple measure of network rearrangement called flexibility, or Ξ, previously defined in [62].The flexibility of a single particle i, ξ i , is defined as the number of times a particle changes in community allegiance across network layers, normalized by the number of possible changes.Mathematically, where g i is the number of times that the particle changes its community and L is the total number of strain steps.
The flexibility of the entire multilayer network is then given by the mean flexibility of all particles where N is the number of particles.
In order to choose a physically relevant value of the interslice coupling, we run 20 optimizations of multilayer community detection on the normal force network A n for each packing, for several values of ω between 0 and 1, in steps of ∆ω = 0.01.Note that ω = 0.01 corresponds to a coupling which is equal to 1/100 of the mean edge weight in a given layer, due to the normalization procedure described in Sec.III B 3. Due to the stochastic nature of modularity maximization [55], we avoid an interlayer coupling close to zero, as the community structure is more likely to be influenced by noise in the algorithm.At each ω, we compute Ξ for all optimizations of a given packing, and then average over optimizations to obtain a single value of flexibility for the packing, Ξ opt .In what follows, we will denote averages over optimizations as • opt , and averages first over optimization and then over packings with an overbar.
As shown in Fig. 5, we observe that the average flexibility Ξ decreases smoothly as the inter-layer coupling ω increases.As expected, at ω = 0, the network flexibility is unity, indicative of the fact that there is no consistency in community structure across layers.That is, at each step all particles are assigned to new communities which are independent from the community they were assigned to in the previous step.At ω = 1, which corresponds to an inter-layer coupling of the same magnitude as the mean edge weight in each layer, the flexibility is approximately zero, and there are no changes in the community structure; instead, we observe complete consistency of community structure throughout layers.
Interestingly, at the first ω value away from zero, Ξ sharply drops to ≈ 0.33.This steep decline can be explained by a recent result, which demonstrates that the case ω = 0 is singular in the sense that even when ω = with 0 < 1, there will be at least some persistent community structure [61].We thus take the value of Ξ obtained from the first small step away from 0 at ω = 0.01 to be the upper bound of network flexibility, Ξ max ; from that point forward, Ξ smoothly decreases to zero (see Appendix B for further considerations into this point).To further characterize the behavior of the flexibilty versus ω curve, we assessed whether the trend could be described as exponential decay, such that for each Ξ opt , we have Ξ opt ≈ Ae −ω/ωo .We find that the data can indeed be well approximated by an exponential.The ω o for all packings fall in the range 0.26 ≤ ω o ≤ 0.40, and represent characteristic inter-layer couplings for the system.
At ω values too far into the exponential tail, the communities will not be sensitive to the structure present within a given layer, and at small values of ω, temporal dependence becomes less important.We thus pick the optimal value of interslice coupling, ω * , to be the value such that Ξ(ω * ) is closest to half of the maximum flexibility, Ξ max .This procedure yields a value of ω that balances the trade-off between the importance of intra-layer edges (particle contact forces) and persistent structure across network layers.In Sec.IV C we further validate this choice by comparing the community structure obtained at ω * to three relevant null models, showing that the real network is distinguishable from the null models in each case.
Using the method described above, we find ω * = 0.16, (denoted by the red dot in Fig. 5), and we use this value in community detection for all packings and for both the normal and tangential force networks.The inset of Fig. 5 shows what the distribution of ω * would be if we were to optimize ω for each packing individually.We acknowledge that there are several other methods that could be used to determine an appropriate coupling, but here we have focused on a straightfoward method to choose a physically meaningful ω, which yields intermediate values of network flexibility.

A. Extraction of evolving mesoscale structure from the force network
To extract pressure dependent particle assemblies (communities) from the multilayer force networks of each experimental run, we maximize multilayer modularity (Eq.5).We consider both the normal and tangential force networks separately, with A ijl defined as in Sec.III A, and P ijl given by Eq. 6 (or 7).As determined in the previous section, in both cases we use γ = 1 and ω = ω * = 0.16.For each particle configuration, we carry out 20 maximizations of the multilayer modularity to obtain an ensemble of partitions, each with their respective value of Q multi .
By the nature of modularity maximization and our choice of null model, the resulting communities correspond to spatially localized, mesoscale groups of particles that display collective organization throughout the compression process.In Fig. 6 we show an example of the community architecture at each step above Φ J detected Choosing an optimal inter-layer coupling.The average network flexibility is governed by the value of the inter-layer coupling ω.Each dot is the average of network flexibility over all experimental packings, Ξ, vs. the interlayer coupling, ω.When ω = 0, Ξ = 1 and the results of static community detection are recovered; there is no persistent structure across compression.When ω = 1, Ξ ≈ 0 and the particle contact forces within each layer are overwhelmed by the weight of the inter-layer coupling such that the community structure exhibits complete consistency throughout the compressive steps.At ω = 0.01, Ξmax ≈ 0.3349.The optimal inter-layer coupling ω * = 0.16 is denoted by the red dot, and is chosen such that Ξ(ω * ) ≈ Ξmax.This ensures a balance between intra-and inter-layer weights.The inset shows what the distribution of ω * would be if the identical optimization were performed for each packing individually.
from the normal force network (A) and tangential force network (B ), of a particular experimental configuration.In both cases, communities are colored according to their multilayer modularity value, Q multi .Importantly, the same color at each strain step corresponds to the same community, to provide a visual sense of linking particle assemblies continuously throughout compression.Unlike in the single layer modularity maximization, the communities here can persist across layers, providing a means to directly examine network evolution and reconfiguration.
It is important to point out that the structure uncovered with multilayer community detection is representative of the system as a whole across all layers, and not necessarily of the static force chain structure of individual layers.Because of this, it is not required that a community in a given layer consists of only physically connected particles; rather communities in a given layer are dependent on the structure of the entire temporal network.In this way, the communities embody changing network architecture, and we can thus observe the break up and coalescence of communities across the applied compressive steps.By eye, the mesoscale organization appears to become more compact, and we observe the emergence of force chain-like organization as the packing fraction increases.As expected for a compression process, these structural patterns are most clearly evident for the normal force network.The communities from the tangential network appear visually smaller than those extracted from the normal forces.In the following sections, we study the reconfiguration and physical properties of this evolving architecture, and quantify differences between the normal and tangential mesostructures.

B. Characterization of mesoscale reconfiguration
We characterize the evolution of the multilayer community structure using two diagnostics: network flexibility, Ξ (Eq.9), and community stationarity, ζ c .Recall that the network flexibility quantifies the amount of reconfiguration in the system at the particle level, determined by the fraction of steps over which a particle changes its community allegiance.As a second measure of structural reconfiguration we also consider the community stationarity [56,63], which measures the consistency of particle content in each community throughout compression.
To define stationarity, we begin by writing the autocorrelation J(c l , c l+m ) between a given community at layer l, c l , and the same community at layer l + m, c l+m , as where |c l ∩ c l+m | is the number of particles present in community c at strain step l that are also present in community c at step l + m, and |c l ∪ c l+m | is the number of distinct particles present in community c at strain step l or step l+m.Then if l i is the layer in which community c first appears, and l f is the layer in which it last appears, the stationarity of community c is In this way, communities that experience large changes in their particle content over consecutive compressive steps will have larger values of ζ c than communities with more consistent structure.The average stationarity of a multilayer network is obtained by taking the mean of ζ c over all n c communities: where we exclude singleton communities, which only contain one particle in all layers.In Fig. 7 we show boxplots over the particle configurations of Ξ and ζ for the normal and tangential force networks.We have averaged these diagnostics over the 20 optimizations, to obtain a mean value of Ξ opt and ζ opt for each experimental packing.We test whether the reconfiguration of the normal and tangential force networks can be distinguished by performing non-parametric permutation tests on the flexibility and stationarity values.For the network of normal forces, there are 97 values of the flexibility Ξ n opt and stationarity ζ n opt (one for each laboratory configuration).Repeating the same protocol for the tangential forces yields another set of values Ξ t opt and ζ t opt .We calculate the difference in the means of these two distributions, and test whether that difference is greater than expected in the null distribution created by re-assigning statistics uniformly at random to the two groups: "normal" and "tangential".Using this test, we find significant differences in the means (over optimizations and then packings) of both statistics, with both p-values less than 1 × 10 −3 .In particular, Ξ n < Ξ t and ζ n > ζ t .From this finding, we conclude that at the same value of interlayer coupling, the multilayer network of normal forces tends to exhibit less reconfiguration during compression than the network of tangential forces.This sensitivity to differences in the evolution of two related but distinct force networks suggests that our method may be more broadly applicable.For example, it could be used to test for differences and classify different types of granular systems composed of varying particle materials, shapes, or sizes.

C. Reconfiguration of network architecture during compression: a null model comparison
Perhaps more important than absolute values of network measures is the question of whether or not the network evolution observed in the real, physical system is significantly different than what is expected from relevant null models, and whether or not our model is sensitive to these differences.In this section, we demonstrate that the multilayer community structure of the compressed granular configuration is indeed distinct from three null models with respect to the diagnostics defined previously (Ξ, and ζ).Here, we additionally consider the multilayer modularity, Q multi .(Recall that Q multi is a general measure of how well the network can be partitioned into densely interconnected groups of particles throughout the compression process, with respect to the physi- cally appropriate geographical null model).To ensure a fair comparison, we test the real force networks from each of the 97 experimental realizations against null models that are built using the force information from the same experimental run.Furthermore, we analyze the normal and tangential forces separately.To determine statistical significance, we perform permutation; assignments of statistic values to the two groups, "real network" or "null network", are permuted uniformly at random to construct a null distribution of expected differences be- A comparison of granular network evolution to three null networks.The community structure of the multilayer granular system is compared to three relevant null models (top row), using the quantities Ξ, ζ, and Q multi .For each null model, the boxplots show the normalized difference between the real and null network diagnostics for all experimental configurations, and for the normal (second row) and tangential (third row) components.Significant differences between the real and null networks exist when the boxplots are above or below the zero line.The p-values obtained from permutation testing are shown beneath each box to quantify the significance of results.(If no p-value is shown, the difference is not significant).tween the two distributions.
In the first case, we impose the elimination of steady perturbation on the system, comparing the real multilayer force networks to null models built by setting all layers to be equal (Fig. 8A).In particular, we construct a null network, A repeat , by repeating the real force network at a constant layer s, L times (recall that L is the total number of layers).Carrying out this process for all such constant layers yields a set of L null networks for each packing.We run community detection 20 times for each null network, and compute the network diagnostics (Q null multi , Ξ null , and ζ null ) for each trial.We average quantities first over the optimization realizations, and then over the different networks.Physically, this simple null model serves as a control that allows us to assess the implications and consequences of compression on the detected community structure.Since there are no changes in topological structure or edge weights from one layer to the next in the synthetic networks, we expect all changes in community structure to be due to noise.This baseline network reconfiguration should be much less than the reconfiguration that occurs in the real networks, which encode the compression procedure.Indeed, for both the normal and tangential force networks, we see that the null models have significantly lower flexibility, Ξ null < Ξ real , and higher stationarity, ζ null > ζ real , than the compressed system (see Fig. 8A,D,G).Note also that the modularity of the null model is greater than that of the true networks, which is also expected since the null models will be partitioned into highly temporally consistent community structure.
We next compare the real system to a null model, A temporal , constructed by permuting the temporal layers (strain steps) of the real networks uniformly at random (Fig. 8B ).In the experimental protocol, compression is applied systematically, in small and always increasing steps.The temporal null model considered here effectively eradicates this regularity in the layer ordering.The expectation is thus that the real networks will have less reconfiguration than the scrambled model.For each experimental packing, we create 20 null networks built from different random permutations of the layers, and run 20 optimizations of the Louvain algorithm for all permutations.As before, we compute Ξ null , ζ null , and Q null multi , and average the results first over optimizations and then over permutations.In Fig. 8E,H we compare the real system to the null model, finding that our method is again sensitive to differences in the evolution of the real and synthetic networks.For both the normal and tangential force networks, Ξ null > Ξ real , and ζ null < ζ real implying more steady and regular progression of community structure in the real system.In addition, we observe a slight decrease in modularity in the temporally permuted normal force network, suggesting that the real compression protocol yields stronger multilayer community structure.The modularity of the tangential forces is less affected by temporal scrambling than that of the normal forces.We also find that the real system can be distinguished from the temporal null model with respect to an alternative measure of reconfiguration known as promiscuity, Ψ [64].See Appendix C 2 for a definition of this statistic and the results of the null model analysis.
In the final null model (Fig. 8C ), we consider the spatial distribution of forces throughout the system.We construct a null model, A edges , by permuting the edge weights uniformly at random within each layer while maintaining the original contact topology and ordering of slices (for a related but distinct null model, see [28]).It is known that the organization of inter-particle forces is crucial for the stability of granular packings.This fact is manifested in force chains, branching groups of particles that bear the majority of the load in the system.Therefore, for the multilayer model to be useful, it should not be agnostic to the pattern of forces present in physically realizable systems.For each of the 97 configurations, we form an ensemble of null networks by permuting the forces uniformly at random within each layer 20 different times, and then optimize modularity 20 times for each permuted network.As before, this is done separately for the normal and tangential forces.In this case, we expect the synthetic networks to display more reconfiguration and less modular structure than the physical networks.We observe that the diagnostics of multilayer community structure are highly distinguishable between the real multilayer force networks and the null model (Fig. 8F,I ).In particular, the force-permuted networks exhibit more flexibility, Ξ null > Ξ real , less stationarity ζ null < ζ real , and decreased modularity Q null multi < Q real multi for both the normal and tangential components than expected in the null model.These results confirm our hypothesis that the evolution of network structure is less variable and undergoes less reorganization in the real system, and that the modularity stronger.In the Appendix C 2, we additionally show that these results hold for the promiscuity.
The findings presented in this section crucially demonstrate that the multilayer network model and community detection are sensitive to differences in the structural evolution of the compressed system compared to three relevant null models, with respect to Ξ, ζ, and Q multi .In particular, we can distinguish between the evolution of the granular networks and null models with consistent topology and weights, scrambled temporal layers, or scrambled edge weights (but same topological structure).Even more importantly, the differences in the reconfiguration agree with what is expected from a physical standpoint, making evident the powerful utility of this framework.

D. Physical properties of evolving mesostructures
In the previous section, we examined three measures to quantify the community structure of the multilayer granular network.However, the diagnostics we considered did not directly describe the evolution of physical properties of the mesoscale organization.We turn now to a physical description of the network architecture, and define measures to quantify both the strength and geometry of community structure throughout compression.We then ask whether the physical properties of mesoscale particle assemblies can be related to the measures of structural reconfiguration previously defined.

Evolution of community characteristics throughout compression
We characterize the physical nature of community structure with three measures: size, strength, and sparsity, and examine how these quantities evolve throughout compression.Similarly to the previous section, we consider and compare both the normal and tangential force network.The size s c l of community c at strain step l is simply the number of particles within the community at that step.We denote the strength of a community at layer l as σ c l , and define it to be the average amount of force (normalized by the mean force in the layer) on a particle in community c from intra-community contacts.Mathematically, where A is either the normalized normal or tangential force network.Finally, we consider a measure of spatial compactness, which we term the community sparsity, η c .The sparsity of a community is closely related to the hull ratio, which has been defined and used to quantify the geometric arrangement of compressed particle assemblies [65].The hull ratio h c l of a community at strain step l can be understood as the ratio of the area a of particles in the community to the area of the convex hull of the group a hull , such that where the area of particle i is a i = πr 2 i , with r i equal to the particle radius.We then take the community sparsity at layer l to be With this definition, small values of η correspond to spatially dense groups of particles, and high values of η correspond to sparse particle configurations.We now use each of these physical characteristics to quantify the evolution of the mesoscale community structure throughout the compression process.Given a partition of particles into communities, we compute s c l , σ c l , and η c l at each strain step l, for all communities except those which have s c l = 1 for all l (i.e., they are always singletons).We then define s l , σ l , and η l to be the average over all communities present in layer l, where n c l is the number of communities present in layer l.
Repeating this process for all community detection optimizations and all particle configurations, we form representative curves of each physical quantity as a function of strain step by averaging the measures defined above first over optimizations and then over experimental configurations.We denote the final averaged physical quantities for size, strength, and sparsity as s , σ , and η , respectively.
We first observe that the size of the mesoscale structure increases smoothly as a function of the strain steps, (Fig. 9 A), suggesting that the scale of the mesostructure organization grows with compression.However, the normal force communities are noticeably larger in size, and undergo relatively more growth throughout strain steps than the communities from the tangential network.These differences imply that the network of normal forces exhibits collective organization on a larger scale than that of the tangential force network; and that the normal force network responds differently to increasing pressure.Some of these features are partially recognizable by eye in comparing the community structure in Fig. 6A and Fig. 6B.We also find that the community strength increases smoothly over the applied strain steps, for both the normal and tangential components.This behavior signifies the mesoscale architecture becoming more strongly connected throughout compression, which agrees with the physical expectation.Note that the average community strength σ was computed on the normalized networks (see Sec. III B 3), which explains the similar scale between the normal and tangential curves.Finally, we observe a slight decrease in community sparsity across strain steps, especially during the beginning stages.In addition, the tangential network displays lower sparsity than that of the normal force network, implying more compact tangential community structure.This feature is likely tied to the smaller size of tangential communities, which constrains the set of possible spatial arrangements of the particles within a community.

Trends in physical characteristics are diverse
In addition to quantifying the average behavior of mesoscale architecture, it is also important to investigate the behavior of individual communities throughout compression.Though it is possible that all communities progress similarly to the average behavior of the system (for example, coalescing to create communities of increasing size and strength, but decreasing sparsity), this does not have to be the case.We find, in fact, that the situation is quite the opposite; at the level of single communities, the evolution of physical structure varies greatly.We demonstrate this in a simple way.First, we identify the number of communities that exhibit linear trends with respect to size, strength, or sparsity as a function of strain step.Tables I and II show the results of this analysis for the normal and tangential networks.We observe that the majority of mesostructures do not exhibit consistent and predictable linear trends in terms of their physical properties throughout the compression process.While some linear tendencies are much more likely to occur than others (for example, increasing strength and decreasing sparsity), the behavior of many communities cannot be characterized by a simple linear relationship.This surprising results highlights the important diversity of mesoscale structural evolution.The average community size s of the normal and tangential networks increases with strain step, but the tangential force communities tend to be smaller than the normal force communities.(B) The community strength σ measures the intra-community forces, and also grows as a function of strain step for both the normal and tangential networks.(C) The sparsity η quantifies the spatial density of community structure.The tangential force communities tend to be more dense than the normal force communities.In all plots, error bars correspond to the standard error of the mean over packings.Next, we ask if and how the set of communities which do have linear behavior with respect to a given physical property, are related to each other.In Fig. 10, we plot the slope of the linear regression fit of sparsity vs. the slope of strength, for each community in all optimizations and experimental packings.Again, the scatter plots point to the variation of mesostructure development, as all quadrants (with the exception of the upper left) are significantly filled in.These data support the notion that communities may coalesce, disband, or become increasingly branch-like, and each of these behaviors is observable as the packing rearranges under compression.
To quantify the co-dependence of the two statistics, we first find the number of communities (i.e. the intersection) that fall within each quadrant for each optimization and packing.For example, if σ ↑ are the communities with linearly increasing strength and η ↑ are those with linearly increasing sparsity, then we compute the number of communities that satisfy σ ↑ ∩ η ↑ (upper right quadrant of the scatter plots).Then to determine how often increasing strength (σ ↑ ) occurs with increasing sparsity (η ↑ ), for example, we normalize the intersection by the total number of communities with σ ↑ .Conversely, if we want to know the percentage of communities with increasing sparsity (η ↑ ) that also have increasing strength (σ ↑ ), then we would instead divide the intersection by the number of communities in η ↑ .Table III in Appendix C 1 shows the percentages for each of the possible combinations.The strongest relationship occurs between communities with σ ↓ and η ↓ .In this case, we find that if a community has linearly decreasing strength with strain step, then it is likely to also be more compact (but note that not many communities decrease in strength in the first place).These results may be due to the community losing particles and thus becoming more dense (on average for the normal (tangential) network, 98.6% (96.6%) of communities with linearly decreasing strength and sparsity also have linearly decreasing size).The nearly empty upper left quadrants are consistent with this relationship as well; communities with decreasing strength rarely become more spatially spread out throughout compression.We also find that on average, more than half of the communities with increasing sparsity also have increasing strength, which is likely due to the community gaining particles, thus allowing it to take on configurations which are more spatially extended (on average for the normal (tangential) network, 93.9% (91.7%) of communities with linearly increasing strength and sparsity also have linearly increasing size).linear trend in community strength linear trend in community sparsity FIG. 10.Scatter plots demonstrating the diversity of network evolution.The relationship between the slopes of communities that exhibit linear trends with respect to strength and sparsity.For both the normal (A) and tangential (B) force networks, all but the upper left quadrant are quite populated, pointing to the diversity in the evolution of community structure throughout the compression process.

E. Linking physical properties to network reconfiguration
Thus far, we have independently characterized the evolution of multilayer community structure using notions of particle rearrangement (flexibility and stationarity), and using physical quantities, (size, strength, and sparsity).We now attempt to link these two ideas together, asking whether network reconfiguration can be related to physical aspects of the network.

Local reconfiguration
We first investigate the relationship between particle flexibility, ξ (Eq.8) and the inter-particle force, f .Recall that ξ is a measure of local reconfiguration in that it is defined for a single particle, but it is determined from the mesoscale community structure of the network.For every multilayer community partition (20 for each particle configuration), we compute the flexibility of each particle as given in Eq. 8, and average these values over partitions.This yields a single value of flexibility ξ i for the i th particle in a given experimental run.We do this for the normal and tangential force networks separately.
Our first finding is that flexibility ξ is strongly correlated with the average force on a particle throughout strain steps f φ , as well as the average absolute change in force on the particle |∆f | φ between consecutive strain steps.This result holds for both the normal and tangential components.For the i th particle, the average change in force |∆f i | φ is given by where f i l is the total force on the i th particle in layer l, determined from the adjacency tensor as f i l = j A ijl .In particular, we observe that particles with high flexibility ξ also tend to have high values of average force f φ and average change in force |∆f | φ .In Fig. 11, we plot ξ vs. f φ for each particle using the normal (A) and tangential (C) force networks of one experimental configuration.Panels (B) and (D) show ξ vs. |∆f | φ .We quantify these relationships for all experimental packings using the Spearman's rank correlation, ρ.For the normal forces, the average correlations over packings for ξ vs. f φ and ξ vs. |∆f | φ are ρ f = 0.81 and ρ ∆f = 0.74, respectively, with all p-values satisfying p f < 1 × 10 −174 and p ∆f < 1 × 1 −127 , respectively.(In Fig. 21 A,B of the Appendix, we show the distributions of ρ f and ρ ∆f for all packings).For the tangential forces, ρ f = 0.80 and ρ ∆f = 0.71, with all p-values satisfying p f < 1 × 10 −181 and p ∆f < 1 × 10 −109 .(In Fig. 21 C,D of the Appendix we show the distributions of ρ f and ρ ∆f ).In addition to the flexibility, we also tested the relationship between force and reconfiguration on a more robust measurement of local network rearrangement called promiscuity [64], finding that the relationship still strongly holds.See Appendix C 2 for a description of the promiscuity statistic, an example scatter plot, and correlation values.
To understand these results, first recall that the flexibility of a particle is a measure of how strongly fixed the particle is to its given community; ξ is the number of times a particle changes community normalized by the number of possible changes, so lower values of ξ correspond to particles that have more stable community allegiance across network layers.Our finding thus implies that large forces (or large changes in forces) are associated with particles reconfiguring and shifting communities throughout compression.Force drives local reconfiguration.(A, C) Scatter plots of particle flexibility ξ vs. the average force on a particle across compression f φ for a sample packing.In both the normal and tangential networks, there is a strong, positive Spearman correlation between the two quantities.(B, D) Scatter plots of particle flexibility ξ vs. the average absolute change in force on a particle across compression |∆f | φ for a sample packing.In both the normal and tangential networks, there is a strong, positive Spearman correlation between the two quantities.See Fig. 21 for the distribution of correlations for each packing).Community reorganization is driven by changes in intra-community force.Scatter plots show the community stationarity, ζc, vs. the average absolute change in community strength across compression, |∆σc| φ , for an example packing.For both the normal (A) and tangential (B) networks, there is a strong Spearman correlation between the two quantities (see Fig. 22 for the distribution of correlations for each packing).

Mesoscale reconfiguration
We now examine the relationship between mesoscale reconfiguration as measured by the stationarity, and changes in physical properties of the multilayer community structure as measured by strength and sparsity.In particular, we first ask if reconfiguration at the community level is correlated with changes in the community strength throughout the compression process.We next investigate if changes in community particle composition drive changes in the spatial arrangement of community structure.For a given experimental run and modularity optimization, the stationarity (ζ) for all communities is computed using Eq.11, and the strength (σ) and sparsity (η) of all communities at each layer are computed using Eqs.13 and 15, respectively.We then compute the Spearman correlation between stationarity, ζ c and the average absolute change in strength |∆σ c | φ and sparsity |∆η c | φ across all layers in which a given community exists.The average absolute change in the strength of a community c across applied strain steps is calculated according to and the average absolute change in sparsity is given by Communities that have s c l = 1 for all layers are ignored.We find that the stationarity is significantly anticorrelated with changes in community strength (a measure of intra-community force).For both the normal (A) and tangential (B) networks, there is an increased correlation between the two quantities across packings.
is ρ ∆σ = −0.65 with all p-values less than 0.0033, and for the tangential network, ρ ∆σ = −0.69 with all p-values less than 1 × 10 −13 .Fig. 22 A,C in Appendix C 3 show the distributions of the optimization-averaged Spearman correlations for each experimental packing.
To examine the effect of mesoscale network reorganization on community spatial structure, we show an example of |∆η c | φ vs. ζ c for one packing in Fig. 13.Calculation of Spearman correlation values between the two statistics identifies a strong negative correlation between stationarity and sparsity, suggesting that changes in community particle content (ζ) indeed lead to changes in spatial density (η).We find ρ ∆η = −0.77with all p-values less than 1 × 10 −10 for the normal force network, and ρ ∆η = −0.75 with all p-values less than 1 × 10 −23 for the tangential forces.Fig. 22 B,D in Appendix C 3 show the distributions of the optimization-averaged Spearman correlations for each experimental packing.
The findings presented here demonstrate first that changes in community strength (i.e., in intra-community force), are associated with mesoscale reconfiguration across compression (i.e., the stationarity).This result complements those of Sec.IV E 1, where high values of force and change in force across compression were associated with high values of local reconfiguration -a similar relationship holds for reconfiguration at the intermediate scale as well.We have additionally shown that if a community experiences large changes in its particle content, then the spatial arrangement of the community's particles undergoes restructuring as well, agreeing with physical intuition.It is thus clear that the reorganization of mesoscale network architecture is strongly tied to changes in physical properties of the mesoscale structure that occur due to the compression process.
It is important to point out that different communities might not exist for the same number of strain steps.Therefore, as a more robust measure of the relationship between community strength and stationarity, we repeat the above analysis but consider only those communities that persist throughout the entire compression cycle.If we define the lifetime of a community to be the number of steps in which it exists divided by the total number of steps in the network, then we consider the set of communities that have lifetimes equal to 1. On average over all optimizations and experimental runs, these long-lived communities correspond to fractions of ≈ 0.62 for the normal force network and ≈ 0.54 for the tangential force network.For both the normal and tangential forces, we find that the correlations between community strength and stationarity increases compared to the correlation calculated using communities with all lifetimes.In particular, ρ ∆σ = −0.80 with all p-values less than 1 × 10 −8 for the normal force network (see Fig. 14A for an example) and ρ ∆σ = −0.78 with all p-values less than 1×10 −12 (see.Fig. 14B for an example) for the tangential force network.(Note that the horizontal "stripes" at high values of ζ correspond mostly to communities with only a few particles; in the case of small groups, it is likely for different communities to achieve the same value of stationarity, because only so many reconfigurations are possible).The stronger correlations observed here may be due to the fact that they are inferred only from communities that exist for the same number of steps (and so are more comparable to one another), and also from the fact that these communities have the most data over which to compute averages.Physically, these longer-lived communities are also expected to undergo more continuous re-arrangements with compression, thus offering smooth estimates of stationarity; while shorter-lived communities may experience more extensive rearrangements, leading to their birth or death, and thus may offer more variable estimates of stationarity.Fig. 23 in Appendix C 3 shows histograms of the Spearman correlations for all packings.

F. Sensitivity to the subsystem
A useful feature of our experimental setup is the presence of the low friction particles in the center of the packing.As a final demonstration of the sensitivity of the multilayer network model to the underlying physics, we examine its ability to detect differences in the physical properties of the subsystem (low µ) vs. the bath (high µ).Sensitivity to such variation is crucial to the utility of this framework.
The first step is to consistently define the subsystem.The majority of subsystem particles are identified using UV dye, as described in Sec.II.We use the UV dye image to initially identify particles as low friction or high friction particles at each strain step.Then, after having tracked particles in each frame, particles are labeled as low friction if they are identified as low friction for more than half of the steps.This technique ensures that low friction particles are not mislabeled as high friction particles, and the reverse rarely occurs.
The second step is to define what constitutes a subsystem community vs. a non-subsystem community.Since the communities we consider consist of more than one particle, it is highly unlikely for all particles within a group to either be only inside the subsystem radius or only outside the subsystem radius.We therefore consider a subsystem community, c s , to be a community such that more than 75% of its particles are contained within the low friction area at some point during the compression process.Any community c b that does not meet this condition is considered part of the bath.To ensure a fair comparison between these two groups, we further require that the communities examined in the two groups must be of similar sizes.(Due to the nature of the experimental setup, the subsystem communities are localized within a smaller region of space, and therefore will have a larger constraint on their possible sizes).Within each optimization, we compute the average size s c φ of all communities across all steps in which the communities exist.We then find the largest and smallest average sizes of the subsystem communities, s cs φ,max and s cs φ,min , and find all bath communities that have sizes s c b φ such that s cs φ,min ≤ s c b φ ≤ s cs φ,max .This protocol mitigates the likelihood that observed differences are due to differing community sizes.
We now proceed to investigate if and how the community strength and sparsity evolve differently for the subsystem mesostructures vs. those in the bath.We follow the procedure described in Sec.IV D 1, but here we separate the subsystem and bath communities.In Fig. 15, we show the results of this analysis for both the normal (A-C) and tangential (D-F) force networks.For each of the three statistics previously considered, s , σ , and η , we plot one curve for the bath (shown in black) and one for the subsystem (shown in red).Panels (A) and (D) show the evolution of community size; here we note the scale of the y-axes, which show that we are indeed comparing c s and c b of similar sizes (the curves remain within a size difference of about 1 particle throughout all steps).
We next examine the intra-community force of the high and low friction groups.Panels (B) and (E) show the average community strengths of the subsystem and bath.For the normal force networks, the two curves are very close; in other words, there does not appear to be a large difference in the evolution of the intra-community normal forces between the low and high friction mesostructures.This is not a surprising result, however, as we expect the difference in friction to more significantly effect the tangential component rather than the normal component.Indeed, when we consider the network of tangential forces, we are clearly able to distinguish between the architecture of the subsystem and bath.At the beginning of the compression cycle, the average strength of the low and high friction communities still follow each other closely, but after the third step, we observe that the two groups separate.In addition to the community strengths being different, the way in which they differ Curves showing the average community size of the subsystem (red curve) vs. bath (black curve) communities for the normal and tangential networks.The small differences between the curves for each network ensures that we compare communities of similar size between the two groups.(B, E) The evolution of community strength in the network of normal forces is not significantly different between the subsystem and bath.However, the tangential network displays a clear separation, suggesting that the multilayer network model and community detection are sensitive to differences in friction.More importantly, the direction of the sensitivity agrees with physical intuition; the modular structures found in the low friction bath are characterized by lower tangential intracommunity force.(C, F) There is a small distinction between the community sparsity of the two groups, with the low friction communities being slightly more dense throughout compression.
agrees with physical intuition; the intra-community tangential force of the high friction particle assemblies grows significantly more than the intra-community force of the low friction particle assemblies as compression continues.
We also consider the average community sparsity in order to quantify differences in the spatial arrangement of the mesostructures.We observe that η tends to be smaller for the subsystem communities than those in the bath, although the differences are quite small.It is thus difficult to make a strong conclusion from the sparsity measure alone.
The findings discussed here first indicate that our model is indeed sensitive to underlying differences in the bath and the subsystem.But perhaps more crucial are the physical conclusions about network evolution; if only the network of normal forces is considered, it is difficult to discern the evolution of the low and high friction mesostructures.In this case, both groups display strong community organization throughout all applied strain steps.However, when the tangential force network is considered, the low friction subsystem forms noticeably weaker communities than the high friction bath, and the evolution of the two groups can be distinguished.This result suggests that the examination of both the normal and tangential force networks can lead to a more complete understanding of mesoscale structure in granular systems.In Sec.V, we discuss broader implications of the findings discussed in this section.

V. DISCUSSION AND METHODOLOGICAL CONSIDERATIONS
Although it is thought that mesoscale structures underlie many bulk properties of granular materials, there is a pressing need to develop theoretical tools for identifying these structures and describing their dynamics.Graph theoretic approaches provide a framework in which to address these challenges, via the network of contact forces connecting the grains.In this study, we use a common geometry (biaxial compression) as a case study for developing novel approaches drawn from network theory.
We describe the evolution of the granular system by modeling it as two temporal multilayer networks (one constructed from the normal component of the interparticle forces, and the other from the tangential component).Then, we extract inherent mesoscale structure from each network using multilayer community detection techniques with a physically grounded null model.Unlike in a static network, the particle assemblies we detect persist across layers, allowing a direct study of the evolution of mesoscale organization.The particle assemblies are quantified with measures of network reconfiguration, as well as with measures of physical architecture.We find that both the normal and tangential components display community structure, though the network of normal forces shows organization on a larger scale (bigger communities), and exhibits less restructuring during the sequence of applied strain steps.We also demonstrate that the community architecture of the real system is significantly distinguishable from three physically-relevant null models.
In addition to characterizing the evolution of the granular network in terms of force and the spatial arrangement of particle assemblies, we show that network reconfiguration can be tied to these physical descriptors.This is a crucial finding, as it bridges the gap between seemingly abstract measures of network reorganization and physical properties of network components.One re-sult in particular that merits discussion is that particle flexibility is positively correlated with the average force on a particle across layers.Although the notion that force drives local reorganization is intuitive from a physical standpoint, it is interesting to consider this result in the context of a different type of network where flexibility has also been studied.Since we consider particles to be nodes and forces to be weighted edges, there is a direct analogy between the total force on particle i in a given static network, and the weighted degree of node i, k i = j A ij in a general undirected network with adjacency matrix A. In a study on the dynamic community structure of functional human brain networks during learning, a seemingly opposite relationship between flexibility and node strength exists [66].It was found that nodes from densely connected brain regions exhibited little change in community allegiance throughout learning, whereas weakly connected nodes had higher values of flexibility.This points to an interesting contrast that exists between the evolution of the granular network and the much less physically constrained functional brain network.It is reasonable to think that the source of the opposing relationship between node strength and flexibility is due to the physical nature of the granular system, and the interpretation of node strength as force, which drives reconfiguration.But perhaps more importantly, this finding motivates the general need to develop more physically informed network measures to study these mechanical systems.The use of the geographic null model (vs. the traditional Newman-Girvan model) is an important step.
In the final section, we consider the presence of the low friction subsystem, asking if our model is sensitive to the subsystem, and more importantly, how.We find that from the normal force network alone, it is difficult to distinguish the two groups.However, consideration of the tangential network shows a clear difference in the evolution of the intra-community forces of the subsystem versus the bath, with the low friction modules characterized by noticeably smaller tangential strength.These conclusions illustrate that low-friction systems can have different community organization, and that the both the normal and tangential networks contain useful information.An interesting direction for future work would be to examine the evolution of sheared systems using the multilayer network approach, where we might expect to see heightened sensitivity to differences in friction, and perhaps altered reconfiguration patterns.
It is our hope that the method and results presented in this work will inform future studies on granular materials.For example, an interesting question is whether or not measures of network reconfiguration can be related to other notions of failure and deformation, like soft spots [67] and force chain buckling [68,69].Better yet, can we predict these rearrangement events from structural features of the network before deformation?Investigation of the latter will likely require the formulation of additional physically informed statistics of network dynamics, and may moreover be complemented by machine learning techniques.The fact that the multilayer network framework is sensitive to differences in the reconfiguration of null models as well as to differences in friction, further suggests that the framework may find applicability beyond this study.One can imagine using this formulation to not only characterize the evolution of a single system, but also to distinguish, quantify, and classify dissimilarities between several systems that are slightly different.For example, compressed systems vs. sheared systems, or experiments vs. simulations.One could also use this framework to examine how varying particle shapes and sizes effect the evolution and structure of the force network at the intermediate scale, or to inform the design of materials that exhibit tailored physical properties under external perturbations.Finally, the ability to identify and quantify mesoscale architecture is important in other physical systems more broadly.The brain [56,62,[70][71][72] and biomaterials [73][74][75][76][77] are a few examples of systems where intermediate structure is crucial to function.
In order to solidify the utility of the model presented here, it will be necessary to carry out a deeper exploration of the methods.The investigation of different null models is one place to start, especially in regard to the tangential network.In this study we took the absolute value of the tangential forces in order to use the geographic null model and directly compare to the normal forces, but in different situations (sheared systems, for example) the sign of the force is important and should not be discarded.However, to keep this information will require a null model that can correctly deal with signed adjacency matrices without losing physical meaning.A more thorough sampling of the phase space formed by the two resolution parameters, ω and γ, would also be useful.It will be interesting to consider the robustness and clarity of results with respect to changes in these quantities.Additionally, it may be more physically reasonable to move away from a scalar temporal coupling and examine more complex (but perhaps more informative) ways of linking layers together.For example, particle properties or similarity measures between the individual static networks could be used to determine the magnitude of coupling strengths in the multilayer formulation.Finally, it will be important to continue developing relevant measures of multilayer network reconfiguration that are informed by the physical questions we are attempting to answer.
In this study, we have utilized multilayer networks to model and characterize the structural reconfiguration that occurs in a compressed granular material.The framework introduced here is sensitive to the important inhomogeneity present in granular media, and allows for the direct extraction of evolving, intermediate scale structure that has not previously been probed or examined.We conclude by noting that in addition to particulate systems, the methods developed and exercised in this work can be applied more broadly across physical, biological, technological, and social systems to understand how dynamic interactions between system components organize into collective structure, which in turn constrains bulk properties and system performance.On the following page, we provide additional examples of community structure.Fig. 16 shows the results of community detection on the tangential force networks when the inter-layer coupling is zero; in this case, we recover the results of single layer community detection.As was the case with the normal force communities at ω = 0 (Fig. 3), the communities at each step are completely independent of one another.This fact motivates the generalization of the network framework to the multilayer regime, which allows for direct observation of evolving mesostructures throughout the compression cycle.Fig. 17 shows an example of community structure detected at a high inter-layer coupling (ω = 1), for the network of normal forces.In this case, the particle assemblies are almost completely consistent across layers, and we thus do not probe network reconfiguration.
Appendix B: Effect of small omega on network flexibility In Sec.III B 4 of the main text, we began the sweep over inter-layer couplings at ω = 0.01, which corresponds to a weight equivalent to 1 100 of the mean edge weight in each layer.Although this value is a relatively small step away from ω = 0 (where the average network flexibility is Ξ = 1), we find that even at such a small coupling, the flexibility drops significantly, to Ξ ≈ 0.3349.Here, we consider whether or not the flexibility can be brought closer to unity by bringing the coupling arbitrarily close to zero.
To test the nature of network dynamics at small coupling, we carry out 20 optimizations of modularity maximization for each packing at ω = 1 × 10 −200 (using the network of normal forces).The flexibility of each optimization is then calculated according to Eq. 9, and a single value for each configuration is obtained by averaging the values over the 20 realizations of community detection.Finally, we compute the average flexibility values over packings to obtain Ξ.Interestingly, we find that in the case of these granular networks, the mean flexibility over experimental configurations still remains much less than one, with Ξ(ω ) = 0.37.Fig. 18 shows the distribution of the optimization-averaged flexibility values, Ξ opt for each packing.This finding suggests that even arbitrarily small values of the inter-layer coupling do not significantly increase the amount of reconfiguration observed in the resulting community structure, and solidifies the validity of our choice of initial ω at 0.01.In Sec.IV D 2 of the main text, we examine communities that exhibit linear trends with respect to strength, σ and sparsity η, finding that communities with decreasing strength very often also have decreasing sparsity, and about half of communities with increasing sparsity also have increasing strength.There are many other possible combinations one could consider; here, we report percentages that quantify the likelihood that any two such combinations of linear trends occur together in a given community (III).To understand the following table, first recall that the quantity σ ↑ ∩ η ↑ , for example, is the number of communities with increasing strength and increasing sparsity.Then, the quantity (σ ↑ ∩ η ↑ )/σ ↑ gives the fraction of all communities with linearly increasing strength that also have linearly increasing sparsity.
The percentage of communities that progress together in terms of their linear trends with respect to strength σ and sparsity η.The first column denotes the relationship that we consider.The second and third columns are the average percentages over optimizations and packings for the given relationship, for the normal and tangential networks, respectively.Errors are the standard errors of the mean over packings.

Promiscuity as a more robust measure of local reconfiguration
Throughout this work, we use the notion of flexibility (Sec.III B 4) as a measure of local reconfiguration that is determined from the underlying mesoscale community structure in the system.As a more robust measure of local network reorganization, we also consider the particle promiscuity.The promiscuity ψ i of particle i is defined as the fraction of all communities in which the particle participates at least once, across all network layers.In the context of the current study, the promiscuity clarifies whether a particle is simply bouncing back and forth between the same two communities (which would yield high ξ but low ψ), or truly participating in many different communities throughout compression (a situation that leads to both high ξ and high ψ).As with the flexibility, we define the network promiscuity Ψ to be the average over all particles, We first test whether the promiscuity differs between the multilayer force networks of the real system, and the temporal (A temporal ) and force-scrambled (A edges ) null models.In both cases, we follow the analysis described in Sec.IV C, and obtain results in line with those found using the flexibility.In particular, for both the normal and tangential force networks, the promiscuity values of the real system are lower than the promiscuities of both null models.In Fig. 19 we show boxplots over the experimental configurations of (Ψ real − Ψ null )/Ψ real .(For each configuration, the value we consider is the average over community detection optimizations).The promiscuity values of the real and null networks are significantly different when the boxplots do not cross the zero-line.Boxplots showing that Ψ real < Ψ null for the normal and tangential networks, where the null model, A edges , is constructed by scrambling the edge weights within each layer while preserving contact structure.In all cases, statistical significance is determined by permutation testing.
We also investigated whether or not the strong positive correlation between inter-particle force and flexibility also holds for inter-particle force and promiscuity.Carrying out the same analysis as described in Sec.IV E 1, we find that the result is robust: high values of the average force or absolute change in force across compression are positively correlated with high values of ψ.In particular, we find ρ f = 0.81 and ρ ∆f = 0.74, with all p-values satisfying p f < 1 × 10 −173 and p ∆f < 1 × 10 −128 for the network of normal forces, and ρ ∆f = 0.79 and ρ ∆f = 0.70, with all p-values satisfying p f < 1 × 10 −169 and p ∆f < 1×10 −107 for the network of tangential forces.We show an example scatter plot of these relationships for the particles in a single experimental packing in Fig. 20, and show the distribution of Spearman correlations for each packing in Fig. 24.

Spearman correlation distributions
Throughout this work, we quantify relationships between network reconfiguration and physical properties of network structure using Spearman correlations and example scatter plots.In the following figures, we show histograms of the Spearman correlation values for these relationships for all experimental packings.In particular, we show flexibility vs. force (Fig. 21), stationarity vs. community strength (Fig. 22 A,C and Fig. 23), community sparsity vs. stationarity (Fig. 22 B,D), and promiscuity vs. force (Fig. 24).In all cases, we show the correlation values for both the normal and tangential networks.
FIG. 1. Experimental setup.(A) A schematic of the apparatus showing two walls bi-axially compressing an array of disk-shaped particles composed of an outer subsystem (black, high µ) and an inner subsystem (red, low µ).(B) The first image is taken with unpolarized white light and is used to locate particle positions.(C) The network of the photoelastic stress pattern, visualized with a second image taken with polarized light.This image allows for the calculation of the normal and tangential forces at each contact.A third image (inset) taken with ultraviolet light is used to identify the subsystem particles, which are tagged with fluorescent ink.

FIG. 2 .
FIG.2.Evolution of the force network.Photoelastic images of the system, taken at a sequence of eight steps above ΦJ ≈ 0.7825.The corresponding network of (normal) forces is overlaid on top of each picture.The collection of these static networks may together be represented as a single multilayer network that captures the evolution of the granular system throughout compression.
FIG. 4.A schematic of a multilayer network with layer-dependent community structure.Each layer represents a static granular force network in which nodes (particles) are connected to one another via intra-layer weighted edges.These weights can be either the normal contact forces or the absolute value of the tangential contact forces.Additionally, the same particle in consecutive layers is linked to itself with an inter-layer coupling, ω.For clarity, we only show two such couplings, but these inter-layer edges exist between all particles (and across all layers).Evolving communities are extracted from the multilayer network to study the mesoscale organization in the system, and to understand how it reconfigures due to the compressive cycle.The community structure can be determined from the network by maximizing the multilayer modularity, Q multi (Eq.5).In this schematic, the particles belonging to different communities are labeled by different colors.Note that the same community can persist across all layers and reconfigure in terms of particle content and strength throughout compression.
FIG. 5.Choosing an optimal inter-layer coupling.The average network flexibility is governed by the value of the inter-layer coupling ω.Each dot is the average of network flexibility over all experimental packings, Ξ, vs. the interlayer coupling, ω.When ω = 0, Ξ = 1 and the results of static community detection are recovered; there is no persistent structure across compression.When ω = 1, Ξ ≈ 0 and the particle contact forces within each layer are overwhelmed by the weight of the inter-layer coupling such that the community structure exhibits complete consistency throughout the compressive steps.At ω = 0.01, Ξmax ≈ 0.3349.The optimal inter-layer coupling ω * = 0.16 is denoted by the red dot, and is chosen such that Ξ(ω * ) ≈ Ξmax.This ensures a balance between intra-and inter-layer weights.The inset shows what the distribution of ω * would be if the identical optimization were performed for each packing individually.

FIG. 6 .
FIG.6.Evolving community structure of a compressed granular packing.For both the normal (A) and tangential (B) networks, particles are shown at their actual locations in physical space, and are colored according to their community assignment.Redder colors correspond to communities with higher multilayer modularity.The packing fraction increases from left to right, and only the steps with Φ > ΦJ are shown.Visually, the normal force communities appear to become more compact throughout compression, and the emergence of force chain like structure is observed.The tangential network forms noticeably smaller groups of particles.We quantify these physical differences, as well as differences in network reconfiguration, throughout the main text.

F
FIG. 8.A comparison of granular network evolution to three null networks.The community structure of the multilayer granular system is compared to three relevant null models (top row), using the quantities Ξ, ζ, and Q multi .For each null model, the boxplots show the normalized difference between the real and null network diagnostics for all experimental configurations, and for the normal (second row) and tangential (third row) components.Significant differences between the real and null networks exist when the boxplots are above or below the zero line.The p-values obtained from permutation testing are shown beneath each box to quantify the significance of results.(If no p-value is shown, the difference is not significant).(A) A null model formed by repeating the same layer for all steps, A repeat .The values Ξ, ζ, and Q multi are statistically different between the real and null networks for the normal (D) and tangential (G) components.(B) A null model, A temporal , constructed by permuting the temporal layers uniformly at random.The values Ξ, ζ, and Q multi are statistically different between the real system and null model for the normal forces (E), and Ξ and ζ are statistically different between the real system and null model for the tangential forces (H).(C) The third null model, A edges , is built by permuting the edge weights uniformly at random, while maintaining the original contact topology and temporal ordering.The values Ξ, ζ, and Q multi are statistically different between the real and null networks for both the normal (F) and tangential (I) components.
FIG. 8.A comparison of granular network evolution to three null networks.The community structure of the multilayer granular system is compared to three relevant null models (top row), using the quantities Ξ, ζ, and Q multi .For each null model, the boxplots show the normalized difference between the real and null network diagnostics for all experimental configurations, and for the normal (second row) and tangential (third row) components.Significant differences between the real and null networks exist when the boxplots are above or below the zero line.The p-values obtained from permutation testing are shown beneath each box to quantify the significance of results.(If no p-value is shown, the difference is not significant).(A) A null model formed by repeating the same layer for all steps, A repeat .The values Ξ, ζ, and Q multi are statistically different between the real and null networks for the normal (D) and tangential (G) components.(B) A null model, A temporal , constructed by permuting the temporal layers uniformly at random.The values Ξ, ζ, and Q multi are statistically different between the real system and null model for the normal forces (E), and Ξ and ζ are statistically different between the real system and null model for the tangential forces (H).(C) The third null model, A edges , is built by permuting the edge weights uniformly at random, while maintaining the original contact topology and temporal ordering.The values Ξ, ζ, and Q multi are statistically different between the real and null networks for both the normal (F) and tangential (I) components.

FIG. 9 .
FIG.9.The evolution of physical characteristics of community structure throughout compression.Community structure is characterized by size, strength, and sparsity for the normal (blue) and tangential (orange) networks.(A) The average community size s of the normal and tangential networks increases with strain step, but the tangential force communities tend to be smaller than the normal force communities.(B) The community strength σ measures the intra-community forces, and also grows as a function of strain step for both the normal and tangential networks.(C) The sparsity η quantifies the spatial density of community structure.The tangential force communities tend to be more dense than the normal force communities.In all plots, error bars correspond to the standard error of the mean over packings.
FIG. 11.Force drives local reconfiguration.(A, C) Scatter plots of particle flexibility ξ vs. the average force on a particle across compression f φ for a sample packing.In both the normal and tangential networks, there is a strong, positive Spearman correlation between the two quantities.(B, D) Scatter plots of particle flexibility ξ vs. the average absolute change in force on a particle across compression |∆f | φ for a sample packing.In both the normal and tangential networks, there is a strong, positive Spearman correlation between the two quantities.See Fig.21for the distribution of correlations for each packing).
FIG. 12.Community reorganization is driven by changes in intra-community force.Scatter plots show the community stationarity, ζc, vs. the average absolute change in community strength across compression, |∆σc| φ , for an example packing.For both the normal (A) and tangential (B) networks, there is a strong Spearman correlation between the two quantities (see Fig.22for the distribution of correlations for each packing).

7 FIG. 13 .
FIG.13.Changes in spatial arrangement across compression relate to community reorganization.Scatter plots show the average absolute change in community sparsity across compression, |∆ηc| φ , vs. community stationarity, ζc, for an example packing.For both the normal (A) and tangential (B) networks, there is a strong Spearman correlation between the two quantities.(See Fig.22for the distribution of correlations for each packing).

1 FIG. 14 .
FIG.14.Stationarity vs. absolute change in community strength for long-lived communities.Scatter plots show the community stationarity, ζc, vs. the average absolute change in community strength across compression, |∆σc| φ , for communities that exist for the entire compression process.For both the normal (A) and tangential (B) networks, there is an increased correlation between the two quantities across packings.

FIG. 15 .
FIG. 15.A comparison of the physical characteristics of the subsystem and bath communities.(A, D)Curves showing the average community size of the subsystem (red curve) vs. bath (black curve) communities for the normal and tangential networks.The small differences between the curves for each network ensures that we compare communities of similar size between the two groups.(B, E) The evolution of community strength in the network of normal forces is not significantly different between the subsystem and bath.However, the tangential network displays a clear separation, suggesting that the multilayer network model and community detection are sensitive to differences in friction.More importantly, the direction of the sensitivity agrees with physical intuition; the modular structures found in the low friction bath are characterized by lower tangential intracommunity force.(C, F) There is a small distinction between the community sparsity of the two groups, with the low friction communities being slightly more dense throughout compression.
Appendix A: Additional figures of community structure

FIG. 18 .
FIG. 18.Histogram of flexibility values over experimental packings for ω = 1 × 10 −200 .Even at a very small inter-layer coupling, the flexibility values of the granular networks do not increase by a significant amount compared to the values obtained at ω = 0.01.

FIG. 19 .
FIG.19.Network promiscuity values differ between the granular and null models.Significant differences between the reconfiguration of the granular network and the null models exist when the boxplots do not cross the zero line.(A) Boxplots showing that Ψ real < Ψ null for the normal and tangential networks, where the null model, A temporal , is constructed by scrambling network layers.(B) Boxplots showing that Ψ real < Ψ null for the normal and tangential networks, where the null model, A edges , is constructed by scrambling the edge weights within each layer while preserving contact structure.In all cases, statistical significance is determined by permutation testing.
FIG. 20.Inter-particle forces predict a more robust measure of local reconfiguration.The relationship between local network reconfiguration and the average force on a particle across network layers holds for a more robust measure of reorganization known as promiscuity.(A, C) Scatter plots of particle promiscuity ψ vs. the average force on a particle across all compressive steps f φ for a sample packing.For both the normal and tangential networks, there is a strong Spearman correlation between the two quantities.(B, D) Scatter plots of particle promiscuity ψ vs. the average absolute change in force on a particle across compression |∆f | φ , for a sample packing.For both the normal and tangential networks, there is a strong Spearman correlation between the two quantities.
The red line then denotes the median over packings and the edges correspond to the 25 th and 75 th percentiles.For the same inter-layer coupling, the normal and tangential networks have statistically different behavior in terms of these measures of network rearrangement.
FIG. 7. Measures of evolving community structure for the normal and tangential force networks.The structure and reconfiguration of the community structure of the normal and tangential force networks is characterized by network flexibility Ξ (A) and network stationarity ζ (B).The boxplots are constructed by averaging statistic values over optimization for each experimental packing.