Bob Behringer was one of the main players in the field of granular matter and contributed more than many to better understanding of these systems. When I first saw the beautiful images of force chains, produced by his experiments on photo-elastic particles [1–5], I remember thinking that this person had just tore a window in the wall of mystery surrounding the poorly understood granular systems. By devising a way to visualise force propagation paths in such media he made it possible to peer inside them and start to sort out the confusion surrounding them at the time. Indeed, his experiments gave rise to a new generation of predictive modelling.

And then I met the man and realised that there was much more to Bob than his experiments. Bob impressed me immensely well beyond his science. The combination of his levelheadedness and his love of life was inspiring. Very notable was his default ‘cool, calm and collected’ mode, with a pen between his teeth and a glint of humour in his eye. But that mode would last only until I described to him a potential new idea. The cool mode would disappear and he would react with the enthusiasm of a child. Indeed, his boundless curiosity made it easy to get him excited with good science. Beyond the science there were his abundant generosity and kindness. It was a pleasure to discuss with him anything, from fundamental science, through music and swimming to the meaning of life in general. Each such discussion felt like embarking on a voyage of discovery and almost always resulted in some more dots being connected and a better understanding emerging. Often, our discussions would meander into bantering, at which he was very good. It was very easy to forget how important his contribution had been to the development of the field and simply love him as a human being. Bob’s death left a large gap in the field, as well as a hole in the hearts of many. The field is not the same without him.

Granular materials (GMs) form an important and ubiquitous state of matter, playing a major role in our everyday life. They support and transmit stress significantly differently than conventional solids, often in the form of non-uniform stress chains [6–14]. Having a fundamental theory to predict the stress that develops in such media under given boundary loads is extremely significant for a wide range of applications in fields beyond science and engineering. All continuum static stress theories are based on balance conditions of force and torque:1∇·σ+Fext=0σ=σT,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&{\boldsymbol{\nabla}}\cdot \varvec{\sigma }+ {\mathbf{F}}_{\mathrm{ext}} = 0 \nonumber \\&\varvec{\sigma }= \varvec{\sigma }^T , \end{aligned}$$\end{document}where σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\sigma }$$\end{document} is the stress tensor, σT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\sigma }^T$$\end{document} its transpose, and F\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{F}}$$\end{document} includes all external and body forces on the medium. These equations must be supplemented by one constitutive relation in two dimensions (three in three dimensions). In elasticity theory, this is achieved by imposing compatibility conditions on the strain and relating the strain to the stress. This procedure results in elliptic equations whose solutions cannot be the observed stress chains. Alternative approaches, such as Critical State theory, popular in the engineering literature, rely on closures that relate the stress components when the medium is on the ‘yield surface’, which is a surface in the space spanned by the stress components. This approach gives a hyperbolic set of equations, yielding stress chains, or slip lines, but it is phenomenological in nature and lacks predictive power [15].

At marginal rigidity, GMs are minimally connected, possessing the lowest possible mean number of force-carrying contacts per particle, at which the structure is at mechanical equilibrium. In this marginally rigid, or isostatic, state, the intergranular forces in the assembly can be determined from balance conditions alone. This makes the strain, and hence stress–strain relations, redundant as input at the grain level. Since the continuum stress field is a coarse-grained representation of the spatial distribution of those forces then continuum strain-based constitutive properties are also redundant. For this reason, such theories, which includes elasticity, are inadequate for GMs.

Isostaticity stress theory has been developed to resolve the problem for ideally statically determinate, or isostatic, media [16–24]. Although most GMs are not precisely isostatic, on the verge of yielding they are sufficiently close to this state, containing isostatic regions, whose sizes depend on the proximity to marginal rigidity [21]. This makes isostaticity theory relevant to general GMs. Real systems can be generated close to marginal rigidity, e.g. by careful preparation or by generating low density shear bands [22, 23].

The irrelevance of strain for marginally rigid systems leaves structure as the only usable constitutive property for closing the stress equations and isostaticity theory closes the equations by one stress–structure relation in 2D (three in 3D). The most general form of the relation in 2D is [16–19]2∑ijqijσij≡Q:σ=0,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \sum _{ij} q_{ij}\sigma _{ij} \equiv {\mathbf{Q}}: \varvec{\sigma }= 0 , \end{aligned}$$\end{document}where Q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{Q}}$$\end{document} is a symmetric fabric tensor, quantifying the local structure in a specific way [25]. Isostaticity theory is the name given to the stress field equations (1) and (2). These equations are hyperbolic, resulting in solutions of stress chains extending along characteristic paths. While friction does not seem to appear explicitly in the equations, it affects the dynamics, that give rise to the static structure and, therefore, the statistics of the fabric tensor Q. Following empirical [16, 17] and mean field [19] proposals for the form of Q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{Q}}$$\end{document}, a first-principle theory derived it directly from the local grain-scale microstructure [20]. Specifically, around a grain g3Qg=12ϵ-1·Cg+CgT·ϵ,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{\mathbf{Q}}}_g = \frac{1}{2} {\varvec{\epsilon }}^{-1}\cdot \left[ {{\mathbf{C}}}_g + \left( {{\mathbf{C}}}_g\right) ^T\right] \cdot \varvec{\epsilon }, \end{aligned}$$\end{document}where ϵ=01-10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\epsilon }=\left( {\begin{matrix} 0 & 1 \\ -1 & 0\end{matrix}}\right) $$\end{document} is a π/2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi /2$$\end{document}-rotation matrix and4Cg=∑cCcg=∑crcg⊗Rcg,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{\mathbf{C}}}_g = \sum _c {{\mathbf{C}}}_{cg} = \sum _c {\mathbf{r}}_{cg} \otimes {{\mathbf{R}}}_{cg} , \end{aligned}$$\end{document}with the sum running over the loops c that surround grain g. The quadrilateral, whose diagonals are the vectors rcg\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{r}}_{cg}$$\end{document} and Rcg\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{R}}_{\mathrm{cg}}$$\end{document} (see Fig. 1), is called ‘quadron’, and it is the structure’s most basic volume element [27]. For brevity, I index the quadrons in the following cg by q rather than cg. This structural tensor, which has a 3D equivalent [28], is useful for most GMs except when the grains are extremely non-convex [29]. The quadrons surrounding grain g are said to belong to it and their number equals its number of contacts with neighbouring grains. This description allows us both to quantify the disordered local structure unambiguously and close the stress equations with a stress–structure coupling tensor Q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf{Q}}}$$\end{document}.

It can be readily verified that Tr{Qg}=∑q∈g\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{{{\mathbf{Q}}}_g\}=\sum _{q\in g}$$\end{document}Tr{Qq}=∑q∈grq·Rq\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{{{\mathbf{Q}}}_{q}\}=\sum _{q\in g} {\mathbf{r}}_{q}\cdot {{\mathbf{R}}}_{q}$$\end{document}, where q are the sum is over the quadrons belonging to grain g. Therefore, Tr{Qg}\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{{{\mathbf{Q}}}_g\}$$\end{document} quantifies the deviation of the quadrons around grain g from a kite form and, consequently, this is a measure of the net rotation of grain g relative to a global mean, which must be zero [20]. In other words, this quantity describes a local rotational fluctuation. It is constructive to consider5Pg=ϵQgϵ-1=∑qsqsq-tqtq/2,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\mathbf{P}}_g = {\varvec{\epsilon}}{{\mathbf{Q}}}_g {\varvec{\epsilon }}^{-1} = \sum _q \left( {{\mathbf{s}}}_{q}{{\mathbf{s}}}_{q} - {{\mathbf{t}}}_{q}{{\mathbf{t}}}_{q}\right) /2 , \end{aligned}$$\end{document}where sq\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf{s}}}_{q}$$\end{document} and tq\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf{t}}}_{q}$$\end{document} are shown in Fig. 1. The antisymmetric part of Cq\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf{C}}}_{q}$$\end{document} can be written as A{Cq}=Aqϵ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{A}}\{{{\mathbf{C}}}_{q}\}=A_{q}\varvec{\epsilon }$$\end{document}, with Aq=12|rq×Rq|\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_{q}=\frac{1}{2}|{\mathbf{r}}_{q}\times {{\mathbf{R}}}_{q}|$$\end{document} being the area of the quadron q. The area associated with grain g is Ag=∑qAq=A{Cg}\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_g=\sum _q A_{q} = {\mathcal{A}}\{{{\mathbf{C}}}_g\}$$\end{document} and the total area of the system is Asys=∑gAg\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_{sys}=\sum _g A_g$$\end{document}.Fig. 1

The loop sides of the quadron q, whose diagonals are rq\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{r}}_{q}$$\end{document} and Rq\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf{R}}}_{q}$$\end{document} and whose area is Aq\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_{q}$$\end{document} (shaded), are the vectors sq\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf{s}}}_{q}$$\end{document} and tq\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf{t}}}_{q}$$\end{document}. The symmetric tensor Pg\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{P}}_g$$\end{document} can be written in terms of these vectors, Pg=12∑qsq⊗sq-tq⊗tq\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{P}}_g=\frac{1}{2}\sum _q\left( {{\mathbf{s}}}_{q}\otimes {{\mathbf{s}}}_{q}-{{\mathbf{t}}}_{q}\otimes {{\mathbf{t}}}_{q}\right) $$\end{document}, with q being all the quadrons belonging to grain g. Note that tq=-sq′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf{t}}}_{q}=-{{\mathbf{s}}}_{q^{\prime }}$$\end{document} if q′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q^{\prime }$$\end{document} is adjacent to q. Therefore, the sum over Pq\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{P}}_{q}$$\end{document} inside a loop c vanishes identically irrespective of the loop shape

The first-principles derivation and, in particular, the ability to derive Qg\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf{Q}}}_g$$\end{document} from local structural characteristics [20] were an important development, but a new problem emerged. The equations of a complete stress theory should have the same form when upscaled, which requires that Qg\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf{Q}}}_g$$\end{document} can be coarse-grained systematically to arbitrarily large length-scales. This is problematic because the volume-average of Qg\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf{Q}}}_g$$\end{document} over any region of space vanishes except for contributions from the region’s boundary. This renders the closure relation (2) a useless trivial identity in the continuum limit. To see this, consider a region Γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document}, of boundary ∂Γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial \Gamma $$\end{document}, area AΓ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_{\Gamma }$$\end{document}, and NΓ∼L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{\Gamma }\sim L^2$$\end{document} grains. We have QΓ=⟨Qg⟩Γ=ϵ⟨Pg⟩Γϵ-1=PΓ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf{Q}}}_{\Gamma } = \langle {{\mathbf{Q}}}_g\rangle _{\Gamma } = {\varvec{\epsilon }}\langle {\mathbf{P}}_g\rangle _{\Gamma } {\varvec{\epsilon }}^{-1} = {\mathbf{P}}_{\Gamma }$$\end{document} and the volume-average of P\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{P}}$$\end{document} is6PΓ=12AΓ∑qsq⊗sq-tq⊗tq.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\mathbf{P}}_{\Gamma } = \frac{1}{2A_{\Gamma }} \sum _{q} \left( {{\mathbf{s}}}_{q}\otimes {{\mathbf{s}}}_{q} - {{\mathbf{t}}}_{q}\otimes {{\mathbf{t}}}_{q} \right) . \end{aligned}$$\end{document}Let us partition the sum over neighbouring pairs in contact, g and g′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g^{\prime }$$\end{document}, surrounding loop c. These terms cancel in pairs inside the loop because tq∈c=-sq′∈c\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf{t}}}_{q\in c}=-{{\mathbf{s}}}_{q^{\prime }\in c}$$\end{document}. It follows that the contribution of loops, fully enclosed within Γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document}, to the sum (6) vanishes—only that of the quadrons along the boundary ∂Γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial \Gamma $$\end{document} survives because those are not cancelled by neighbour quadrons. The vanishing of this sum is independent of the structural geometry, topology and grain shapes. It follows that PΓ∼1/L→0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{P}}_{\Gamma }\sim 1/L \rightarrow 0$$\end{document} as L increases, rendering (2) a trivial identity.

The method developed below consists of three steps. The first is by mapping the system into one with staggered order (SO), to be defined below. The second is coarse-graining the stress equations to the two-grain scale, with a fabric tensor that does not vanish identically to zero. The third step then proceeds as for any conventional coarse-graining method, by volume-averaging Q over increasingly large volumes until the desired continuum limit is reached.

1. Packs with staggered order (SO)

The resolution of this problem is based on the observation that Eq. (2), which was developed for frictional systems, can be rewritten in terms of only half the degrees of freedom in systems possessing a SO. SO in granular packs means that grains can be labeled + and −, such that any + grain is in contact only with − grains and vice versa. This is tantamount to the condition that each loop is surrounded by an even number of grains [20]. In such systems, the upscaling procedure is in the three steps detailed below: partitioning the assembly into pairs of + and − in contact; coarse-graining the constitutive equation to the two-grain scale by writing the stress equations for unit pairs in terms of half the degrees of freedom; upscaling the constitutive equation by conventional volume averaging.

The general force moment on grain g is Fg=∑hρgh×fgh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{F}}_g = \sum _{h}{\varvec{\rho }}_{gh}\times {\mathbf{f}}_{gh}$$\end{document}, where h are the neighbours of g, fgh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{f}}_{gh}$$\end{document} the force that grain h exerts on g, and ρgh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\rho }}_{gh}$$\end{document} are the position vectors of the contacts between g and h. The local stress on grain g is then σg=Fg/Ag\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\sigma }}_g = {{\mathcal{F}}}_g / A_g$$\end{document}. The mean and differential stresses of a +/− pair are, respectively, σm=(F++F-)/Apair\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\sigma }}_m = ({\mathcal{F}}_+ + {\mathcal{F}}_-)/A_{pair}$$\end{document} and σd=(F+-F-)/A±\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\sigma }}_d = ({\mathcal{F}}_+ - {\mathcal{F}}_-)/A_{\pm }$$\end{document}, where A±=A++A-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_{\pm } = A_+ + A_-$$\end{document} is the area associated with the pair. Averaged over a large region, Γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document}, ⟨σm⟩Γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {\varvec{\sigma }}_m\rangle _{\Gamma }$$\end{document} converges to the traditional continuum stress field, with the local stress fluctuations ⟨σd⟩Γ∼1/L→0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {\varvec{\sigma }}_d\rangle _{\Gamma } \sim 1/L \rightarrow 0$$\end{document}, with L the the system linear size. We combine the constitutive equations (3) for the pair, rearrange:7(Q++Q-):σm+(Q+-Q-):σd=0(Q+-Q-):σm+(Q++Q-):σd=0,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} ({{\mathbf{Q}}}_+ + {{\mathbf{Q}}}_-):{\varvec{\sigma }}_m + ({{\mathbf{Q}}}_+ - {{\mathbf{Q}}}_-):{\varvec{\sigma }}_d= & 0 \nonumber \\ ({{\mathbf{Q}}}_+ - {{\mathbf{Q}}}_-):{\varvec{\sigma }}_m + ({{\mathbf{Q}}}_+ + {{\mathbf{Q}}}_-):{\varvec{\sigma }}_d= & 0 , \end{aligned}$$\end{document}and average over Γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document}. The first equation becomes a trivial identity. The correlation between the two small quantities, ⟨σdQ++Q-⟩Γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {\varvec{\sigma }}_d\left( {{\mathbf{Q}}}_+ + {{\mathbf{Q}}}_-\right) \rangle _{\Gamma }$$\end{document} also decays with size, which leaves ⟨(Q+-Q-):σm⟩=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle ({{\mathbf{Q}}}_+ - {{\mathbf{Q}}}_-):\sigma _m\rangle = 0$$\end{document}, in which the subscript Γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document} was dropped for brevity. The averages of Q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{Q}}$$\end{document} and σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\sigma }$$\end{document} decouple as Γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document} increases:8⟨(Q+-Q-):σm⟩=⟨Q+-Q-⟩:⟨σm⟩.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \langle ({{\mathbf{Q}}}_+ - {{\mathbf{Q}}}_-):\sigma _m\rangle = \langle {{\mathbf{Q}}}_+ - {{\mathbf{Q}}}_-\rangle :\langle \sigma _m\rangle . \end{aligned}$$\end{document}Since ⟨Q++Q-⟩=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {{\mathbf{Q}}}_+ + {{\mathbf{Q}}}_-\rangle = 0$$\end{document} then ⟨Q+⟩=-⟨Q-⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {{\mathbf{Q}}}_+\rangle =- \langle {{\mathbf{Q}}}_-\rangle $$\end{document} and we can substitute this into (8) to obtain the upscaled relation9⟨Q+⟩:⟨σm⟩=0.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \langle {{\mathbf{Q}}}_{+} \rangle :\langle {\varvec{\sigma }}_m\rangle = 0 . \end{aligned}$$\end{document}This constitutive relation is a consistent coarse-grained version of the original one. Moreover, it has the significant advantage that the average of the constitutive quantity Q+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbf{Q}}}_+$$\end{document} does not vanish identically on coarse-graining like the original equation. Thus, this resolves the coarse-graining problem in SO systems.

The method presented next aims to coarse-grain the fabric tensor in packings of frictional grains. Packings of ideally frictionless grains have been shown in [20] to be mappable naturally to systems of frictional grains that possess SO. This is done by first assuming that the grains are frictional and then introducing infinitesimally small ball bearings at the contact points. It has been shown that the resultant system remains isostatic. Thus, for frictionless systems, Eq. (9) can be coarse-grained straightforwardly by volume averaging and the following procedure is not required.

2. Extension to general GMs

In general GMs, there are almost always loops surrounded by odd numbers of grains (OL), which ‘frustrates’ the SO, as at least two same-sign grains must be in contact.Fig. 2

The frustration around an OL surrounded by even-edged loops. Labelling the grains around it alternately by + and −, leaves two neighbouring ‘frustrated’ grains with the same sign (g0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_0$$\end{document} and gf\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_f$$\end{document}). Starting from a neighbour of g0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_0$$\end{document} and labelling the next shell (thick lines), leaves another frustrated pair lying adjacent to the first one. Continuing this procedure shell by shell around the original OL, results in a line of such frustrated grain pairs emanating from it. This line is capped when it is incident on another OL

In summary, the problem of upscaling the stress equations of isostatic GMs is unique in that the volume average of the constitutive fabric tensor, which couples to the stress to close the field equations, vanishes in the large-scale limit. This has led to imposing phenomenological or empirical fabric tensors. Particular examples of such closures are the yield conditions, such as Mohr–Coulomb, in plasticity-based theories. This problem is resolved here from first-principles by developing a specialised upscaling method. The method is based on the observation that the constitutive relation can be written in terms of only half the degrees of freedom in ideally unfrustrated granular packs, whose volume average over Q need not vanish identically in the continuum limit. Using this observation, the closure equation is first upscaled to the two-grain scale in such systems and from then on coarse-grained to the continuum conventionally by volume averaging.

Since most granular structures are frustrated, a ‘de-frustration’ method has been developed to transform any planar granular structure into an unfrustrated one. The method is based on joining frustrated grain pairs to lift the local frustration. The procedure leads to a renormalisation of the local fabric tensor and an example of such a calculation for a simple system including two defects in an otherwise honeycomb-like structure has been illustrated. It should be noted that the ‘de-frustrated’ structure need not be, and in most cases is not, marginally rigid. However, this does not pose a difficulty because the original physical system is marginally rigid and therefore isostaticity theory applies regardless of the mathematical manipulation of the structure.

The unusual aspects of this upscaling procedure are a direct consequence of the vanishing of a straightforward volume averaging of the constitutive quantity Q, a feature that is not common in any other coarse-graining procedure. This difficulty necessitates an upscaling in several stages: (1) de-frustrating the system into one with SO; (2) upscaling to the two-grain scale by writing the closure equation for pairs of grains and using half the degrees of freedom; (3) conventional volume averaging of the renormalised Q over increasing lengthscales.

It would be interesting to test the method on systems in which both the structure and the forces can be visualised, such as the many sample systems produced in the lab of Bob Behringer [1–5]. Moreover, the problem is still outstanding for three-dimensional systems, which have not been discussed here. Work on extension of the method to such systems is ongoing.