Stress Response of Granular Systems

TL;DR

This paper introduces a framework for understanding stress response in two-dimensional granular media, highlighting force localization and the influence of network structure and constraints on stability.

Contribution

It develops a graph Laplacian-based approach to analyze force response and localization in granular systems, incorporating microscopic force balance and gauge degrees of freedom.

Findings

Force localization is demonstrated through spectral analysis.

Numerical evidence from soft disk packings supports the theory.

Constraints like torque balance affect packing stability.

Abstract

We develop a framework for stress response in two dimensional granular media, with and without friction, that respects vector force balance at the microscopic level. We introduce local gauge degrees of freedom that determine the response of contact forces between constituent grains on a given, disordered, contact network, to external perturbations. By mapping this response to the spectral properties of the graph Laplacian corresponding to the underlying contact network, we show that this naturally leads to spatial localization of forces. We present numerical evidence for localization using exact diagonalization studies of network Laplacians of soft disk packings. Finally, we discuss the role of other constraints, such as torque balance, in determining the stability of a granular packing to external perturbations.

Equations

$$\sum _c{\vec{f}}_{g,c}=0,$$

$${\vec{f}}_{g,c}=-{\vec{f}}_{g^{\prime },c},$$

$${\vec{f}}_{g,c}={\vec{h}}_{v^{\prime }}-{\vec{h}}_v,$$

$$\hat{\sigma }$$

$${\hat{\sigma }}_g$$

$$\sum _c{\vec{r}}_{g,c}\times {\vec{f}}_{g,c}=0,$$

$$\sum _{(g,g^{\prime })}{\vec{r}}_{g,g^{\prime }}=0,$$

$$\sum _c{\vec{f}}_{g,c}=-{\vec{f}}_g^{\mathrm{body}}.$$

$${\vec{f}}_{g,c}={\vec{h}}_{v^{\prime }}-{\vec{h}}_v+{\vec{\varphi }}_{g^{\prime }}-{\vec{\varphi }}_g.$$

$${\square }^2\overrightarrow{{\varphi }_0}={\vec{\varphi }}_{g_1}+{\vec{\varphi }}_{g_2}+{\vec{\varphi }}_{g_3}+{\vec{\varphi }}_{g_4}-4{\vec{\varphi }}_{g_0}.$$

$${\square }^2|\vec{\varphi }⟩=-|{\vec{f}}^{\mathrm{body}}⟩,$$

$${\overrightarrow{\stackrel{~}{f}}}_{g,c}={\vec{f}}_{g,c}-\square \vec{\varphi }.$$

$${\square }^2=\sum _{i=1}^{N_G}{\lambda }_i|{\lambda }_i⟩⟨{\lambda }_i|,$$

$${\lambda }_1=0,\text{leavevmode}\text{leavevmode}|{\lambda }_1⟩=(\mathrm{111...1}).$$

$$\underset{({\square }^2{)}^{-1}}{\underbrace{\left(\sum _{i>1}\frac{1}{{\lambda }_i}|{\lambda }_i⟩⟨{\lambda }_i|\right)}}{\square }^2=\mathbb{I}-|{\lambda }_1⟩⟨{\lambda }_1|.$$

$$-({\square }^2{)}^{-1}|{\vec{f}}^{body}⟩=|\vec{\varphi }⟩-|{\lambda }_1⟩⟨{\lambda }_1|\vec{\varphi }⟩$$

$$=|\vec{\varphi }-\frac{1}{N}\sum _{i=1}^N\vec{\varphi }\text{leavevmode}\text{leavevmode}⟩.$$

$$|{\vec{\mathcal{F}}}^{body}⟩=|{\vec{f}}^{body}⟩-|{\lambda }_1⟩⟨{\lambda }_1|{\vec{f}}^{body}⟩$$

$$=|{\vec{f}}^{body}-\frac{1}{N}\sum _{i=1}^N{\vec{f}}^{body}⟩.$$

$$|{\vec{\mathcal{F}}}^{body}⟩=-{\square }^2|\vec{\varphi }⟩,$$

$$|\vec{\varphi }⟩=-({\square }^2{)}^{-1}|{\vec{\mathcal{F}}}^{body}⟩.$$

$$V({\vec{r}}_{g,g^{\prime }})=\frac{1}{2}{\left(1-\frac{|{\vec{r}}_{g,g^{\prime }}|}{{\sigma }_{g,g^{\prime }}}\right)}^2\Theta \left(1-\frac{|{\vec{r}}_{g,g^{\prime }}|}{{\sigma }_{g,g^{\prime }}}\right),$$

$$E_G=\frac{1}{N_G}\sum _{(g,g^{\prime })}V({\vec{r}}_{g,g^{\prime }}),$$

$$q^{-1}(\lambda )=\sum _j{e}_{j,\lambda }^4$$

Full text

††journal: Journal of Statistical Physics

*PACS numbers : 72.15.Rn ; 02.50.-r *

Contents

@starttoc

toc

1 Introduction

Force transmission in granular materials is a well-studied problem with wide ranging applications behringer_review_2016 ; majumdar_coppersmith_science_1995 ; majumdar_coppersmith_pre_1996 ; mueth_jaeger_pre_1998 ; cates_wittmer_prl_1998 ; bouchaud_claudin_crp_2002 ; geng_prl_2001 ; erpelding_epl_2010 and several properties of stress transmission in granular media continue to be the subject of active research. Several recent studies have addressed the question of the stress response of granular packings goldenberg_prl_2002 ; goldenberg_nature_2005 ; goldenberg_prl_2006 ; goldenberg_pre_2008 ; yan_arxiv ; reydellet_Clement_PRL2001 ; atman_epje_2005 ; atman_jphyscm_2005 ; atman_epl_2013 ; gland_epje_2006 , showing that the usual elasticity theories for homogeneous materials do not apply to materials with granular constituents, and therefore new frameworks need to be developed to deal with such systems. One of the most important characteristics that has emerged from these studies is the inhomogeneous nature of this stress propagation bouchaud_epje_2001 ; socolar_epje_2002 . A striking aspect of granular systems is that forces are primarily carried by a sparse, tenuous network of contacts that have become known as “force chains”. Experiments using photoelastic beads provide clear evidence of this phenomenon geng_prl_2001 , and force chains have emerged as the defining characteristic of granular solids. Yet, at present we do not have a theory of how forces are localized in space and the role played by network disorder on the character of this spatial localization geng_prl_2001 .

Static granular media is controlled by the constraints of mechanical equilibrium, which is the main ingredient in models of such systems. Several theoretical frameworks have been proposed to explain how granular materials respond to external forces bouchaud_review_2002 . These include lattice-based models such as the $q$-model majumdar_coppersmith_science_1995 ; majumdar_coppersmith_pre_1996 , and its extensions narayan_pre_2000 ; tighe_pg_2009 , and continuum models that posit some history-dependent relation between the components of the stress tensor bouchaud_cates_claudin_jpI_1995 ; vanel_howell_pre_1999 . These frameworks lead to Partial Differential Equations (PDEs) that are elliptic, hyperbolic or parabolic depending on these closure relations bouchaud_claudin_crp_2002 . The $q$-model, which incorporates scalar force balance in a model of granular piles in a gravitational field majumdar_coppersmith_science_1995 ; majumdar_coppersmith_pre_1996 , successfully accounts for the distribution of contact forces in the large force limit. In the continuum limit, the model reduces to the diffusion equation, predicting a horizontal spread of force-bearing contacts that grows as the square-root of the depth in the granular pile. Such a spread is observed in experiments under certain conditions of preparation geng_prl_2001 . Models with hyperbolic PDEs and disorder predict a wave-like propagation, which is similar to experimental observations in ordered arrays of grains geng_prl_2001 . Thus, the relationship between the stress response of granular materials and the underlying disorder of the contact network remains to be understood.

One of the important open questions in the field of granular systems is, how does one account for vector force balance on a disordered network of contacts? In this paper we address this problem by studying the response of a two dimensional granular material to an externally applied force. We show that the inhomogeneous propagation of stress at the grain level can be linked to the inherent randomness in the underlying fabric of contacts. We connect the problem of stress transmission in this system to that of diffusion on the disordered planar graph formed by the contacts between the constituent grains. This allows us to develop a theoretical framework to describe the response of such a granular system to an imposed body force, and thereby probe the origin of the force localization properties in such systems. In this work, we do not address the question of what the proper continuum level description is of stress propagation in granular media, rather, we demonstrate that the response at the granular scale is controlled by the spectral properties of the graph Laplacian describing the disordered contact network.

The constraints of mechanical equilibrium, that are necessarily satisfied in a static granular packing, lead to a gauge potential formulation of forces and stresses satake_mom_1993 ; ball_blumenfeld_prl_2002 . We show that the response of the internal forces to an applied body force can be described by an additional set of gauge potentials, which satisfy equations involving the Laplacian of the contact network. By introducing these auxiliary fields we can account for the change in the local stress tensor induced by an external perturbation. In addition to introducing the formalism, we present results of exact diagonalization studies of network Laplacians produced using numerical simulations of soft disk packings. These clearly illustrate the localization properties of force transmission in such materials. Although our formulation is valid for a general frictional granular packing, in our simulations we focus specifically on linear spring potentials of fricitionless soft disks. Finally we discuss the role of additional constraints such as torque balance in determining the stability of granular packings to external perturbations.

2 Local Constraints and Gauge Potentials

In this section we discuss the constraints that need to be satisfied at the local level in a static granular packing. The local nature of these constraints allows a description of the system in terms of gauge potentials. Local constraints are also crucial in determining the response of jammed packings to external perturbations, and give rise to deviations from linear elasticity.

2.1 Grains and Voids in Jammed Configurations

Granular systems are inherently porous with voids interspersed between grains in contact, and are inhomogeneous on the granular scale. This non-isotropic fabric of the underlying material leads to interesting properties of stress transmission between the grains in the system. A challenge in granular mechanics is to understand how this “graininess” kadanoff_review_rmp_1999 affects the bulk behavior of granular solids. In order to associate local quantities to granular packings with complex internal structures such as non-convex voids (see Fig. 1), it becomes necessary to associate well defined regions of space to specific parts of the packing. In two dimensional granular packings, the plane can be decomposed into polygonal regions belonging to grains as well as to voids (see Fig. 1). The grain polygons are formed by connecting the contact points on the boundaries of each grain in a cyclic manner. The void polygons are formed in a similar manner by cyclically connecting the contacts associated with each void. The grain and void polygons together tessellate the entire space. This construction allows a decomposition of the space into well defined polygonal regions and provides a way of probing the spatial structure of granular materials ramola_chakraborty_jstat_2016 ; ramola_chakraborty_arxiv_2016 . The two graphs formed by adjacent grains and adjacent voids are dual to each other. As we show in the next sections, this construction allows one to construct local gauge fields associated with these polygons which encodes the local force balance conditions.

2.2 Force Balance

The stability of granular materials stems from the fact that each individual configuration is in mechanical equilibrium. The internal stresses in jammed packings are mediated via contact forces between the constituent grains which are pairwise within the system. The mechanical equilibrium condition translates to the fact that the forces acting on every grain sum to zero. For packings with only contact forces within the system, this condition can be represented by the equation

[TABLE]

where $\vec{f}_{g,c}$ represents the force acting on the grain $g$, through the contact $c$, including both normal and tangential, frictional forces. The sum is taken over all the contacts $\{c\}$ for a given grain $g$. Next, Newton’s third law dictates that

[TABLE]

at each contact $c$ between the grains $g$ and $g^{\prime}$. These two equations (Eqs. (1) and (2)) can be used to construct an alternative representation of the forces in a two dimensional granular packing in terms of vector height fields which we discuss in the next section.

2.3 Height Fields

The force balance condition combined with Newton’s third law naturally leads to an alternative representation of the forces in the system, parametrized by a vector height field that lives on the edges, or the faces of the void polygons satake_mom_1993 ; ball_blumenfeld_prl_2002 ; henkes_chakraborty_pre_2009 ; degiuli_thesis ; degiuli_pre_2011 . In this work we choose to place these fields on the faces of the void polygons (see Fig. 2). As seen from Fig. 2, each contact in a two dimensional granular packing is shared by two void polygons. We define a set of height vectors on the voids $\{\vec{h}_{v}\}$, where $v$ indexes the voids in the system. The force at each contact $c$ is given by the difference of heights on the two voids associated with each contact

[TABLE]

with $v$ and $v^{\prime}$ ordered cyclically around the centre of the grain $g$. This naturally leads to Eq. (2), as the force on grain $g^{\prime}$ through the contact $c$ is simply $\vec{f}_{g^{\prime},c}=\vec{h}_{v}-\vec{h}_{v^{\prime}}$. It is interesting to note that although the contact forces are associated with the network of contacts between grains, the height vectors are associated with the dual lattice of the contact network i.e. the network of voids $\{v\}$. The cyclic convention allows one to visualize the heights as vector currents circulating in the counterclockwise direction on the edges of a void polygon ball_blumenfeld_prl_2002 . Given a set of contact forces, the definition of $\{\vec{h}_{v}\}$ is unique, modulo a choice of origin: they are gauge potentials for the stress tensor. This uniqueness of the height representation is ensured by the force balance condition (Eq. 1). Traversing a path around any non-trivial loop within the system returns to the same value of the height field. This can be easily seen by circulating around each grain, which form the basic loops in the system. As an example we consider a grain $g_{0}$ which has four neighbours $g_{1,2,3,4}$ in contact through contacts $c_{1,2,3,4}$ (see Fig. 2). Using the definition of the heights in Eq. (3), the forces acting on the grain $g_{0}$ are given by

[TABLE]

The left hand side represents the force on every contact of the grain $g_{0}$ which sums to zero, while the right hand side represents the difference in the heights at each contact. The sum is taken over all adjacent pairs of voids $(v,v^{\prime})$ surrounding the grain $g_{0}$. Starting with a value $\vec{h}_{v_{1}}$ on the first void, the values of the heights around the grain $g_{0}$ can be cyclically constructed using the contact forces. The force balance condition thus necessitates that the heights around the loop return to the same value $\vec{h}_{v_{1}}$. The non-trivial nature of the loops that can occur in higher dimensions is the main obstacle in extending this simple height construction beyond two dimensions.

Given the forces within the system, the stress tensor for a given packing is then defined as

[TABLE]

Here $\vec{r}_{g,c}=\vec{r}_{c}-\vec{r}_{g}$, with $\vec{r}_{c}$ being the position of the contact $c$, and $\vec{r}_{g}$ being the position of the centre of the grain $g$. $V$ represents the volume of the entire system. The vectors $\vec{e}_{g,v}$ define the vectors connecting the two contacts on grain $g$ (cyclically) that form an edge of the void polygon $v$ (illustrated in Fig. 1). Using the representation of $\hat{\sigma}_{g}$ in terms of the height vectors, it is easy to show that the total force moment tensor reduces to a boundary term henkes_chakraborty_pre_2009 , which is the discrete version of Stokes’ theorem. From a continuum perspective, the height fields can therefore be viewed as the gauge potential that enforces $\nabla\cdot\hat{\sigma}=0$ ball_blumenfeld_prl_2002 .

Finally it is important to note that a well defined height field at the local level can only be constructed in situations where pairwise forces are the only forces present in the system, since both Eqs. (1) and (2) are essential to the construction. Therefore several cases of interest are excluded from this framework, most notably systems where body forces, such as gravity, act on the grains. Similarly granular systems that are perturbed by an external force are also not amenable to a height description. It is therefore important to extend the gauge field treatment to such cases. In Section 3 we extend this height construction to such cases, which allows us to study the transmission of stress within two dimensional granular materials.

2.4 Additional Constraints

In addition to local force balance, a mechanically stable configuration must also satisfy torque balance at the grain level. These are represented by the set of equations

[TABLE]

where the sum is taken over all the contacts $\{c\}$ for a given grain $g$. In continuum, the torque balance constraint translates to the symmetric property of the stress tensor. Finally, in order for the packing to be valid, it must also satisfy additional constraints in real space. For example, the vector sum of inter particle distances taken over every closed trajectory in the system, must sum to zero. For the planar contact network in two dimensional packings, these constraints can be parameterized by the set of basic loops that enclose each void. These can then be represented by the set of equations

[TABLE]

where $\vec{r}_{g,g^{\prime}}=\vec{r}_{g^{\prime}}-\vec{r}_{g}$ is the interparticle distance vector between two adjacent grains located at $\vec{r}_{g^{\prime}}$ and $\vec{r}_{g}$ respectively. The sum is taken over all adjacent pairs $(g,g^{\prime})$ surrounding each void $v$. This illustrates the fact that for systems where the forces are related to the inter particle distances through a force law, the distrbution of stress in the system is intimately linked to the underlying real space network and its constraints.

3 Response to a Perturbation

In this section, we use the gauge potential formalism to analyse the response of granular packings to external perturbations such as an imposed body force while respecting the local force balance constraint. In the presence of body forces, the continuum equation of mechanical equilibrium is $\nabla\cdot\hat{\sigma}=-\vec{f}^{\rm body}$. We generalize the height field construction, which imposes mechanical equilibrium at the discrete grain level, to this situation. At the granular level, the mechanical equilibrium condition with an imposed body force $\vec{f}^{\rm body}_{g}$ on each grain $g$ now becomes

[TABLE]

with the sum taken over all the contacts $\{c\}$ of the grain $g$.

3.1 Auxiliary Fields

In order to account for the general case where contact forces and body forces are present within the system, we introduce an additional auxiliary vector field $\{\vec{\phi}_{g}\}$ associated with the grain polygons. This serves to generalize the height construction of Eq. (3). Since each contact has two associated grains and two associated voids, the contact forces are now given by the difference of height variables on the voids (as before) in addition to a difference of auxiliary fields on the grains as

[TABLE]

Here the contact $c$ is shared between the grains $g$ and $g^{\prime}$ and $v$ and $v^{\prime}$ are once again ordered cyclically around the centre of the grain $g$. It is easy to see that this definition satisfies Newton’s third law (Eq. 2). In the absence of body forces, the $\vec{\phi}$ field is identically zero and the above definition reduces to Eq. (3). As an illustrative example we consider the previous configuration (Fig. 2) with grain $g_{0}$ which is in contact with four neighbours $g_{1,2,3,4}$ through the contacts $c_{1,2,3,4}$, with additional body forces $\{\vec{f}^{\rm body}_{g}\}$ acting on each grain (see Fig. 3). The forces acting on the grain $g_{0}$ are given by

[TABLE]

The sum of the forces on the grain $g_{0}$ sum to the negative of the body force $\vec{f}^{\rm body}_{g_{0}}$ on the grain due to the mechanical equilibrium condition (Eq. (12)). The difference of height fields around the grain $g_{0}$ sum to zero as before. The summation on the right involving the $\vec{\phi}$ fields on the grains $g^{\prime}$ is simply the network Laplacian defined, on grain $g_{0}$ as

[TABLE]

It is straightforward to show that this Laplacian equation is valid for every grain. It is important to note that both the body forces and the auxiliary fields in the above equation are associated with the grains. In general, one can write the relationship between $\{\vec{f}^{\rm body}\}$ and $\{\vec{\phi}\}$ in a vectorial notation as

[TABLE]

where $|\vec{\phi}\rangle$ represents the vector $(\vec{\phi}_{g_{1}},\vec{\phi}_{g_{2}},...\vec{\phi}_{g_{N_{G}}})$ and $|\vec{f}^{\rm body}\rangle$ is the vector of body forces $(\vec{f}^{\rm body}_{g_{1}},\vec{f}^{\rm body}_{g_{2}},...\vec{f}^{\rm body}_{g_{N_{G}}})$. Here $N_{G}$ is the total number of grains in the system. The Laplacian operator $\square^{2}$ is now a matrix acting on these states. Eq. (20) represents our main result, and we can use it to analyse the response of a mechanically stable configuration to imposed body forces. Given a set of body forces and the contact network, we can invert this equation to obtain the auxiliary fields $\{\vec{\phi}\}$. The changes in the contact forces that develop as a response to an applied body force is then simply given by the difference of the $\vec{\phi}$ field at each contact (Eq. (18)).

Once we have determined the $\vec{\phi}$ field, we have local fields that incorporate the effect of the body forces at each contact. We can subtract the difference of these $\vec{\phi}$ fields, from the original contact forces to obtain “effective” contact forces ${\vec{\tilde{f}}}_{g,c}$ which satisfy the constraints of local force balance (Eq. (1)) and Newton’s third law (Eq. (2)) with

[TABLE]

Here $\square\vec{\phi}=\vec{\phi}_{g^{\prime}}-\vec{\phi}_{g}$ represents the gradient of the $\vec{\phi}$ field on the network. These effective forces can then be used to determine the height fields on the voids, using the construction described in the Section 2.3.

In the case where a body force is imposed on a granular packing already in mechanical equilibrium, the difference between the effective contact forces and the original ones represents the response of the granular packing to this perturbation. In the most general case, these new effective forces would induce changes in the real space network to satisfy the other constraints of mechanical equilibrium: torque balance and, in the case of frictional grains, the Coulomb condition of static friction. In many existing treatments of stress transmission and response in granular materials majumdar_coppersmith_pre_1996 ; bouchaud_review_2002 , structural changes are not allowed. The argument being that for rough, rigid grains, there is an indeterminacy at the contact level that allows for multiple force configurations to be consistent with a given real-space contact network. Experiments have also analyzed responses to external forces that occur without any changes in the network geng_prl_2001 . It should be remarked that allowing for network reorganization in response to external perturbations has been shown to have a significant effect on the coarse-grained description of stress transmissiongoldenberg_nature_2005 .

In this paper, we focus on the solutions for the effective contact forces and their localization properties in the absence of rearrangements of the contact network. For frictionless grains, which obey a given force law, a change in the forces in a given packing is necessarily accompanied by a structural change. For disks interacting via one-sided linear spring potentials, these real-space displacements can be obtained in closed form. In this case, the change in $\vec{\phi}$ fields directly represents the change in displacements of the particles, satisfying the required additional constraints of Section 2.4.

3.2 Inverting the Body Forces

The network Laplacian is the adjacency matrix of the graph representing the contact network with an added diagonal matrix whose entries are the number of contacts of the grain corresponding to that row, and has several well known properties. In our case the network is a disordered planar graph. $\square^{2}$ has the eigenfunction expansion

[TABLE]

with $N_{G}$ being the dimensionality, i.e. the number of grains in the system. $\square^{2}$ has one zero eigenvalue, with eigenvector

[TABLE]

This trivial zero mode is a consequence of the “conservation” law: $\sum_{j}\square^{2}_{i,j}=0$ for every row $i$ of the matrix representing the network Laplacian.The rest of the eigenvalues are all negative. We can next use the above eigenvalue expansion to invert the body forces and obtain the $\vec{\phi}$ fields in Eq. (20). We define a restricted inverse $(\square^{2})^{-1}$ of the Laplacian operator by projecting out the zero mode. We then have

[TABLE]

Next, using Eq. (24) in Eq. (20), we obtain the inversion

[TABLE]

3.3 Centre of Mass Frame

The inversion of the above equations is more natural when one considers the center of mass frame of reference. We can define

[TABLE]

Here $|\vec{\mathcal{F}}^{body}\rangle$ represents the vector of body forces in the center of mass frame. This then leads to a more symmetric formulation of Eq. (20)

[TABLE]

along with the inversion equation

[TABLE]

The above equations provide a unique solution to $|\vec{\phi}\rangle$ for a given network and a given set of body forces. If this solution fails to satisfy the other constraints of mechanical equilibrium such as torque balance, the network has to necessarily rearrange and indicates an instability of the network to this perturbation. Our current treatment, which focuses only on force balance, cannot address these questions of instability. In the figures illustrating the inhomogeneous response (Figs. 4 and 5), we have used body forces that are much smaller than the average force (and therefore the overlaps) between the grains, leading to very small changes in the contact forces, keeping the connectivity of the network unperturbed in the process.

As an illustrative example of the stress response within this framework, in Fig. 4 we plot the changes in contact forces that develop as a response to localized body forces in a jammed packing of soft frictionless disks. The body forces (represented by red arrows) act at the centers (i.e. centers of mass) of three grains separated by a small distance. This illustrates the effect of a localized stress perturbation to the packing. The reason for perturbing three grains is to create a non-trivial local perturbation that leaves the entire system in force balance. The changes in the contact forces in response to these body forces are obtained by solving Eq. (20) for a given initial jammed packing. The inhomogeneous nature of the stress response is clearly illustrated. As a more dramatic example, we plot this response for a sheared packing of soft frictional disks vinutha_unpublished in Fig. (5). The response provides characteristic signatures of the emergence of “force chains” in this case. To construct these responses, numerically simulated packings of frictionless and frictional grains were used to construct the network Laplacian. The above set of equations were then used to calculate $|\vec{\phi}\rangle$ for a network, which was then used to calculate the changes in the contact forces resulting from the imposed $|\vec{\mathcal{F}}^{body}\rangle.$ Below, we consider the relation of the response of an ensemble of packings to the spectral properties of an ensemble of graph Laplacians.

4 Spectrum of the Laplacian

We note that our formulation is closely related to the diffusion equation for stresses obtained by going to the continuum limit of the $q$-model describing stress transmission in a granular pile created under gravity bouchaud_review_2002 . In that formalism, the disorder is represented by the distribution of the $q_{i,j}$s that specify how the weight of a grain is split between different contacts. This disorder renormalizes the diffusion constant in the equation for the Greens function for stress propagation. Stresses transmit essentially vertically downwards with a small spread bouchaud_review_2002 . The network Laplacian formulation is able to take care of the underlying randomness in the system through the disorder in the contact network. This framework also demonstrates that the $q_{i,j}$s cannot be considered as independent, random stochastic variables. Instead, these variables that specify the redistribution of body forces are determined by the underlying network with its associated randomness.

4.1 Stress Localization

The connection to the problem of diffusion in the gauge potential formulation is immediately obvious as diffusion on the random planar graph is described precisely by Eq. (20). In addition, the discrete Laplacian matrix is important in many physical situations and several physical models on graphs rely heavily on it. For example, the kinetic energy term for hopping models, the dynamics of random vibrational networks, ferromagnetic $O(n)$ models as well as models of non-interacting bosons on graphs invoke the network Laplacian hastings_prl_2003 ; burioni_prl_2000 . Several intriguing connections to the problem of localization follow immediately. We can therefore use the already sophisticated machinery developed in the field of localization to study the phenomenon of force localization in granular systems.

A natural question to consider is then the contribution of the different eigenvectors of the Laplacian (Eq. (22)) to the stress response of the system. Fig. 4 provided an illustration of the total change in contact forces as a response to imposed body forces for a single configuration of jammed frictionless disks. In Fig. 6, we plot the response of the contact forces for the same configuration using just the largest negative and the smallest negative eigenvalues of the network Laplacian. We find clear signatures of localization in the higher end of the spectrum. To better understand these localization properties, we study packings of soft disks numerically.

4.2 Matrix Ensemble

It is clear from the formulation in Section 3, that the spatial localization, or lack thereof, of the $\vec{\phi}$ fields is determined by the nature of the eigenfunctions of the network Laplacian. The eigenvectors that possess a large overlap with $|\vec{\mathcal{F}}^{body}\rangle$ will contribute overwhelmingly to this response, and therefore it is important to understand how localized the eigenvectors of the Laplacian matrix are. In usual two dimensional localization problems, the disorder is manifest in the system as a quenched randomness, either explicitly in the interactions or in the spatial motion, which is then averaged over. In the case of granular packings, the disorder in the network of a given configuration plays the role of such a quenched variable. The ensemble of adjacency matrices can be thought of as an ensemble of random matrices with the randomness entering through the connectivities of the particles. The extent of stress localization is thus controlled by the ensemble of random matrices that represent the network Laplacians of disordered granular packings created through some protocol, and the nature of the perturbation.

5 Numerical Simulations of Frictionless Grains

In this section, we probe the spatial localization properties of the stress response in frictionless granular media using numerically simulated packings. We emphasize that our framework is not restricted to frictionless packings, we use frictionless packings as a model system to illustrate the application of our theoretical framework. To do this, we consider the paradigmatic example of a system of frictionless soft disks interacting via one sided linear spring potentials. To study the stress response of this system, we create a jammed packing of soft frictionless disks and perturb this with spatially localized body forces ellenbroek_pg_2005 . We then measure the statistics of the forces that develop as a consequence of the imposed body forces using the formulation developed in Section 3.

5.1 Jammed Packings of Soft Disks

We simulate a system of soft disks interacting via linear spring potentials of the form

[TABLE]

where $\Theta$ is the Heaviside function and $\sigma_{g,g^{\prime}}=\sigma_{g}+\sigma_{g^{\prime}}$ is the sum of the undistorted radii of disks $g$ and $g^{\prime}$. The total energy per grain of the system is given by

[TABLE]

where the sum is taken over all pairs $(g,g^{\prime})$, with $g\neq g^{\prime}$. $N_{G}$ is the total number of grains in the packing. We create jammed packings in mechanical equilibrium using a conjugate gradient minimization of Eq. (30). The number of grains that are part of the rigid structure of the contact network varies between different configurations, i.e. $N_{G}\equiv N_{G}-N_{R}$, where $N_{R}$ is the number of “rattlers”, particles that are not in contact with any of the others. This crucially decreases the dimensionality of the Laplacian matrix, making it singular, and therefore rattlers need to be removed from the system before any numerical procedure is implemented. We simulate systems of particles with a $50:50$ mixture of disks with diameter ratios $1:1.4$, at varying global energies between $E_{G}=10^{-15}$ and $10^{-5}$ ramola_chakraborty_jstat_2016 . The number of grains in our simulations vary between $N_{G}=512$ and $2048$.

5.2 Exact Diagonalization

We next exactly diagonalize the Laplacian matrix associated with the contact network of each configuration and measure the statistics of their eigenvalues and eigenvectors. We measure two characteristic signatures of localization in our system, namely, the density of states of the eigenvalues, and the inverse participation ratio slanina_epj_2012 ; clark_maestro_arxiv_2015 . The network Laplacian is a $N_{G}\times N_{G}$ real symmetric matrix with eigenvalues $\lambda_{i},i=1,...,N_{G}$ and corresponding normalized eigenvectors $|\lambda\rangle\equiv(e_{1,\lambda},e_{2,\lambda},...,e_{N_{G},\lambda})$. We measure the density of states $\rho(\lambda)$ of the eigenvalues of the Laplacian at different global energies. The density of states for an ensemble of $N_{G}=1024$ grains at different global energies is illustrated in Fig. 7. We find single isolated states within the lower spectrum of eigenvalues and a continuum of states with higher eigenvalues.

We next compute the inverse participation ratios of the eigenvalues for different energies. The Inverse Participation Ratio (IPR) corresponding to an eigenvector is defined as

[TABLE]

For a localized mode the IPR would be of $O(1)$ and for a delocalized mode this quantity would be of $O(1/N_{G})$. The IPR for an ensemble of $N_{G}=1024$ grains at different global energies is illustrated in Fig. 8. We find that a large part of the spectrum is in fact localized, with a small number of delocalized modes. This is similar to what one would expect in two dimensional disordered models abrahams_review_2010 , where states are in general localized. However, the nature of the disorder in granular systems still remains to be elucidated. This would require a detailed study of the properties of the random networks that arise in frictional granular packings, and is a fruitful direction for future research.

6 Force Tiles with Body Forces

In this section we discuss another application of the gauge potential framework developed in this paper, namely the construction of “force tiles” for systems with body forces. The condition of mechanical equilibrium associated with a static granular system with purely pairwise contact forces allows one to construct a useful representation of the forces in the system known as a force tile. This is constructed as follows. Using Eq. (1), the “vector sum” of the forces on each grain, i.e. the force vectors associated with the contacts of each grain arranged head to tail (with a cyclic convention), form a closed polygon. This forms a “force polygon” associated with this grain (see Fig. 9). Since the sum is taken cyclically over the contacts for each grain, we obtain convex polygons for frictionless systems. In frictional systems, the force polygons can be non convex and even self-intersecting making the graph non planar. Next, Eq. (2) imposes the condition that every force vector in the system, has an equal and opposite counterpart that belongs to its neighboring grain. This leads to the force polygons being exactly edge-matching, and one can then use this fact to tile these polygons next to each other (see Fig. 9). This construction produces a network known as the force tile network, with the edges representing the forces in the system. The two representations: force tilings and height fields defined in Section 2.3 are related. It is easy to see that the positions of the vertices of the force tiles represent the height vectors starting from an arbitrary origin, since the vector distance between these vertices provides the forces (in correspondence with the definition of the heights). Since each face of the force tiling graph is uniquely associated with a grain, the adjacency of the faces is the adjacency of the grains in a packing. The adjacency of the vertices are simply the adjacency of the voids (since the heights are associated with the void polygons). The force tiling representation has several intriguing properties tighe_jstat_2011 , an important one being that the distances between the vertices represents a measure of the amount of stress between two points in a packing. The extent of the tiling provides the total amount of stress in the packing. Since force tilings provide information about the nature of the stress distribution within a system, they provide sensitive measures with which to characterize stress induced transitions in granular media sarkar_prl_2013 ; sarkar_pre_2016 . In the presence of body forces, however, the force tile construction fails tighe_pg_2009 . This is because the force polygons no longer close as the contact forces do not sum precisely to zero. This makes constructing force tiles for granular piles and for systems with hydrodynamic drag impossible.

The Laplacian framework developed in this paper can be used to extend the construction of force tiles for systems where contact forces are not the only forces in the system. Given the contact and body forces in the system, by using the network of contacts one can construct the $\vec{\phi}$ field as detailed in Section 3. This then allows for a unique construction of the height field $\{\vec{h}\}$ which are the vertices of the force tiling. As an illustrative example, in Fig. 10 we show the force tiling associated with a dense suspension of $2000$ soft disks with pairwise hydrodynamic interactions and high drag forces mari_unpublished . In this case, the viscous drag acts as a body force on each grain. The positions of the vertices were computed by inverting Eq. (20).

7 Discussion

In this paper we have discussed how stresses in granular packings are transmitted through the network of contact forces in response to external perturbations. Our formalism allowed us to construct force-balanced solutions on disordered networks that respect vector force balance at the microscopic level. This new potential formulation opens several interesting avenues. The Laplacian framework can be easily extended to systems where contact forces are not the only forces present in the system. This occurs frequently in systems with hydrodynamic forces, where viscous drag plays a major role, and granular piles in a gravitational field.

Our construction assumes that there is enough indeterminacy in the forces at the contact level depending on how a packing is created. This is particularly true for systems with frictional grains where an exact force law is not applicable. Even so, the constraints that are left out of the analysis such as the Coulomb constraint and torque balance would change the effect of the response. If we assume, as in the $q$-model, that roughness at the grain level would lead to an indeterminacy in the actual position of the contacts, then one can find a torque balanced solution as long as the perturbation is small enough, and the network does not rearrange. Our construction therefore accurately describes systems near the infinitely rigid limit for which there is a large separation scales between forces and displacements vanhecke_prl_2004 . Generalizations of the framework presented in this work to account for the other constraints of mechanical equilibrium, and allowing for network rearrangements would be an interesting avenue for future research.

8 Acknowledgements

This work has been supported by NSF-DMR 1409093 and the W. M. Keck Foundation.