The unusual problem of upscaling isostaticity theory for granular matter

Isostaticity theory (IT) provides stress field equations for marginally rigid packs of non-cohesive particles, interacting through hard-core forces. Its main advantage over strain-based theories is by closing the stress equations with stress-structure, rather than stress-strain, relations, which enables modelling the stress chains, often observed in experiments and simulations. The usefulness of IT has been argued to extend beyond its applicability at marginal rigidity. It has been shown to be derivable from first principles in two-dimensions, with the structure quantified by a fabric tensor that couples to the stress field. However, upscaling IT to the continuum is done currently empirically by postulating convenient closure equations. The problem is that a volume average of the fabric tensor vanishes in the continuum limit, trivialising the closure equation. This poses an unusual upscaling problem, necessitating a new approach. Such an approach is developed here, resolving the problem for planar granular assemblies. The new method is developed initially for idealised 'unfrustrated' packs by coarse-graining first to the two-grain scale, after which a conventional coarse-graining can be used. It is then extended to general realistic systems, by introducing an intermediate `de-frustration' procedure. The applicability of the method is illustrated with a tractable example.


About Bob:
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], 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.

Introduction:
Granular materials (GMs) form an important and ubiq-uitous 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 [2]- [10]. 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: where σ is the stress tensor, σ T its transpose, and F 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 [11]. At marginal rigidity, GMs are minimally connected, possessing the lowest possible mean number of forcecarrying 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 [12]- [20]. 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 [17]. 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 [18,19].
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 [12]- [15] ij where Q is a symmetric fabric tensor, quantifying the local structure in a specific way [21]. 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 [12,13] and mean field [15] proposals for the form of Q, a first-principle theory derived it directly from the local grain-scale microstructure [16]. Specifically, around a grain g where ǫ = 0 1 −1 0 is a π/2-rotation matrix and with the sum running over the loops c that surround grain g. The quadrilateral, whose diagonals are the vectors r cg and R cg (see Fig. 1), is called 'quadron', and it is the structure's most basic volume element [23]. For brevity, I index the quadrons in the following cg by q rather than cg. This structural tensor, which has a 3D equivalent [24], is useful for most GMs except when the grains are extremely non-convex [25]. 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.
It can be readily verified that Tr{Q g } = q∈g Tr{Q q } = q∈g r q · R q , where q are the sum is over the quadrons belonging to grain g. Therefore, Tr{Q g } quantifies the deviation of the quadrons around The loop sides of the quadron q, whose diagonals are rq and Rq and whose area is Aq (shaded), are the vectors sq and tq. The symmetric tensor Pg can be written in terms of these vectors, Pg = 1 2 q (sq ⊗ sq − tq ⊗ tq), with q being all the quadrons belonging to grain g. Note that tq = −s q ′ if q ′ is adjacent to q. Therefore, the sum over Pq inside a loop c vanishes identically irrespective of the loop shape.
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 [16]. In other words, this quantity describes a local rotational fluctuation. It is constructive to consider where s q and t q are shown in Fig. 1. The antisymmetric part of C q can be written as A{C q } = A q ǫ, with A q = 1 2 |r q × R q | being the area of the quadron q. The area associated with grain g is A g = q A q = A{C g } and the total area of the system is A sys = g A g .
The problem: The first-principles derivation and, in particular, the ability to derive Q g from local structural characteristics [16] 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 Q g can be coarse-grained systematically to arbitrarily large length-scales. This is problematic because the volume-average of Q g 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 Γ, of boundary ∂Γ, area A Γ , and N Γ ∼ L 2 grains. We have Q Γ = Q g Γ = ǫ P g Γ ǫ −1 = P Γ and the volume-average of P is Let us partition the sum over neighbouring pairs in contact, g and g ′ , surrounding loop c. These terms cancel in pairs inside the loop because t q∈c = −s q ′ ∈c . It follows that the contribution of loops, fully enclosed within Γ, to the sum (6) vanishes -only that of the quadrons along the boundary ∂Γ 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 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 volumeaveraging Q over increasingly large volumes until the desired continuum limit is reached.
The upscaling method: I. 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 [16]. 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 F g = h ρ gh × f gh , where h are the neighbours of g, f gh the force that grain h exerts on g, and ρ gh are the position vectors of the contacts between g and h. The local stress on grain g is then σ g = F g /A g . The mean and differential stresses of a +/-pair are, respectively, σ m = (F + + F − )/A pair and σ d = (F + − F − )/A ± , where A ± = A + + A − is the area associated with the pair. Averaged over a large region, Γ, σ m Γ converges to the traditional continuum stress field, with the local stress fluctuations σ d Γ ∼ 1/L → 0, with L the the system linear size. We combine the constitutive equations (3) for the pair, rearrange: and average over Γ. The first equation becomes a trivial identity. The correlation between the two small quantities, σ d (Q + + Q − ) Γ also decays with size, which leaves (Q + − Q − ) : σ m = 0, in which the subscript Γ was dropped for brevity. The averages of Q and σ decouple as Γ increases: Since Q + + Q − = 0 then Q + = − Q − and we can substitute this into (8) to obtain the upscaled relation 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 + does not vanish identically on coarse-graining like the original equation. Thus, this resolves the coarsegraining 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 [16] 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 coarsegrained straightforwardly by volume averaging and the following procedure is not required.

II. 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.
The extension of the above method to general systems is by renormalising the structure to lift the frustration. It is convenient to focus first on systems without 3-edge loops -the extension to include these will be described below. Consider a single OL, within a region of only even-edged loops, such as sketched in Fig. 2. Label the grains around it as + and -in the clockwise direction, alternatingly. The last grain, g f , is in contact with the first one, g 0 , and they are both −. Starting from a neighbour of g 0 , which does not belong to the OL, label similarly the first shell of loops surrounding the OL. With this shell's loops being even-edged, it must contain exactly one frustrated pair of grains and these, which is adjacent to the pair g 0 -g f . Repeating this process shell . The area of the super-grain is the sum of the areas of its constituents, which is also the coefficient of the antisymmetric part of the renormalised geometric tensor C pair . Since the renormalised fabric tensor P is the sum over the tensors of the frustrated grains, their common vectors s and t between the joined pair (thick black lines in (a)) are also eliminated in (b), disposing of the effect of the contact point in the structural description. The stress field is unaffected by this because these vectors cancelled out anyway in the original structure. Therefore, the joining procedure is self-consistent.
by shell outwards, a continuous line of same-sign pairs emanates from the single OL, as sketched in Fig. 2. This line extends to the system boundary. However, if the line is 'incident' on another OL, it ends, as sketched in Fig.  2. This observation is the key to renormalising the structure for lifting the frustration. Isolating same-sign grain pairs into lines and regarding each such pair as a rigid super-grain, we recover a SO structure! The procedure is the following. Firstly, identify all the OLs in the network. Next, partition the OLs into nearest pairs, minimising the overall number of contacts between them. Thirdly, identify a 'frustration' line between each OL pair, using the above method and avoiding crossing of these lines. Finally, declare each frustrated pair of grains a supergrain and compute its fabric tensor, Q as the sum of the fabric tensors of its constituents. The result is a system of rigid objects, some the original grains and some the super-grains, possessing a SO. Finally, use the above coarse-graining procedure to upscale the constitutive relation to regions of required size, namely, containing sufficient numbers of loops each.
Although this renormalisation changes the local topology, it does not alter the forces on the original grains. Turning two grains into a super-grain eliminates the contact between them and consequently some quadrons, which are the descriptors and basic elements of the local structure. Since the structure tensor of the super-grain is the sum of those of its constituents, C ± = C + + C − , then the area associated with the super-grain is A ± ǫ = A{C ± } = (A{C + } + A{C − }) = (A + + A − ) ǫ, as expected. Another feature of the procedure concerns the physics of the stress field. Considering eq. (4), eq. (6) and Fig. 3, note that the vectors s and t, which originally extended between the eliminated contact point and the centres of the loops flanking it, cancel out and they do not play any role in the renormalised tensor Q ± that results from the pair of the eliminated grains. Since the stress-structure coupling can be mediated only by grain contacts, the removal of the effect of the eliminated contact point by the local renormalisation is significant -it underpins the self-consistency of the procedure with its local structural and geometrical interpretation.
3-edge loops are special because joining any two of their grains annihilates the loop altogether. Such loops are very sparse in isostatic frictional systems, whose mean number of loop edges is 6 (the result of Euler's theorem). These loops are pre-processed before the upscaling to the two-grain scale, by merging the loop's three grains into an initial super-grain. In the rare occasion that clusters of such 3-edge loops occur the entire cluster is merged into a super-grain. The de-frustrastion procedure is then applied to the resulting structure, which is devoid of 3-edge loops.
Once the 3-edge loops have been eliminated, the procedure described in the paper works for any concentration of odd loops and, in particular, when the entire system is made of them. Such systems are highly unlikely to occur but, if they do, the upscaling procedure would result in the two grains, shared by every pair of neighbouring loops, being joined into a super-grain. This loses an edge for each cell, making them all even-edged, and reduces the number of loops to a half.

Example:
To illustrate the method, consider first the granular pack shown in Fig. 4. The underlying structure is ordered on a honeycomb-like lattice, possessing SO. A local defect has been introduced by expanding one loop into an octagon, generating two pentagonal neighbours. The latter are OLs, introducing two sources of local frustration. The structure tensor C is purely antisymmetric in the regular hexagonal regions and Q g vanishes. The frustration defects, highlighted by ellipses in the figure, introduce a local finite fluctuation in the fabric tensor and we wish to calculate the renormalised effect on the constitutive relation after renormalising to SO. Carrying out the above joining procedure on the frustrated grains, and summing over the resulting Q tensors of the positive grains, we obtain the renormalised contribution: Substituting into (9), the local constitutive relation around the defect is then Generally, adding a uniformly random distribution of defects at area density ρ to a defect-free isostatic system, of an initial fabric tensor Q 0 , the coarse-grained fabric tensor of the effective medium is Q em = (1 − ρ)Q 0 + ρQ def ect . Using this constitutive relation together with (11) provides the stress equations at any desired length scale. It is useful to illustrate the issue with a system that is partly uniform and partly with defects. Suppose the uniform medium is the above honeycomb, deformed slightly to perturb the lattice's symmetry, which trivialises the diagonal of Q 0 . Placing particles in contact on the vertices makes a Kagome structure. The deformation can be chosen to yield a range of fabric tensors of the form Q 0 = Ka 2 1−ξ −(1+ξ) −(1+ξ) 1−ξ , with K a constant and ξ < 1 [26]. For simplicity, I choose a deformation that yields ξ = 1/4. In the uniform medium, this gives rise to characteristics oriented at gradients λ 0 1 = 1/3 and λ 0 2 = 3 (A and B in Fig. 5). Let this medium fill the entire half-plane −∞ ≤ y ≤ ∞ and 0 ≤ x ≤ ∞. Now fill the the region x 0 ≤ x ≤ ∞ with the aforementioned defects at density ρ = 10%, the fabric tensor in this region is The resulting characteristics in this region are at gradients λ em 1 ≈ 0.17 and λ em 2 ≈ 26.70. Applying at the origin the stress σ 0 = sxx 0 0 0 gives rise to stresses along the left two characteristic paths, A and B: σ A = −3sxx 1/3 1/9 . These chains reach the boundary of the defect-filled region at points (x 0 , x 0 /3) and (x 0 , 3x 0 ), whereupon they act as boundary loads for the medium in x > x 0 . From each of these load sources emanates a pair of stress chains at gradients λ em 1 and λ em 2 . The stresses along these paths is calculated using the method outlined in [20]. tions 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 'defrustration' 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 'defrustrated' 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: (i) de-frustrating the system into one with SO; (ii) upscaling to the two-grain scale by writing the closure equation for pairs of grains and using half the degrees of freedom; (iii) 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]. 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.
Compliance with ethical standards: I declare that I have no conflict of interest.