Discrete models of force chain networks

A fundamental property of any material is its response to a localized stress applied at a boundary. For granular materials consisting of hard, cohesionless particles, not even the general form of the stress response is known. Directed force chain networks (DFCNs) provide a theoretical framework for addressing this issue, and analysis of simplified DFCN models reveal both rich mathematical structure and surprising properties. We review some basic elements of DFCN models and present a class of homogeneous solutions for cases in which force chains are restricted to lie on a discrete set of directions.


1.
Introduction. Cohesionless granular materials exhibit a range of behavior that has fascinated physicists, applied mathematicians, and engineers from Coulomb, in the late 18 th century, to Schaeffer, in the late 20 th and early 21 st . It is therefore somewhat surprising that perhaps the simplest question one can ask about a granular material remains unanswered even now. Given a box full of sand, if a marble is placed on top of the sand, how is the weight of the marble distributed on the bottom of the box? The problem is not just quantitative; we do not even know the qualitative form of the pattern of vertical force on the bottom generated by the marble. We do not know, for example, whether the vertical force is largest directly underneath the marble or has a maximum on a ring. 1 Experiments suggest that qualitative features of the response may depend on the way in which the sand was put into the box or on the geometric and frictional properties of the grains. [14,15]. Experiments on 2D layers with depths up to approximately 12 grain diameters show both two-peak and single-peak responses, depending on the geometry of the grains and/or degree of disorder in the packing. [14] Experiments on 3D layers have probed the response in layers as deep as 300 grain diameters, and single peaks directly beneath the applied load were observed. [15] A theoretical framework that allows a unified interpretation of these results would be of great interest.
Let us formulate a thought experiment that focuses on the central features governing the stress pattern at the bottom of the sandbox. The fundamental problem has to do with how stress is transmitted through the bulk of the material, so it is natural to eliminate the side walls from consideration by imagining a slab of sand that is infinite in horizontal extent. It is also tempting to eliminate the effect of the bottom boundary by imagining it to be at infinity, where all stresses are infinitesimal, and calculating the stress as a function of depth far from the bottom boundary. We will see, however, that in the context of an interesting class of models the bottom boundary condition plays an essential role.
We can also eliminate from the problem the complicating effects of gravity. We consider a slab of sand between two plates that are being pushed together. The density of the sand is assumed to be high enough that the material rigidly resists the force applied by the plates. The material is then subjected to a localized force created by poking a needle through the top plate and applying a specified force to it. Note that both plates and the confining pressure are essential in this setup. Without them, some grains at the surfaces of the slab would likely be subject to outward forces that could not be balanced, so the material would not support any force at all on the needle.
A third simplification is to assume that the bulk granular material is isotropic. That is, we assume that the geometry of the grain packing does not distinguish any particular direction in space. While the local environment of any particular grain is clearly anisotropic, we assume that this anisotropy, as measured for example by the fabric tensor which characterizes the degree to which there are preferred directions for "bonds" joining the centers of grains in contact, vanishes when coarse-grained over sufficiently large volumes. This is a strong assumption, as construction history of a real granular material may well distinguish favored directions for contacts. [16] Nevertheless, there are highly nontrivial structures to be uncovered even in isotropic models and it is surely appropriate to understand them before attempting to include the complications and additional parameters associated with intrinsic anisotropies.
Finally, for simplicity, we can work with a two-dimensional system. The top of the slab is a line defined to be at z = 0 and the bottom is defined to be z = d. The slab is infinite in the horizontal x direction. We are interested in determining the response to a localized force applied at x = 0 and z = 0.
Treating the sand as an ordinary (linear) elastic medium, a standard treatment [1] indicates that σ zz (x, z) should have a single peak centered on x = 0 with a width that grows linearly with depth. 2 This treatment assumes, however, that the stress field is governed by equations derived using the continuity of a well defined displacement field and that there is a unique stress associated with each strain configuration. Neither of these is true for a granular material consisting of rigid particles that support frictional forces, so a new approach is needed.
The problem of deriving macroscopic stress equations from known microscopic or grain scale physics has proven quite difficult. An indication of just how difficult this might be can be found in almost any experiment or numerical simulation that generates images of the intensity distribution of stresses on scales of several grain diameters. [5,14,18,19] In such images, one sees immediately that the stress is concentrated on filamentary structures called force chains that support compressive stress. These chains appear to be relatively straight on scales of up to 10 grain diameters or so and give the visual impression of splitting and fusing at a variety of angles. The presence of such structures suggests that passage from the grain scale to the macroscopic stress will involve two distinct steps: we must understand how grain scale physics favors the formation of chains, and also how the interactions among chains determine the macroscopic stress field. While these two tasks must ultimately be facets of a unified theory, it may be useful to approach them separately.
Models of directed force chain networks (DFCNs) have recently been proposed to bridge the gap between the scale of individual chains and the macroscopic stress, leaving open the issue of how grain scale physics promotes chain formation. [20,21] (This type of model is also referred to as a Y Y -model. [20]). As detailed below, a "Boltzmann equation" governing the densities of chains with specified strengths and orientations can be obtained, assuming only that whether a chain splits at a given point is determined only by the local environment around that point and that environments that lead to splitting are homogeneously distributed throughout the material. An essential feature of this equation is that the unstressed solution is unstable; perturbations of it lead to divergent responses. Thus in order to calculate physically meaningful response functions, it is necessary to perturb around a prestressed state. The only case for which response functions have been calculated is the special one in which all force chains support the same magnitude of stress and lie on a 6-fold symmetric set of vectors. In that case, the response consists of two peaks whose centers diverge linearly with depth but whose widths grow like the square root of depth. In other words, the only solved case of response functions for DFCNs suggests that stress propagates along characteristic directions determined by the homogeneously pre-stressed state, in marked contrast to expectations from standard elasticity theory. [21] This propagation of stress is a feature of noisy hyperbolic equations and is therefore called a "hyperbolic response." [22] It is worth noting here that hyperbolic response is known to occur in models of isostatic packings of frictionless grains [12] and that simple constitutive relations leading to hyperbolic stress equilibrium equations have been shown to explain nontrivial features of experiments on sandpiles and granular columns. [23] One must also note, however, that direct measurements of the response seem to indicate that single-peaked response is associated with strongly disordered packings, and hyperbolic response occurs only in ordered systems or at shallow depths. [14,15] The question immediately arises as to whether the hyperbolic response of the 6-fold DFCN model is an artifact of the confinement to a discrete set of directions, or perhaps of the unique property that all chains in the 6-fold DFCN have the same strength. To investigate the robustness of the 6-fold DFCN results, one would like to find other analytically tractable cases. Here we address the first step in this direction -the identification of models for which stable solutions can be found for homogeneous macroscopic stresses.
In this paper we present homogeneous solutions to the DFCN Boltzmann equation for special cases in which chains lie along a discrete set of directions with 8-fold, or more generally 4N-fold, symmetry in two dimensions. Even the homogeneous solutions -the pre-stressed states about which we might contemplate calculating response functions -have nontrivial features. The actual calculation of the response functions these solutions is beyond our present scope.
This paper also provides corrections to the presentation fo the general Boltzmann equation of Ref. [21]. The corrections do not affect any of the calculations on discrete models reported in that paper, but do clarify some conceptual points crucial for future work on continuum models or models permitting multiple splitting angles. 2. The DFCN model and Boltzmann equation.

Definitions.
A DFCN is a collection of line segments with associated compressive forces assigned such that force balance is achieved at every vertex. Figure 1 shows a simple example that illustrates several features. First, each segment is assigned a strength, or force intensity, that is indicated in the figure by the thickness of the line. Every chain is assumed to be under compression, so the force exerted on a vertex due to a given chain is directed along the chain toward the vertex, and the vector sum of the forces at any given vertex must vanish. Second, each segment is assigned a direction, indicated by the arrow, which distinguishes its "beginning" from its "end." The physical distinction lies in the local grain configurations at either end of the chain. Certain configurations do not allow propagation of a chain through them in certain orientations and therefore give rise to force chain splitting. The relevant configurations of this type for the network shown in the figure are indicated with large circles. Other local configurations would not necessarily require a splitting event, but do permit two chains to fuse when they intersect. The relevant configurations of this type are indicated with small circles.
The directions of chains in the entire network are determined by the specification of the applied boundary forces. One must think, for now, of the force being applied at point A as coming from a needle being poked through the top plate and held at constant force. This ensures that a chain beginning at A must be present. At point B, a local configuration exists that requires the chain to end. To balance the force at B, two chains must be initiated. Hence the arrows on these to chains must point away from B.
By tracing the chains initiated by applied forces at the boundary and taking into account the type of local environments at the chain endpoints, each chain direction can be uniquely assigned, with the exception of rare cases in which fusions happen to occur precisely at places that might be mistaken for splittings. When all directions are assigned, one is likely to find some chains that end, rather than begin, on a boundary. The boundary forces required to balance these forces, shown as dashed arrows in Figure 1 must be interpreted as a response to the applied forces.
Finally, there is the possibility that a single chain can appear to have inconsistent requirements on its direction, as occurs for the chain passing through point E in Figure 1. Here two splitting events have initiated chains along the same direction that meet head on. Formally, this corresponds to a fusion event in which the outgoing chain has zero strength, which thus appears as the annihilation of two chains, and the circle marking this fusion can be imagined to lie anywhere along the chain. The probability of such configurations is determined by the grain size, or, more precisely, by the ratio of the grain size to the typical separation between force chains of opposite directions. In the following, we will take this probability to be zero. (But keep in mind that a nonzero probability can be invoked to cut off the divergence in the density of weak force chains that arises in the 8-fold model.) Figure 2 further illustrates the difference between chain splitting and chain fusion in the context of a hexagonal array of grains. In panel (a), two forces are applied at the top boundary and the forces propagate along chains, crossing at the central grain. In panel (b), a similar situation is shown, but this time one of the chains is assumed to be specified by fixing its position at the bottom boundary instead of the top. In this case, a fusion of the two chains at the central grain is permitted (though not required) and the resulting configuration may be different. Panel (c) shows a local configuration (a missing disk) that requires the splitting of an incoming chain.
Thus each chain is characterized by an intensity, f , and a directionn, or, equivalently, a vector f . Force balance at a vertex requires that f for incoming chains equals f for outgoing chains. Positive values of f correspond to compressive stress along a chain. negative values to tensile stress. To model noncohesive materials, which do not support tensile stress, we take f to be always positive definite. Thus there is no ambiguity introduced by using f to denote the pair (f,n). When convenient, we will use the angular variable θ to indicate the direction ofn. We will use the term " f -chain" to refer to a chain with the given strength and direction.
2.2. The Boltzmann equation. The following discussion supersedes the treatment given in Ref. [21]. In that paper, mistaken reasoning was used in writing the continuum equation. The mistake does not affect the bulk of that paper, since the forms obtained there for networks of chains confined to a discrete set of directions were the same as those derived below. Nevertheless, the corrections made here may be important for future work on continuum models and are necessary for conceptual clarity. Note especially that the definition of P has different dimensions here from that used in Ref. [21]. The theory is developed here explicitly for the case of two dimensions. Generalization to higher dimensions is straightforward.
Let P (f, θ, r) represent the density (in an ensemble average) of force chains of intensity f and direction θ passing through the spatial point r. In other words, is the number of chains with intensity between f and f + δf that cross a unit length line segment passing through r and perpendicular toû. With this definition, P ( f ) is defined with respect to a uniform measure in the 2D space of possible forces.
Given P ( f , r), the local macroscopic stress tensor is given by [20] We wish to construct an equation that describes the ensemble average of the spatial variations in the chain densities for a system subject to specified boundary conditions. To do so, we consider the probability that a chain f will exist at a point r + ǫn, assuming we know P (f,n, r).
A given material will be characterized by a two scalar parameters, λ and Y , and two angular functions, φ s and φ f .
• λ is the average distance in any specified direction between the point where a chain begins and the nearest point that will cause it to split; i.e., the mean length chains would have if there were no fusions.
is probability that a given splitting of a f 1 -chain results in chains with strengths and directions given by f 2 and f 3 , normalized to unity in the sense defined below.
is relative probability that a f 2 -chain and a f 3 -chain will fuse when they intersect and thereby form a f 1 -chain, normalized to unity in the sense defined below. • Y is an overall efficiency with which two intersecting chains will fuse. That is, the total probability that two intersecting chains will fuse is Y φ f φ s and φ f must both include delta functions that enforce force balance at each intersection and Heaviside functions that guarantee all forces are compressive.
We neglect, for now, the finite size of individual grains. That is, we assume that λ is large compared to a grain diameter and treat the granular packing underlying the force chain network as a continuous medium. The equation governing the spatial variation of P as follows (for notational convenience, we have dropped the r from the argument of all of the P 's): Each term on the right hand side represents a type of event that can alter the density P ( f ). We discuss each in turn: • The first term represents the loss of f -chains due to splitting. The sin factor is exhibited explicitly for future convenience.
Here P ( f ) may be taken out of the integral. To ensure that λ is the mean distance between a chain will propagate before splitting (in the absence of interactions with other chains), we must therefore have Note that φ s has dimensions of 1/f orce 4 . • The second term represents the creation of f -chains due to the splitting of chains in other directions. The splitting function here must be the same as in the first term, but with arguments exchanged, since the rate at which chains chains appear due to splitting must match the rate at which the parent chains split. The factor of 2 counts the identical integral arising from the exchange of the prime and double-prime labels. • The third term, which is nonlinear in P , represents the loss of f -chains due to their fusions with f ′′ -chains (or f ′ -chains, hence the factor of 2). The quantity |sin(θ − θ ′′ )| P ( f )P ( f ′′ ) is the density of intersections of f and f ′′ -chains, and φ f is specifies the probability of fusion when two such chains meet. Here again P ( f ) can be taken out of the integral. We choose to normalize φ f such that the remaining integral has a maximum value of unity over the set of for all values of f 0 . In other words, we normalize φ f according to the type of intersection most likely to produce a fusion if all chain densities were identical. In equation form: or, equivalently, where the maximization is over all possible values of f and f 0 . The parameter Y then provides an absolute measure of the fusion efficiency. φ f has dimensions 1/f orce 2 , in contrast to φ s .
• The fourth term represents the creation of f -chains due to the fusion of chains in other directions. For consistency, the fusion function must be the same as in the third term.
For materials composed of perfectly rigid grains, a rescaling of all of the forces in a given DFCN yields another perfectly acceptable network. Hence, Eq. (3) must be invariant under a rescaling of all the force intensities. This is verified by straightforward dimension counting (unlike the form of the Boltzmann equation suggested in Ref. [21]).
Consider further the general forms of φ s and φ f . For an isotropic material, φ can depend only on the differences between angles. The general form is where F a is a combination of the force intensities that provides the correct dimensions for φ a , and ψ a is a symmetric function of its arguments. (One might imagine functions of the ratios of force intensities multiplying the arguments of ψ a , but these will always reduce to functions of the angles alone due to the δ function.) In the case of fusions, dimension counting in Eq. (4) implies F f is dimensionless and hence equal to unity, up to a constant that can be absorbed in Y . In the case of splitting, on the other hand, dimension counting in Eq. (4) implies F s has dimension 1/f orce 2 . The relations between the f 's enforced by the δ function guarantee that any combination of f , f ′ , and f ′′ with the right dimensions is equally valid; the differences between them can simply be absorbed into ψ s . The natural choice is F s = 1/(f ′ f ′′ ). With this choice a constant ψ s = 1/π 2 corresponds to an equal probability for every possible splitting configuration.
It is clear that λ can be scaled to unity without loss of generality by choice of the unit of length. Though Y is a dimensionless parameter with physical significance, from a mathematical point of view it plays a trivial role in Eq. (3). To solve for P for any given value of Y , we solve the case Y = 1, then simply divide all P 's by Y . From here on, we take λ = 1 and Y = 1.
2.3. Specialization to discrete directions. We now make two simplifying assumptions to arrive at a set of ordinary differential equations that permits analytical solution. First, we assume that φ s and φ f are nonzero only for vertices of the forms shown in Figure 3. This means that all chains lie in the discrete set of directions θ n = (n − 1 2 )π/(2N ) (measured from the positive z direction, which is downward), and that the strong force at a given vertex is always related to two weak ones by a factor of ξ = 2 cos(α). Defining f m ≡ f 0 ξ m , the continuous function P (f, θ) becomes a set of discrete functions P (f m , θ n ), which will be denoted P n,m δ( f −f mnn ). In this notation, it is always assumed that n is taken modulo 4N . Note that P n,m is simply a number per unit length; the dimensions of force are taken care of by the two-dimensional delta function. We will treat the case α = π/(2N ) explicitly here. Generalization to any α that is an integer multiple of this is straightforward.
This first assumption corresponds to the following forms for the angular parts of the splitting and fusion functions: The difference in character between φ s and φ f stems from the fact that φ f appears in integrals that are quadratic in P . When P is a sum of two-dimensional delta functions, the fusion integrals have two more delta function factors than the splitting integrals. The same can be seen in the normalization conditions of Equations (4) and Equations (6). Second, we assume that the boundary conditions of interest are uniform across the top and bottom of the slab. Translational symmetry in the x-direction then dictates that all horizontal gradients of P n,m vanish.
Inserting these assumptions into Eq. (3) and integrating both sides over a small volume of force space in the vicinity of f mnn , we obtain the following ordinary differential equations for the P n,m 's: (cos θ n ) ∂ z P n,m = −P n,m + P n+1,m+1 + P n−1,m+1 (10) +P n+1,m−1 P n−1,m−1 − P n,m P n−2,m − P n,m P n+2,m .
The result is straightforward, but one must be careful to account for all the trigonometric factors arising from integrating over the delta functions, including the integral on both sides that is required to isolate P n,m .
The different terms on the right hand side of Eq. (10) correspond to the processes described above that can create or destroy chains. In the present case, the first term represents the decay of P n,m along the direction θ n due to the splitting of an existing (n, m) chain. The next two (linear) terms account for the addition to P n,m due to the splitting of other chains in adjacent orientations. Note that if a chain is to split and create a contribution to P n,m , which must have strength f 0 ξ m , it must begin with strength f 0 ξ m+1 . The first nonlinear term accounts for events in which two chains of strength f 0 ξ m−1 fuse, and the last two nonlinear terms to events in which chains contributing to P n,m are deflected due to fusions with other chains.
In Reference [21], Eq. (10) was studied in detail for the case of 6-fold symmetry. In that case we have ξ = 1, so the hierarchy of equations indexed by m collapses. In addition, cos θ 2 = cos θ 5 = 0, so the entire system reduces to four ODEs and two algebraic relations. If we assume mirror symmetry about the vertical axis, these are reduced to two equation that can be solved completely.
In order to observe the variations of P with force strength, we must consider symmetries other than 6-fold. To avoid technical complications arising from the presence of strictly horizontal chains having cos θ = 0, we consider only the cases of 4N -fold symmetry for positive integers N . The simplest case is that of 8-fold symmetry, and example of which is illustrated in Figure 1.
Note that the discreteness of the possible orientations of force chains does not imply any discreteness in the possible positions of the splitting or fusion vertices. Note also that the system remains isotropic in the sense that Eq. 10, including all coefficients, is identical for all n, except for the geometric factor on the left-hand side that results from restricting attention to variations in z and not x.
3. Homogeneous networks for symmetric splittings and fusions. General solutions of Eq. (10) are not yet known. We study here the important special case of homogeneous solutions; i.e., the fixed points, for which all gradients vanish. These solutions are both nontrivial and crucial for understanding the response function. Consider the trajectory of the vector P n,m as z is varied. As detailed in Reference [21] for the 6-fold case, trajectories that do not pass very close to a stable (or marginally stable) fixed point inevitably lead to divergences at z >> λ that are reflected in unphysical negative values of some P n,m . Trajectories that avoid these divergences do so by staying in the vicinity of a fixed point over most of the range of z. As also noted in Reference [21], the trivial fixed point P n,m = 0 is unstable and hence cannot serve as a reference state for a linear response theory. Thus the response function we ultimately seek will be dominated by the DFCN structure corresponding to some nontrivial fixed point.

Isotropic solutions.
It is instructive to begin with a search for isotropic solutions; i.e., solutions for which P n,m = P m for all n. The fixed point equation derived from Eq. (10) reads To simplify this recursion relation, define Eq. (11) says 2d m+1 − d m = 0, which implies For d 0 = 0, Eq. (12) gives For d 0 < 0, this recursion relation either produces a negative numbers or a divergence at large m, neither of which is physically acceptable. If P 0 is sufficiently small, the repeated squaring drives P m toward zero at least as fast as p −2 m with p > 1. The negative contribution from the d 0 term can only speed up this decay as long as P m is positive. But this means that P m decays faster than 2 −m , so the d 0 term eventually makes P m negative. Now if P 0 is large enough to avoid P m decaying toward 0, repetitive squaring causes it to diverge like p 2 m with p > 1; the d 0 term becomes irrelevant at large m. For d 0 > 0, on the other hand, problems arise for large, negative m. We write P m−1 = √ P m − d 0 2 −m and see that complex values will be generated unless P m is always greater than d 0 2 −m . But this is impossible because P m−1 is strictly less than P m .
Thus we are left with d 0 = 0 as the only physically relevant case. When d 0 = 0, Eq. (10) is satisfied in a special way: both −P m + P 2 m−1 = 0 and 2P m+1 − 2P 2 m = 0 are satisfied simultaneously. These two relations are actually identical up to a shift in m, and admit a one-parameter family of solutions where we must have p > 1 to avoid divergence at large m. Despite its apparent simplicity, it is instructive to examine this solution in a bit more detail. We consider first the behavior at large negative m, then the behavior at large positive m. Note that P m → 1 for large negative m. This means that the densities of all chains of intensity f 0 ξ −m for large m are the same. In other words, the total density of force chains, m P m diverges. Though this may appear troublesome, it does not yield a divergence in the pressure: m P m f m remains perfectly finite. The divergence in chain density is an artifact of the restriction to a fusion function φ f that does not allow weak chains to fuse with stronger ones. Removal of this artifact requires working with a continuous distribution of force chain directions and intensities and it beyond the scope of this work.
For large positive m, we can compute the probability distribution for contact forces, a quantity that has been shown in a variety of experiments and numerical simulations to decay exponentially for large forces. For the homogeneous, isotropic solution, we have where f m = f 0 ξ m and P m = p 2 m . The denominator in this expression is just the spacing between points with successive m's, the required conversion factor from the discrete density P m to the density per unit force. This can be rewritten as Recall that p > 1 and ξ < 2, so β > 1 and P decays faster than exponentially for large f . For the 8-fold case, we have β = 2 and hence a Gaussian decay. But smaller splitting angles give ξ closer to 2, and for a splitting angle as large as 2α = 60 • , we have β ≈ 1.26. In this case, the full distribution (including the 1/f term) can look surprisingly like a simple exponential over several decades, as shown in Figure 4. Experience with the 6-fold model suggests that the isotropic fixed points of Eq. (10) are non-generic. In the 6-fold model, linearization of the theory in the vicinity of a fixed point generally leads to hyperbolic response functions with peaks that are sharper for more strongly anisotropic fixed points. The isotropic case is pathological in that the peak width is divergent [21] and additional cancellations conspire to produce single peak, but one that obeys different scaling laws from the predictions of elliptic models. [24] It is therefore important to consider solutions of Eq. (10) that are not constrained to be isotropic.

Anisotropic solutions.
To find fixed points of Eq. (10), we begin with the following ansatz : where p is a real parameter and q m , x n , and ζ are to be determined. (This guess combines features of the isotropic solution above and the 6-fold solutions discussed in Reference [21].) Substituting into Eq. (10), we have Taking a cue from the isotropic case, we consider the possibility that the cancellation required by this equation occurs term by term in the following way: 3 Note that all three of these equations are identical up to shifts of n and m. Since the q terms are independent of n, we obtain two recursion relations: Eq. (23) implies the free parameters being p, q 0 , and x 0 . Recall that the force intensity associated with the chains whose density is P n,m is f m = f 0 ξ m and note that ξ = ζ 1 . Note also that for j = 0, we have ζ = 2, so that the x 0 and q 0 terms can be combined. This leaves us with precisely the same isotropic solution discussed above. Eq. (28) requires further examination to determine the range of different structures it represents. First, each three-parameter family is really only a one-parameter family of distinct physical networks. There are ways of rescaling the parameters that have no effect on the physical system begin described: The first of these corresponds to a redefinition of the force scale, f 0 → f 0 ξ ν , which cannot affect the physics. Strictly speaking, the system is exactly invariant under this transformation only if ν is an integer. Arbitrary values of ν, however, merely shift the positions of the discrete set of P n,m 's along a single smooth curve for each n, with integral ν's being the values that shift the entire set into itself. For our purposes, the differences between solutions related by non-integer ν are unimportant details. (See Figure 5.) The second transformation leads to a mathematically identical solution. (The case q 0 = 0 is not physically relevant, as it leads to divergence of the total stress from large m force chains.) By performing the first transformation with ν = log(x 0 /q 0 )/ log(2/ζ), followed by the second, we can adjust q 0 and x 0 both to unity. Second, note that we cannot take arbitrary linear combinations of the solutions to Eq. (24), since they would not satisfy the nonlinear Eq. (10). We can, however, take linear combinations of eigenvectors that have degenerate eigenvalues. Since ζ 4N −j = ζ j for 1 ≤ j ≤ 2N − 1 we can generate solutions from arbitrary linear combinations of the degenerate pair of eigenvectors y (j) and y (4N −j) . Thus j = 0 and j = 2N each provide a single continuous family of solutions parameterized by p, while all other j's come in pairs that provide families parameterized by p and by the relative weights of the two y (j) 's.
Consider now the family spanned by y (j) and y (4N −j) , for some specific j. By forming linear combinations, we can construct basis vectors e j that are symmetric or antisymmetric under reflection through a vertical axis. The symmetric combination will give solutions for which P n,m = P 4N +1−n,m , as might be expected if the average confining force on the plates is purely normal. The antisymmetric combination will give solutions that break this symmetry, though they do not yield any negative values of P n,m . These might be important for describing systems subjected to shear as well as compression. The explicit forms of the symmetric combinations (including the trivially symmetric isotropic case) can be expressed as where the normalization is || e j || = 4N for all j.
The antisymmetric combinations, normalized to 4N , are given by In the following, we focus on the symmetric cases only. The e (0) solution is the isotropic case discussed above. The e (1) solution describes a stress configuration with a dipole-like anisotropy. At each generation m, the angular variation P n,m has a single maximum (for p > 1) at n = 1 and a minimum n = 2N . There are more downward-pointing than upward-pointing chains for each undirected orientation. Moreover, because ζ 1 is positive, the decay of P n,m with m is monotonic both at large m, where it decays to 0, and at large negative m, where it decays to 1. Figure 5 shows the form of P n,m for this case. Note that the term "dipole-like" applies here to the P 's. Since the stress field does not distinguish between force chains in opposite directions, the anisotropy in the stress is "quadrupole-like" -stronger in the vertical directions than in the horizontal. A striking feature of the solutions is the occurrence of peaks for some n, which become more pronounced for larger p. Peaks are evident for all n when the stress, f m P n,m , is plotted rather than the chain density. (See Figure 6.) Differentiating Eq. (28) with q 0 = x 0 = 1, we find that the peak in P n,m occurs at m = log y (j) n + log log |ζ n | − log log 2 log 2 − log |ζ n | , which is independent of p. (The peak in f m P n,m does depend weakly on p.) The force scale represented by this peak must have its origin in the boundary conditions, since a linear rescaling of all forces in the bulk has no effect on the DFCN structure. Recall that part of the process required to collapse the solutions onto the oneparameter family indexed by p was to rescale f 0 , which corresponds to a rescaling of the strengths of the force chains injected at the boundaries. There are subtleties lurking here, however, as will be discussed further in Section 3.3. A full calculation of the single contact force distribution P(f ) is slightly more complicated than for the isotropic case, as it requires a sum over the different chain orientations. It is not hard to extract the dominant term at large f , however. For the symmetric j = 1 case with q 0 = x 0 = 1, Eq. (28) can be written as where β = log 2/ log ξ as above, and γ = log ζ 1 / log ξ = 1. Since β > 1, the dominant term at large f will be P 1,m = p −f β m +f . It is important to note that the anisotropy of the solution is not a consequence of any material anisotropy. The bulk properties of the material we are describing are perfectly isotropic. The boundary conditions on the system will generally not be isotropic, however, so there is no reason to expect that the pre-stressed state about which the response function should ultimately be calculated should have isotropic stress. The generic situation described by the discrete DFCN theory presented here is one in which an intrinsically isotropic material is subjected to loading forces that generate anisotropic distributions of stress chains. Because different degrees of anisotropy correspond to different fixed points of the full nonlinear theory and hence to different linearizations, the pre-stressed state may have a profound influence on the response function. This is in contrast to an ordinary elastic medium, in which the response function is not affected by an anisotropic background stress.
The solutions for higher j have features that seem a bit unlikely from a physical standpoint, though not obviously unacceptable. For 2 ≤ j < N , we have higher orders of anisotropy, with P n,m exhibiting additional maxima and minima as a function of n at each m, but still with monotonic decay at large m. For N < j ≤ 2N , the eigenvalues are negative, which causes ζ m and hence P n,m to oscillate as m varies. Whether or not these solutions are relevant hinges on the question of which boundary conditions are associated with generic physical systems.
3.3. Remark on boundary conditions. Given several families of homogeneous solutions, it is natural to ask which solution will be selected for a given set of boundary conditions. From a mathematical perspective, this is a straightforward question. As described in detail for the 6-fold case in Reference [21], the only way to specify boundary conditions on Eq. (10) that are sure to be self-consistent and not produce negative chain densities is to specify the densities of only the ingoing chains at each boundary. In general, the densities specified will not match perfectly any of the homogeneous solutions P n,m for all m and n, even if horizontal translational symmetry is assumed. Based on experience with the 6-fold case, we might expect the trajectory in the space of P n,m as the slab is traversed to quickly approach a fixed point and stay near it over a substantial portion of the slab, then make a rapid transition to the vicinity of another fixed point, and finally move off to a point consistent with the bottom boundary condition. The details of the transition regions and selection of the fixed point remain to be worked out in the present case.
From a physical perspective, the situation is not so clear. Consider the following two possible choices of boundary conditions for the top surface. In both cases, assume that P n,m = 0 for all ingoing (n, m) except the ones specified as follows: (BC1) P 0,5 = P 1,5 = 1; and (BC2) P 0,1 = P 1,1 = ξ 5 . Assume also that the boundary conditions on the bottom slab are simple reflections of these through the horizontal. Thus we have boundary conditions with the same total stress stipulated, but in one case it is injected via a sparse set of strong forces, and in the other via a denser set of weaker forces. These conditions will not necessarily lead to the same fixed points in the bulk. Since the chain densities are measured with respect to the intrinsic length scale λ, BC1 and BC2 are not simply related by a scale transformation. We therefore expect them to lead to fixed points corresponding to different p's, which implies that they will require different rescalings of f 0 . This in turn indicates that the force scales fixed by the two different boundary conditions differ, in spite of the fact that they correspond to specifying the same overall force.
The resolution of this apparent paradox lies in the fact that the boundary conditions do not specify the densities of the outgoing chains. Ultimately, the force scale must be determined by the zz component of the stress, which must be constant throughout the slab; i.e., σ zz must determine the prescribed rescaling of f 0 . But the densities of outgoing chains produced by BC1 and BC2 will differ, which will cause σ zz to differ in the two cases. In order to determine σ zz , which is the total pressure applied to the top surface, we must first determine the full solution, or at least which fixed point the trajectory approaches.
The converse situation, in which the stress is specified but not the individual chain densities, is more familiar, but just as subtle. One might expect it to be sufficient to specify some components of the stress tensor and its spatial derivatives on the top boundary and others on the bottom, but the DFCN equations require independent specifications of densities of chains of many directions and strengths. In an experiment in which a granular material is put under compression between to flat plates, it is not at all clear how the chains densities at the boundaries are determined. Both the distribution over directions n and strengths m are relevant. An important topic of future research will be to develop a method of assigning directions to force chains observed in experimental images or numerical models, and to learn from them how the boundary conditions appear to be imposed. 4. Conclusion. In this paper we have exhibited for the first time solutions of the Boltzmann equation for the chain densities in a model that permits the interaction of chains of different intensities. These solutions have highly nontrivial structure that raises a number of interesting questions. Most importantly, they permit calculations of both stresses and single particle force distributions in a unified way. Refinement of the model is clearly necessary: we would like to avoid divergences in the density of weak chains; and we would like to generalize the splitting and fusion functions, which might bring the single contact force distribution P(f ) into closer agreement with experiment. Nevertheless, there is much to be gained from further study of the models presented here, or minor modifications of them.
The availability of closed form analytic solutions suggests that a complete theory of the stress configurations and response functions for these particular models is within reach. In addition, it appears possible that the linearized theory in the vicinity of the fixed points presented here could lead to a new derivation of constitutive relations for closing the equilibrium equations for the macroscopic stresses. Finally, the solutions presented here may provide a basis for extensions to similar discrete models in which more than one type of vertex structure is allowed. All of these topics are currently under study. [25] The author looks forward to continued collaboration with David Schaeffer on these topics, as he is certain to offer valuable insights and pose questions that lead to fruitful calculations.