Stresses in two-dimensional isostatic granular systems: exact solutions

It has been recognized that the concept of isostaticity holds the key to understanding stress transmission in cohesionless granular media. Here the field equations of isostaticity theory in two-dimensions are studied and solutions are derived. The equations are first decoupled into integro-differential equations for the three independent stress components, highlighting the role of a particular position-dependent fabric tensor. In disordered, but statistically isotropic, systems the fluctuations decay with length-scale and the decoupled equations can be expanded to first order and solved. The zero-order solutions are obtained in closed form and give rise to force chains that propagate along straight characteristic lines. At this order solutions do not attenuate and chains cross one another without interference or scattering. The analysis of the first-order correction reveals emergence of weaker secondary force chains that branch off the zero-order chains and into a ‘cone of influence’ that they span.


Introduction
The physics of granular matter has attracted much attention in recent years due to their ubiquity in nature and the dominant role that they play in everyday life. Examples are soil, gravel, beans, agricultural grain and seeds, pharmaceutical products, and powders, to name a few [1]. Most technological materials pass through a particulate form at some stage during their processing and the handling of particulates is an essential issue to a wide range of industries. From a theoretical point of view, granular matter exhibits a wide range of complex behaviour that is far from fully understood, exhibiting features that are solid-like, liquid-like, gas-like, or a simultaneous combination of these [2].
A granular system is an assembly of macroscopic particles whose large size has two main consequences: one is that it makes thermal fluctuations irrelevant and another is that particles interact practically upon direct contact. Dry granular materials transmit only compressive forces while in wet or charged materials attractive forces are also possible. The following discussion focuses only on compressive interactions.
To understand the mechanical behaviour of granular matter it is essential to make the link between transmission of forces on the granular level and a macroscopic description of a continuous stress field. It is straightforward to write down the equations that govern the discrete intergranular forces, namely, Newton's equations of balance of forces and torque moments. It is the translation of those into a set of differential equations for a representative stress field that has proved elusive. Observations of non-uniform stresses in systems of dry grains [3,4] and in numerical simulations [5] prompted a search for a description that goes beyond conventional elasticity theory. The reason is that the latter cannot explain straightforwardly such stresses without resorting to elaborate mechanisms for long range anisotropic organization of bulk elastic constants. This led to several models of stress transmission, using empirical [6] and statistical [7]- [9] approaches.
It has been recognized that to make progress in the field it is essential to understand first stress transmission in isostatic, or statically determinate (SD) systems [7,8]. The latter are structures whose intergranular forces can be determined from the balance of forces and torque moments alone. To be SD, a macroscopic pack of infinitely rigid grains has to satisfy only one condition -that the mean number of force-carrying contacts per grain has a particular value z c . This value depends on the dimensionality of the system d, on the roughness of the particles and on whether or not they are spherical. It also depends on the ratio of boundary to volume particles, but this dependence tends to zero rapidly with system size and can be neglected for large systems. For packs of rough particles z c = d + 1, for frictionless and arbitrarily-shaped particles z c = d(d + 1), and for smooth spheres z c = 2d [10].
Granular assemblies in nature are not necessarily isostatic and the usefulness of this concept to general particulate systems has been questioned [11]. In particular, two aspects of the isostatic picture came under fire. One is that real particles always have a finite rigidity, however high, and the other is that in most real materials the mean number of contacts per particle in mechanical equilibrium is larger than z c . Both these issues have been addressed recently and it has been concluded that they do not invalidate the relevance of isostaticity to more general granular matter. A detailed discussion of this issue is outside the scope of this paper, but it deserves some elaboration. Firstly, it has been established in [12] that packs of compliant grains of mean contact number z c do transmit stresses as ideal isostatic systems, albeit with corrections that decay rapidly with system size. This issue has also been studied in [13]. Secondly, Blumenfeld [14] introduced the idea that when the mean contact number is larger than z c then elastic regions form and most realistic granular materials are in fact two-phase composites-part isostatic and part elastic. Therefore, understanding the mechanics of pure isostatic systems is an essential step to the modelling of granular systems in general. The purpose of this paper is to develop further a recently-proposed isostaticity theory for two-dimensional (2D) systems [14,15].
The most significant consequence of static determinacy to the present discussion is that it vitiates stress-strain constitutive relations. This can be understood as follows. Being SD, the intergranular forces can be determined irrespective of grain compliance. Since the stress is only a continuous representation of these forces then the elastic properties of the material are redundant to the determination of the stress field. This means that approaches that use strainbased information, including elasticity theory, are ineffective for these materials. This realization has given rise to an extensive search for a new set of stress field equations.
Any continuous static stress theory must include balance conditions. In 2D these consist of two force equations, and one equation of global torque moment, In these equations g i are the component of an external force field g = (g x , g y ), which may include body forces. The closure of these equations requires a relation between the stress field and some constitutive properties. Elasticity theory imposes compatibility conditions on the strain and supplements it by phenomenological or empirical stress-strain relations [16]. In contrast, isostaticity theory closes the equations by the relation p xx σ yy + p yy σ xx = 2p xy σ xy .
Originally, this equation had been suggested on empirical grounds [6], with p ij values designed to fit observations of arches in conical granular piles. More recent work has established this relation from first principles [15]. There the parameters p ij have been shown to be the components of a symmetric rank-two fabric tensorP that characterizes the local structure. Note that equations (1)- (3) contain no information about the 'constitution' of the material and such information is only injected in equation (4). Therefore, the latter can be regarded as a constitutive local stress-structure relation. It is interesting to note that equation (4) represents a rare case where a constitutive relation can be derived rather than postulated or measured.
The values of p ij have been shown to fluctuate on the granular scale around zero mean, regardless of the topology or anisotropy of the structure, and this makes the coarse-graining of relation (4) far from straightforward [15]. This difficulty has been resolved eventually by the construction of a specialized procedure that takes advantage of inherent local anti-correlations between neighbouring grains [17]. A review of the coarse-graining procedure is outside the thrust of this discussion, but one of its consequences is relevant to the following discussion. On the granular level, the procedure averages, almost everywhere, over every second grain. Combined with the anti-correlations, this gives rise to finite mean values of the renormalized p ij and it is this that makes possible the upscaling of equation (4). In its coarse-grained version, this equation is valid for macroscopic scales with p ij a continuous fabric tensor. In a later development [14] it was shown that, under a simplifying assumption, to be pointed out below, the field equations can be solved in closed form. The resulting solutions were shown to consist of stress chains that cross the system. This paper reports an analysis of the field equations beyond the discussion in [14]. The following are the main results.
1. The equations are decoupled into individual integro-differential equations for the stress components σ ij . 2. The stress is expanded to first-order in terms of the deviations of the fabric parameters p ij from their mean values and an equation is derived for the correction field. This shows that the solution obtained in [14] is the zero-order term in the stress expansion. 3. It is found that the zero-order stress chains do not attenuate away from the load source and that, on incidence, they pass through one another unchanged. 4. It is found that the correction field generates secondary side-branches into a cone of influence. The validity and usefulness of the approach is discussed and supported numerically.

The decoupled equations
The coarse-graining procedure of [17] gives a renormalized fabric tensor whose trace, p xx + p yy , is predominantly of a specific sign. Whether positive or negative is immaterial both to the physics of the problem and, as can be seen from equation (4), to the solution of the equations. Let us then consider initially systems whose p ij = 0 everywhere. The special cases, when one of the p ij vanishes locally, are simpler to analyse and are detailed in appendix A. Substituting σ yy from equation (4) into (2) gives where q ij ≡ p ij /p xx . Differentiating equation (1) with respect to y, multiplying the result by q yy and adding to it the derivative of equation (5) with respect to x, gives q yy ∂ yy + ∂ xx + 2q xy ∂ xy σ xy − ∂ xy q yy + ∂ y q yy ∂ x + ∂ x q yy ∂ y σ xx + 2 ∂ xy q xy + ∂ y q xy ∂ x + ∂ x q xy ∂ y σ xy = ∂ x g y + q yy ∂ y g x .
Integrating equation (1) gives where φ(y) is an arbitrary function of y. Substituting now this expression for σ xx in (6), gives an integro-differential equation for σ xy . Upon multiplication by p xx and rearrangement, this equation reads (p xx ∂ xx + 2p xy ∂ xy + p yy ∂ yy )σ xy + 2p xx (∂ y q xy ∂ x + ∂ x q xy ∂ y )σ xy + p xx ∂ y q yy ∂ x + ∂ x q yy ∂ y + ∂ xy q yy x ∂ y σ xy dx = p xx ∂ x g y + p yy ∂ y g x + p xx ∂ y q yy ∂ x + ∂ x q yy ∂ y + ∂ xy q yy φ(y) + x g x dx . (8) The right-hand side of equation (8) does not involve the stress and can be regarded as a 'source' term. The function φ(y) is arbitrary and, for simplicity, we shall put it to zero in the following. Thus, the entire right-hand side is a known function of position. The first operator on the left-hand side of (8) isP :H ≡ L, where H ij = ∂ ij is known as the Hessian. The equations for σ xx and σ yy can be obtained in a similar fashion. It has been shown in [14] that, under the assumption that the p ij are constant, the equations for all the σ ij are identical but for the source terms. This result is reproduced and extended in appendix B.

Analytical solutions
For clarity, the following discussion focuses on σ xy but the analysis and results are also valid, up to quantitative details, for the diagonal terms, which follow similar equations. At first glance, equation (8) may appear too formidable to surrender a closed-form solution. Certainly, it can be solved numerically, for example by an iterative procedure. Starting from an initial guess for σ xy , it can be substituted into the integral on the left-hand side and the integral computed. Regarding this as an additional source term, the resulting differential equation can be solved for σ xy . Treating the solution as the next guess, it is inserted back into the integral on the left-hand side and the procedure is repeated. The hope is that after a reasonable number of iterations the solution and the guess converge to an acceptable accuracy. The disadvantages of this procedure are that it is time-consuming, sensitively dependent on the proximity of the initial guess to the actual solution, and that there is no guarantee that it converges quickly or even at all. Consequently, it is questionable whether this approach is easier than a direct brute-force solution of the original coupled field equations [18]. Another approach is to treat directly the coupled balance equations as a hyperbolic system and this is currently under study [19]. The approach taken here is based on the expectation that, barring long-range correlations, coarse-graining in isotropic systems should smoothen out the grain-scale fluctuations of the fabric parameters p ij progressively with size. Consequently, above some length-scale, which is still well below the size of the system, the fabric tensor can be represented by p ij = π ij (1 + δ ij ).
Here π ij are constant and δ ij 1 are space-dependent fluctuations of zero mean. Moreover, the gradients of p ij will also be ignored. The validity of these approximations will be discussed in the concluding section. Here and in the following the term π ij δ ij implies no summation over i and j. Using this form, let us expand the stress tensor to first order in δ ij , σ xy = σ (0) xy + σ (1) xy , such that σ (1) xy is proportional to a linear combination of the δ ij s. The idea is to construct the equation that governs σ (1) xy and analyse it.

The zero-order equation and its solution
The zero-order equation corresponds to a constant value of the fabric tensor across the entire system. Setting δ ij = 0 in equation (8) gives The left-hand side of (9) is the zero order of the operator L, L (0) =¯ :H and, for consistency, the right-hand side source term is named f (0) xy . In terms of these, equation (9) reads This is the equation obtained in [14]. To solve it, it is convenient to use the following linear transformation of variables where S ≡ −det¯ has been established to be real on macroscopic length-scales [14]. In terms of the new variables, the zero-order equation becomes Since S is real, equations (9) and (12) are hyperbolic, a fact that has a fundamental consequences. The general solution of equation (12) is where xy + B η xy η and xy + B ζ xy ζ are arbitrary functions of η = v − u and ζ = v + u, respectively. These functions are determined by the boundary data and comprise the solution to the homogeneous equation. The lines of constant η and ζ are the characteristic curves of equation (12), along which the solutions 'propagate' in the medium.
Equation (9) is not unique to σ xy . To this order, it applies to all the components σ ij , only with different source terms f (0) ij [14], given explicitly in appendix B. Thus, the zero-order solution for the entire stress tensor can be written in tensorial form, In this expression a = η or ζ, α η ≡ π xy + S, and α ζ ≡ π xy − S.
As is well known from the textbook 1D wave equation, the Green function of equation (12) for an infinite medium is where H(x) is the Heavyside step function whose value is unity when x > 0 and zero otherwise.

Example.
For illustration, suppose that the SD medium occupies the semi-infinite plane x 0 and a localized load is applied to the boundary around x = 0. For simplicity, let us take the force field g to be constant, which results in f (0) ij = 0 (see also appendix B). Because the equations are hyperbolic, care should be taken to avoid ill-posedness in the choice of the boundary conditions and simplest are The form (12) suggests that it is convenient to transform the problem from (x, y) to the (u, v) plane. In the new coordinate system, every point along the boundary (x = 0, y) corresponds to a point (u = 0, v = π xx S y). Furthermore, the boundary data become where U is the derivative of U with respect to its argument y. The solution for the zero-order stress field is then The zero-order stress response to a localized Gaussian load. The response comprises two different-magnitude stress chains that span a 'cone of influence' between them. The chains penetrate the medium at different angles that are determined by the fabric tensor.
Consider the relatively simple case V ij = 0. Changing variables in the integral on the right-hand side to τ ≡ S π xx τ and carrying out the integration, gives The solution represents penetration of the boundary stress σ ij into the system along two characteristic lines, where C η and C ζ are constants, determined by the positions of the localized boundary stresses. An example of this solution is shown in figure 1 for a bell-shaped form of U ij , localized around (x = 0, y = 0). The solution consists of two stress chains penetrating the material along the characteristic lines defined by C η = C ζ = 0, The elevated stresses along the characteristic lines give rise to forces that are significantly larger along the lines than in the medium around them. These are force chains-a phenomenon that has been the centre of attention following experimental [4] and numerical [5] observations. The magnitudes of the forces along the chains is calculated below. From (14) and (20), we observe an interesting feature of the solution-to this order, the magnitude of the stresses do not attenuate along the characteristic lines. Moreover, the linearity of the equation means that chains crossing one another do not interact or scatter. This is similar to the 1D wave equation, where two signals that propagate at different speeds pass through one another without interference. Thus, a stress chain initiating at a point on the boundary, should travel unperturbed and unattenuated through the medium regardless of other chains that it might cross along the way. Since the forces are derived directly from the stresses (see below) they enjoy the same properties. This result is in contrast to a common belief that force chains, even from different sources, may bifurcate or join [20].
Pairs of chains originate from localized stress sources, diverging as they penetrate the medium. Every such pair forms a 'cone of influence'-the analogue of the cone of light in the wave equation. The angle of the cone head can be calculated from the directions of the characteristics as follows. The unit vector tangent to the characteristic line a = η, ζ is where dy a /dx = α a /π xx and α a have been defined above. The angle between the chains can be readily calculated: The forces along the chains can be computed from the stress and the directions of the chains, f a =σ · t a . This gives for the force magnitude along the a chain where the stress in this expression is that given in the solution (14). This calculation makes it evident that the forces are also concentrated along the stress chains. These observations are useful for the analysis of many macroscopic properties, e.g. the statistics of force chains, as discussed in more detail below.

The first-order equation
To obtain the equation for the first-order stress term we need to expand the fabric parameters in equation (8) to linear order in the δ ij . For clarity, this is done in detail in appendix B. Using the notations in appendix B, the equation for σ (1) ij is The right-hand side is a function only of known quantities-gradients of the force field g, values of the fabric parameters and expressions that depend on the zero-order solution σ (0) ij . Therefore, it can be regarded as a known source term and, in analogy to the zero-order equation, it is fit to name it f (1) xy . A glaring feature of equation (26) is that the differential operator acting on σ (1) ij is L (0) -exactly as in equation (10). This sheds light on the range of validity of this approach, as will be discussed below.
As an example, consider a constant fabric tensor with an oscillatory perturbation as follows where 1. Using the results in appendix B and a little algebra, the first-order equation is The source term on the right drives the first-order solution and it gives rise to interesting effects that will be discussed below. It should be noted that, since σ (0) xy is constant along a characteristic, the modulation of the source term is due to the variability ofP. This example will be used below to test the quality of the approximation used in this approach.

Secondary chains and stress leakage
Whether the zero-order forces are physically realizable depends on their stability under the introduction of the first-order term. For the effect of the correction to be observed at all, it must not be washed out by the contribution of the external field and therefore g is taken to be negligibly small in the following.
The first-order solution to equation (26) should satisfy zero boundary conditions, since all the boundary loading is taken care of at the zero-order. σ (1) xy is determined directly by the behaviour of the source term f (1) xy . The spatial distribution of this term is analysed in detail in appendix C and it is found that f (1) xy is concentrated only along the stress chains. Now, the solution to equation (26) can be written in terms of the Green function (16) Therefore, every source point in (26) 'propagates' along the local characteristics η and ζ that emanate from this point. Imagine a pair of zero-order stress chains originating from a point on the boundary (x b , y b ). The values of the constants C η and C ζ that define these characteristics are determined from (21) using this point. Consider a particular point A = (x 0 , y η (x 0 )) along, say, the η-chain ( figure 2). The magnitude of f (1) xy at this point is a source for two secondary chains of diminished magnitudes, linear in the δ ij . The value of C η at A is the same as along the entire main η-chain and therefore the path of the secondary η-chain follows the main η-chain. It follows that this corrects the magnitude of the stress along the main chain downstream from A. Correspondingly, the stress at every point along the η-chain is modified by the superposition of the first-order corrections from all the source points upstream from it. Since the first-order correction is linear in the δ ij , whose mean vanishes, then this correction averages to zero away from the point of origin of the main chain. Nevertheless, the stress fluctuations may give rise to an overall increase or decrease the magnitude along any one particular chain. The extent to which this is significant depends on the spatial correlations of the δp ij . For example, if the fluctuations of the fabric parameters are uncorrelated then the magnitude of the stress at point A is the zero-order value plus the following superposition, Two main zero-order stress chains are shown, spanning a 'cone of influence' between them. These chains act as sources for the first-order stress terms, which propagate into the cone of influence as secondary chains. Two arbitrary source points are shown to illustrate this effect. The magnitudes of the secondary chains may have either sign and they are reduced by a factor linear in the δ ij relative to the zero-order chains (but they are exaggerated here for visualization).
where the a kl are constants. If there are no correlations in the δ kl , the fluctuations of this sum increase as the square root of the length of the chain from its inception point to A. Different correlations in the fabric parameters would give rise to different dependence of the fluctuations on the length of the chain. Thus, the stress magnitudes along it may increase or decrease slowly. The same rationale holds for the main ζ-chain, which is modulated by the secondary ζ-chains that it gives rise to.
The above is only one effect of the first-order correction. The source at point A also gives rise to a secondary chain in the ζ-direction (see figure 2). The value of C ζ along this chain is C ζ = π xx y η (x 0 ) − α ζ x 0 = C η − 2Sx 0 . In this approximation, this secondary chain runs parallel to the main ζ-chain (although in reality it might be only roughly parallel, see discussion below) and, most importantly, it runs into the cone of influence between the two main chains. Point A is but one example-since every point along the η-chain acts as a source then the entire main chain sheds a continuous first-order field into the cone of influence. Similarly, so does the main ζ-chain. Consequently, at every point within the cone of influence the stress is a superposition of two first-order stresses originating at two specific points on the main chains, e.g. point B in figure 2. Since the correction is linear in the δ ij then the stress field within the cone of influence fluctuates around zero average with a small magnitude.
This discussion leads to an interesting conclusion. Since every point along the η-chain sends a secondary chains only into the cone of influence then all the first-order field must be confined to this region. Thus, to this order, the stress outside the cone of influence remains identically zero. Recent analysis [19] shows that this conclusion goes beyond the current approximation and is valid to media with general fabric tensors.

Many-chain and extended solutions
The above analysis offers possibilities beyond the force chains solutions. To illustrate its potential, consider the stress response to an extended boundary load. The corresponding boundary data are functions U ij and V ij that are spread over a significant length of the boundary. The terms localized and extended are relative, depending only on measurement resolution, and it is in this sense that extended loads are considered. The analysis of such boundary data is made straightforward by the linearity of equation (10). Let us represent the functions U ij and V ij as sums over localized contributions (relative to the resolution of the measurement), e.g. in the form where y m are discrete points distributed along the boundary, around whichŨ ij,m (y − y m ) and V ij,m (y − y m ) are localized. The number of such points can be either finite or infinite and, correspondingly, so can be the sum in (31). The linearity of equation (10) means that the solution for the stress field is a superposition of the individual chain pair solutions that emanate from each y m . For example, consider the case study where a ij,m are the coefficient of such an expansion. A load source at y m sends a pair of chains into the medium along characteristics of constant values C ζ m = C η m = π xx y m and, from (22), the equations of these lines are Using expression (20), the solution is the following superposition of non-interfering stress chains Visual evidence for this type of solutions, albeit tentative (see discussion in the concluding section) can be found in several reports [4]. A distinct advantage of having an analytical multi-chain solution is that it opens possibilities for a range of explicit calculations. For examples, it can be used to investigate the mechanical behaviour due to given distributions of boundary loads, it can be used to analyse the statistics of spatial force distributions in the material, etc. For illustration, given a particular form of the boundary data, the stress-stress correlation function is a kl,m a ij,n 1 + π xy S e −(y−y ζ (x)) 2 /2 2 m + 1 − π xy S e −(y−y η (x)) 2 /2 2 m × 1 + π xy S e −(y+y −y ζ (x+x )) 2 /2 2 n + 1 − π xy S e −(y+y −y η (x+x )) 2 /2 2 n .
The form of this function represents a signature of the force chain network, if there is one, and can be used to distinguish between solutions with different statistics and, in particular, between force chain solutions and the more uniform solutions of elasticity theory. Expression (34) also makes it possible to compute force statistics, a quantity that is amenable to direct experimental measurements. For example, consider a section within the material, made parallel to the y-axis at some value of x. The local force fluctuates as a function of y and a useful characteristic of these fluctuations is the force-force correlation function where e x is a unit vector in the x-direction and L y is the length of interest. This function can be evaluated by substituting for the expressions forσ(x, y) andσ(x, y + y ) directly from (34). It represents the force-force correlation along the section. For SD systems this function should be statistically independent of x for all x L y min(α a /π xx ). These examples illustrate the potential of the analytic many-chain solution. A detailed analysis of the stress statistics is beyond the scope of this exposition and will be reported elsewhere.

Discussion and conclusion
To conclude, this paper presents new theoretical developments in 2D SD packs of grains. It has been shown that the field equations, which include a previously derived stress-structure relation, can be decoupled into independent equations for the three stress components. The equations are integro-differential and hence difficult to solve in closed form. Instead, they have been analysed for systems whose coarse-grained fabric tensors can be described by small fluctuations about a macroscopic volume average. For these systems, it is possible to write the equation for the linear correction to the stress field and analyse it in detail.
The zero-order solution, which corresponds to a constant fabric tensor, has been investigated in some detail. To this order, the response to a localized stress source is a pair of stress chains emanating from it. Due to the linearity of the equations, the response to a collection of localized sources is a superposition of such pairs of chains. The magnitude of the stress is concentrated along a chain and it is proportional to the magnitude of the load source. The trajectories that the chains take and the opening angle between members of the pairs have been derived as functions of the fabric parameters. The stress chains give rise to forces that are concentrated along the same chains and it is therefore natural to interpret them as the continuous representation of the force chains that are observed frequently in granular systems. The forces along the chains have been found explicitly in terms of the boundary loads and the continuous stress field.
The stability of the zero-order chain solutions under the introduction of the first-order term has been investigated. It has been shown that the source term of the first-order equation is strongly concentrated along the zero-order chains. One important implication of this result is that the firstorder solution modulates the magnitude of the stresses along the chain in a perturbative manner. The mean of the modulation vanishes for long chains, but the fluctuations from chain to chain can be quite large. Therefore, although the mean stress along a chain is the zero-order value, the stress fluctuations increase or decrease down the chain in a way that depends on the spatial correlations of the fluctuations in the fabric parameters. For example, for uncorrelated fluctuations of the fabric parameters, the fluctuations of the stress down the chain of length l increase as √ l.
Another implication of the spatial distribution of the first-order source term is that every pair of stress chains 'leaks' a first-order fluctuating stress field into the cone of influence between them. The spatial mean of this field vanishes and, unlike the magnitudes along the chains, its fluctuations remain small. It has been shown that, except for effects of the external fields g, the zero-and first-order terms cannot give rise to any stress outside the cone of influence. If this conclusion can be extended to more general systems then it has very significant implication on the spatial distribution of the stress field and its statistical properties.
It is important to discuss the range of validity of the approach taken here. A glaring feature of the first-order equation (26) is that the left-hand side hyperbolic operator is L (0) and therefore that it does not change the direction of the characteristic lines. This is a limitation of the approximation because a local fluctuation in the fabric tensor should also modify the characteristics orientation. Indeed, it is known that a continued expansion in powers of the δ ij does not converge to the correct solution [20]. Nevertheless, for small perturbations, this approach approximates well the corrections to the zero-order solutions in several aspects: it explains the change of the magnitude of the stresses along the characteristics, it explains the 'shedding' of secondary chains into the cone of influence, and it predicts the vanishing of stresses outside these cones. Note that the cones of influence discussed here are not necessarily confined by straight lines but rather by pairs of meandering characteristics. Moreover, the zero average of the fabric parameters on the granular scale ensures that the characteristics fluctuate around the predicted main paths, determined by p ij . Therefore, the approximation to first order, while failing to predict the exact meandering of the characteristics around the mean paths, is still useful for quantitative predictions for small fluctuations.
As a test of these results and their range of validity, an exact numerical solution of the stress equations (2) was carried out by Professor M Gerritsen for the fabric tensor given in (27). Applying a narrow Gaussian loading σ xy at the boundary, the full solution for the stress field has been computed for two perturbations, = 0.1 and 0.25. The solutions are plotted in figure 3. The solution indeed validates the periodicity that the source term of equation (28) gives rise to. Moreover, it establishes that: (i) the characteristics fluctuate about the zero-order straight lines; (ii) that a cone of influence exists outside of which no stress propagate; (iii) the characteristics leak low-amplitude secondary chains into the cone of influence. Even when is as large as 0.25 ( figure 3(b)) the deviations of the path fluctuations from the straight line are relatively small. Thus, the numerical solution gives evidence not only to the usefulness of this approximation but also that some of the predictions hold well beyond the first order.
The analysis presented here clarifies an apparent misconception in the literature. The directions of the stress chains have nothing to do with the directions of the principal axes of the stress tensor. Ultimately, this is because the directions of the characteristics are determined by the local fabric tensor and not by the stresses. The entire boundary stress is supported by the characteristics, not only its principal axes. This can be understood clearly for the case of a constant external field g. Applying a stress at the source such that it is given in principal axes, say σ xy = 0, then σ xy remains zero along the characteristics, regardless of the its direction. This direction need not coincide with any of the original principal axes. Moreover, it can be seen from the solution found above that only a non-constant external field can rotate locally the principal axes, namely, induce a non-diagonal component in the stress tensor through a finite source term f (0) xy = 0. However, even if this happens, there is no reason for this term to align any of the stress axes in the direction of the characteristic. At the zero-order, force chains exhibit two features: (i) forces do not attenuate along chains and (ii) on crossing, chains pass through one another without change in shape or direction, much like single solitons. The first feature starts to deteriorate with introduction of perturbations, as evidenced in figure 3, but the deterioration is gradual. The second feature remains true beyond the first order. This is due to the linearity of the equations, which allows for the superposition of boundary loads. These results provide understanding on the interactions between force chains beyond the common belief in the community that force chains can scatter from one another on incidence or join together [19]. This belief has led to models that assign the force chains random, or quasi-random, directions at every contact point along their paths. Such models give rise wrongly to parabolic chain trajectories-another manifestation of the difference from the hyperbolic chain behaviour.
These results provide a way to test experimentally which description of the stress field is more suitable. If chains scatter from one another then they should do so on every intersection. But scrutiny of the experimental images of visualized force chains, e.g. in [4], lends tentative support to isostaticity theory-force chains do appear at times to pass through one another. Additionally, the experimental observations seem to exhibit small attenuation before they eventually terminate. It is important to remember, however, that the analysis carried out here holds for marginally rigid (isostatic) states, while most measurements are made on granular systems that are not exactly so. Rather, the images of force chains in the literature have been taken from granular systems whose mean contact numbers, z, were normally higher than the marginal rigidity value z c . As discussed elsewhere [14], such systems are probably only partly isostatic, containing non-isostatic regions whose sizes depend on the difference between the mean contact number and z c . In such regions, the hyperbolic stress equations do not apply and one can regard these as elastic 'defects'. Due to the local high connectivity, chains incident on elastic regions disperse and their magnitudes fall below a pre-assigned threshold, in which case they are deemed to have terminated. This leads to a finite typical length of force chains, which gets shorter as z − z c increases. If two chains are incident on a relatively small elastic region then it may seem that they scatter from one another, while in fact they scatter from the non-isostatic defect. Moreover, if the connectivity between the region and its environment is not much larger than the isostatic value then chain dispersion is limited and they emerge from the region seemingly scattered in different directions. Therefore, it is essential that tests of the results found here be done in systems constructed carefully close to isostatic states. This can be done either by consolidating marginally rigid granular structures, as done in [22], or by the application of controlled shears to existing packs to lower the mean coordination number to values that are sufficiently close to z c .
The analytical many-chain solutions obtained here can be useful to analyses of the statistics of force chain networks, a frequently observed phenomenon that has been discussed much in the literature [23], but which has defied a good theoretical understanding (but see also [24,25]). It has been demonstrated that an analytical many-chain solution can be used to obtain explicit expressions for the stress-stress and force-force correlations in the system. An in-depth investigation along these lines, taking into considerations the limitations of the firstorder approach, is under way and will be reported elsewhere [26]. Such statistical analyses can provide a basis for predictions that can be tested experimentally.
The discussion of the many-chain, and more extended, solutions highlights an issue that has no analogue in elastic materials. Pressing with a force F x on a flat plate of length L along the x = 0 boundary of the half-plane material described in the text, it is traditional to specify the boundary pressure F x /L. This is because in conventional elasticity theory the stress field disperses away from the boundary and force fluctuations along the boundary blur out. But in SD systems they do not! The hyperbolic nature of equation (9) propagates such fluctuations undisturbed into the material. This suggests that care should be taken in the specification of the boundary loads. Since the boundary of a granular pack is never flat then the plate presses on protruding grains differently than on their neighbours. These locally elevated forces act as localized load sources and give rise to pronounced chains. If the typical distance between such sources along the boundary is larger than the scale of resolution, or interest, then the particular distribution of chains is more significant than the mean pressure. It appears that many experimental observations of force chains fall within this category, in particular, experiments where measurements are on scales of single grains. A detailed discussion of the effects of measurement scales will be presented elsewhere [26].

Appendix A. The decoupled stress equations for the special cases p ij=0
A.1. p xx = 0, p yy = 0, p xy = 0 In this case, equation (3) relates σ xx directly to σ xy and, upon substitution into equation (1), we obtain a first-order differential equation for either of these variables. The solutions of those equations are straightforward. For example, the equation for σ xy is ∂ x (rσ xy ) + ∂ y σ xy = g x , where r ≡ 2p xy /p yy . The zero-order solution for r = constant is obtained by defining a transformation to new variables u such that ∂ u = r∂ x + ∂ y . This can be done by the substitution Under this transformation where φ(v) is an arbitrary function to be determined by the boundary data. The first-order solution for small δr = r − r is where ψ(v) is again arbitrary. From this it is possible to obtain the solution for σ xx , using equation (3) and the solution for σ yy from equation (2). Now, the stress field equations are symmetric under interchange of x and y everywhere and correspondingly g x with g y . Therefore, upon this transformation, the above analysis becomes also applicable to the special case p yy = 0, p xx = 0 and p xy = 0.
A.2. p xx =0 , p yy = 0, p xy = 0 This case can be regarded as a rotation of the local coordinates so as to align along the principal axes of the fabric tensor. The stress equations now give two coupled first-order differential equations for σ xx and σ yy . One is equation (1) and the other is ∂ x σ xy − ∂ y (q yy σ xx ) = g y , (A.5) where q yy ≡ p yy /p xx . Note that since the determinant of the fabric tensor is negative then q yy < 0. On differentiation of equation (1) with respect to x, differentiation of equation (A.5) with respect to y, and subtraction of the two, we get a second-order differential equation for σ xx ∂ xx + ∂ yy q yy σ xx = ∂ x g x − ∂ y g y . (A.6) To zero-order in the fluctuations of the fabric tensor, this equation can be solved straightforwardly, yielding two characteristic lines x ± y/ √ −q yy . Stress solutions propagate along these lines, exactly as discussed in the text for σ xy . Substituting the solution for σ xx into equation (1) gives a first-order differential equation for σ xy , which can be solved straightforwardly. Substituting both these solutions into equation (4) then gives the solution for σ yy .
Using the notations in the text, let us define the following operators.
In the first of these expressions,H is the Hessian, defined in the text. In terms of these, equation (8) for σ xy in the text reads where the function φ(y) is an arbitrary constant of integration that is introduced here for completeness, but which plays no part in the analysis. Let us expand these operators in the small fluctuations of the geometric parameters, The zero order of these operators are, respectively, L (0) = ij π ij ∂ ij =¯ :H, M (0) = 0, (B.4) N (0) = 0, K (0) = π xx ∂ x g y + π yy ∂ y g x .
Expanding now both sides of (B.2) to first-order in the δ ij , the zero-and first-order equations for The source terms f (n) xy depend only on known quantities: the fabric parameters, the gradients of the external fields and, in the case of f (1) xy , also on σ (0) xy . For example, the zero-order source terms are f (0) xx = p xx ∂ x + 2p xy ∂ y g x − p xx ∂ y g y , f (0) yy = p yy ∂ y + 2p xy ∂ x g y − p yy ∂ x g x . (B.10) Repeating this analysis for the other two stress components, it can be shown that they are governed by exactly the same hierarchy of equations, only with different source terms.