# Force distribution in a scalar model for non-cohesive granular material

## I Introduction

Attempts to understand stress distribution in static, non-cohesive, granular materials have uncovered a rich structure of force chains that currently is poorly understood . Although mechanical systems such as bead or granular packings subject to uniaxial compression seem straightforward, they harbor fundamental theoretical challenges concerning the connection between the micro and macro scales. Such systems are far from thermal equilibrium since thermal energy scales are negligible compared to the potential energies stored in elastic deformation of the grains. In addition, the force between grains is strongly nonlinear, as it vanishes identically when the grains are not in contact.

In this paper, we report results and analyses of a toy model for the compaction of a granular material. Our model consists of “grains” placed on a square lattice and connected by vertical springs. The grains are constrained to move only in the vertical direction, so the model treats forces and stresses as scalar quantities. (See Fig. 1.) Geometric disorder is introduced in the model by assigning to each spring an equilibrium length drawn randomly from a square distribution. Though we cannot hope to capture the true physics of tensorial stresses in this manner, we can study the validity of arguments that have been applied to real granular systems and should equally well apply to our toy model. In particular, we compare the distribution of inter-grain forces generated by our model to the distribution predicted by the $`q`$-model of Coppersmith et al , and we compare the stress-strain power law generated by our model to that expected on the basis of a mean-field argument. Our results provide a warning: in neither case are the theoretical predictions borne out by the numerics.

## II Definition of the model

The structure of our model is shown in Fig. 1. Grains are placed at the even-parity sites of a $`2N_x\times (N_y+1)`$ cubic lattice with lattice constant $`b`$. We refer to such a system as an “$`N_x\times N_y`$” system. Excluding the sites on the lower boundary, it contains $`N_xN_y`$ grains. Grains are constrained to move only in the vertical direction, as if sliding on fixed, frictionless wires (shown as thin lines in the figure). We let $`u_{i,j}`$ denote the vertical displacement of the grain at lattice site $`(i,j)`$ from its nominal height $`jb`$, with $`u_{i,j}>0`$ for upward displacements. Between every two nearest-neighbor grains there is a vertical spring. All of the springs respond linearly under compression, with spring constant $`k`$, but they produce no force under extension. The energy of the spring between grains $`i,j`$ and $`i\pm 1,j1`$ is taken to be

$$E_{ij}^\pm =\{\begin{array}{cc}\frac{1}{2}k\left(u_{i,j}s_{ij}^\pm u_{i\pm 1,j1}\right)^2\hfill & \text{if}(u_{i,j}s_{ij}^\pm u_{i\pm 1,j1})<0\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(1)   
where $`s_{ij}^+`$ and $`s_{ij}^{}`$ are independent, quenched random variables designating the different uncompressed lengths of the springs. We take $`s_{ij}^+`$ and $`s_{ij}^{}`$ to be distributed uniformly over the interval $`(ϵ,ϵ)`$. Note that when $`s=0`$, the uncompressed spring has length $`b`$, the lattice constant. One may imagine the grains to have arbitrary shapes and the uncompressed spring length to be the vertical spacing between grain centers when the grains are just barely in contact.

The system is taken to satisfy periodic boundary conditions in the horizontal direction. All grains at the top are constrained from above by a rigid “ceiling” and therefore satisfy $`u_{i,N_y}0`$. On the bottom row, grains are constrained by a rigid “floor” whose height $`U`$ is an independent variable; thus $`u_{i,0}U`$. Let $`U_0`$ denote the maximum value of $`U`$ such that no springs are compressed. Typically $`U_0`$, which is a function of the random $`s_{ij}^\pm `$’s, is negative and of the order of $`N_yϵ`$ in absolute value. (It cannot be less than $`N_yϵ`$.) Below we measure the floor height in terms of the incremental variable $`\mathrm{\Delta }UU_0`$.

Though we use the terms “floor”, “vertical”, etc., for convenience, there is no gravity in the model and there is complete statistical symmetry between the up and down directions. The algorithm discussed below for determining equilibrium configurations explicitly breaks this symmetry, though the configurations it produces do not.

### A Meaningful parameters

Our model apparently depends on the parameters $`ϵ`$ and $`k`$, but these merely set the scales for distance and force. For the rest of this paper, numerical values for distances — spring-length variations $`s_{ij}^\pm `$, grain displacements $`u_{i,j}`$, and boundary displacement $`\mathrm{\Delta }`$ — will be quoted in units of $`ϵ`$; and forces — individual forces $`f_{ij}^\pm `$ and the total force $`F`$ on the ceiling — in units of $`kϵ`$.

Although the width of the distribution of spring lengths scales out of the problem, its shape could conceivably make a difference. Nevertheless, we do not expect the precise shape to be important as long as the probability of having very large uncompressed lengths decays sufficiently rapidly. In any event, we consider only the square distribution.

Let us contrast the alternative of random spring constants with the random equilibrium lengths we study here. If the $`k`$’s were random but all equilibrium lengths were equal, then as $`\mathrm{\Delta }`$ increased all the springs would begin to be compressed simultaneously at $`U=0`$ and the system would be perfectly linear for all $`\mathrm{\Delta }>0`$; the force-vs.-displacement law would be $`F(\mathrm{\Delta })=K\mathrm{\Delta }`$ for $`\mathrm{\Delta }>0`$ and $`F(\mathrm{\Delta })=0`$ for $`\mathrm{\Delta }<0`$, where $`K`$ is a constant determined by the random $`k`$’s. By contrast, with random spring lengths we capture some of the geometric disorder of real granular materials, which inevitably introduces nonlinearities associated with the formation of new contacts during compression.

### B General features of equilibrium configurations

We will refer to a spring under compression as an “active bond”, and to a site connected to such a spring as an “active site.” An equilibrium configuration for $`\mathrm{\Delta }`$ positive, but not too large, has a mixture of active and inactive sites and bonds. For a given set of $`s_{ij}^\pm `$’s, the equilibrium configuration of active sites and bonds is unique, although the exact positions of inactive sites are not determined. This uniqueness follows from the strict convexity of the total energy as a function of the lengths of active bonds.

Fig. 2 shows three equilibrium configurations with the same set of $`s_{ij}^\pm `$’s on a $`40\times 40`$ lattice, for different amounts of compression. Only the active bonds are shown, with line thickness proportional to the force transmitted. Also shown are two bars representing the ceiling and the floor.

Let $`F(\mathrm{\Delta })`$ be the total force exerted by the system on the ceiling. Fig. 3 shows a typical force curve for a $`40\times 40`$ system, together with the fractions of sites and bonds that are activated. Elementary considerations reveal the limiting behaviors of $`F`$ for large and small $`\mathrm{\Delta }`$. For very small $`\mathrm{\Delta }`$, the force on the ceiling is due to a single chain of activated springs that reaches from top to bottom. This chain, which follows the directed path $`i(j)`$ through the lattice that yields the highest value of $`_{j=1}^{N_y}s_{ij}`$, produces a force linear in $`\mathrm{\Delta }`$:

$$F(\mathrm{\Delta })=\frac{1}{N_y}\mathrm{\Delta }\text{(}\mathrm{\Delta }\text{ small)}.$$
(2)   
When $`\mathrm{\Delta }`$ is very large, all of the springs in the system are compressed. The system is then again entirely linear: because all of the springs have the same stiffness,

$$F(\mathrm{\Delta })=\frac{2N_x}{N_y}\mathrm{\Delta }+C\text{(}\mathrm{\Delta }\text{ large)},$$
(3)   
where $`C`$ is a (negative) constant that depends on the random $`s_{ij}^\pm `$’s and is difficult to compute.

In general, as $`\mathrm{\Delta }`$ is increased from zero, additional chains are activated and $`F(\mathrm{\Delta })`$ increases. In fact it can be shown that $`F`$ must increase monotonically with $`\mathrm{\Delta }`$. (This is why $`C`$ must be negative.) Surprisingly, however, the slope of $`F(\mathrm{\Delta })`$ may decrease with increasing $`\mathrm{\Delta }`$. Occasionally bonds that are active at small values of $`\mathrm{\Delta }`$ become inactive for larger $`\mathrm{\Delta }`$, thus reducing the effective spring constant of the network. The simplest example of this behavior is illustrated in Fig. 4. The figure shows a sequence of configurations in which the central bond is originally active, but becomes inactive due to the effective stiffening of other parts of the network when additional bonds are activated.

## III The numerical algorithm

The algorithm we employ for generating equilibrium configurations of our model relies heavily on the fact that for a given set of active bonds, the system response is linear. By solving for the times at which inactive bonds are activated during compression, the algorithm generates the entire curve $`F(\mathrm{\Delta })`$, starting from $`\mathrm{\Delta }=0`$ and explicitly visiting every configuration. To construct the initial, uncompressed configuration, we begin with $`u_{i,0}=U_0`$ for all $`i`$ and march up layer by layer applying the rule

$$u_{i,j}=\mathrm{max}\{(u_{i1,j1}+s_{ij}^{}),(u_{i+1,j1}+s_{ij}^+)\}$$
(4)   
This configuration can be thought of as the packing that would result from an infinitesimal gravitational force acting on the system. (Strictly speaking, we set $`u_{i,0}=0`$, march upward, determine $`U_0`$ from the height reached, and translate all displacements accordingly.)

A typical initial configuration is shown in Fig. 5. Sites joined by line segments are separated precisely by the equilibrium length of the spring between them. Springs that are not long enough to connect their two sites are not drawn. Already one can see that there is a rich geometric structure hidden in this model, quite similar to that described by Roux et al . The heavy lines in the figure highlight structures we call “trees”, consisting of all sites that are connected to the same site on the first layer.

We refer to the distance between the upper end of an inactive spring and the site to which it will eventually connect as a “gap”. Note that sizes of the gaps reflect the structure of our initialization algorithm. We have moved all the grains to their lowest possible position, concentrating all of the small gaps that might be present in a chain into a single large gap at the top of a branch of a tree.

To generate the force curve, we maintain a list of active sites and their active bonds, and we keep track of changes in the tree structure of the uncompressed chains as the floor is raised. For any fixed set of active bonds, the $`u_{i,j}`$’s increase linearly with $`\mathrm{\Delta }`$, since all of the active springs are linear. It is therefore a straightforward linear problem to solve for the rate, $`du/d\mathrm{\Delta }`$, at which each $`u_{i,j}`$ advances. Advancing $`\mathrm{\Delta }`$ through one stage involves solving for $`u_{i,j}`$ for some arbitrary increase in $`\mathrm{\Delta }`$, for which we employ an iterative biconjugate gradient routine , using the current configuration of $`u_{i,j}`$’s as an initial guess for the solution. The solution is used to determine the rates of advance of all the active sites. The rate of advance of each inactive site is equal to that of the active site that supports it through an uncompressed chain. The rates are used to determine the value of $`\mathrm{\Delta }`$ at which the network topology changes. The $`u_{i,j}`$’s are then updated according to the calculated rates, the list of active sites or the tree structure is updated, and the process is repeated.

The updating of the network topology is necessitated either by the closure of a gap or by the breaking of a bond. The closure of a gap can generate one of two types of events. First, it can cause an additional set of sites to be activated, which increases the size of our active site list (and also increases the size of the linear problem that must be solved on the next iteration). Second, it can result in a “push” event, in which an inactive branch of one tree is simply transferred to another tree without becoming active. Pushes dominate the behavior of the system in the early stages of compression, becoming more and more rare as compression continues and the number of inactive sites decreases.

The breaking of a bond, as mentioned in Section II B, is the final possibility for changing the network of active bonds. Breaks are generally rare events. During full compression of a $`40\times 40`$ system, which undergoes approximately 2000 events (pushes, chain additions, breaks), there are typically about 5 breaks.

The algorithm reaches completion when all bonds have been activated. As explained above, further compression would be homogeneous with all $`f_{ij}^\pm `$’s increasing at precisely the same rate. To compress one configuration on a 40x40 lattice, generating one complete force curve, requires approximately one hour of computation on a typical 200MHz workstation. We have accumulated data for various lattices with at most a few thousand sites. Future efforts to optimize the algorithm should permit investigation of substantially larger lattices.

Despite our reliance on linear methods for evolving the system between gap events, it must be emphasized that the absence of tensile forces introduces a strong nonlinearity for larger increments in $`\mathrm{\Delta }`$. An alternative approach to following the entire evolution of the system would be to solve directly for the configuration of the system at an arbitrarily chosen value of $`\mathrm{\Delta }`$ by minimizing the nonlinear energy function of Eq. (1). Such an approach might speed up the calculation even with linear springs, and it would be essential if the springs were nonlinear (e.g. Hertzian) under compression.

## IV The macroscopic force as a function of displacement

### A General remarks and a “mean-field” prediction

Fig 6 displays numerically computed curves of $`F`$ vs. $`\mathrm{\Delta }`$ for a few rectangular systems of various sizes. At a cursory level, the force curves appear to have the form $`F(\mathrm{\Delta })\mathrm{\Delta }^\nu `$, which correspond to lines of slope $`\nu `$ on our log-log plots. Interestingly, this exponent depends strongly on system geometry.

Closer inspection of Fig. 6 reveals both expected and unexpected behavior. For any particular realization of the disorder, there is a region, sometimes substantial, on the log-log plot for small values of $`\mathrm{\Delta }`$ in which the force curve appears nearly linear because it is dominated by a single chain. There is also a crossover to linear behavior for large $`\mathrm{\Delta }`$ since, as the fraction of active bonds approaches unity, the system approaches the linear limit described in previous sections. The intermediate regime is the one of interest, and the behavior there is rather complex. In many runs an exponent $`\nu >1`$ in the intermediate regime can be readily identified. However, in other runs (see, e.g., Fig. 6c), small-system statistics tend to obscure the phenomena. It appears that completely reliable measurements of the exponents in the various regimes will require bigger systems, beyond the reach of our current numerical codes. Nevertheless, we believe that the results for systems of a few thousand grains support the conclusions drawn below. (See Ref  for a discussion of a closely related model.)

In studying force curves for bead packs, several authors have proposed a mean-field argument suggesting that for various granular systems one should observe $`\nu =\alpha +1`$, where $`\alpha `$ is the exponent of the single-contact force law. Thus for our model the mean-field theory predicts $`\nu =2`$. (For the reader’s convenience, we summarize a version of the theory in the Appendix. The treatment there is similar to that of References  and .) The argument is based on the assumption that the rates at which pairs of nearest neighbor sites approach each other may all be taken to be equal to the average rate of compression.

In Subsection C below we present numerical results for our model showing that $`\nu `$ can be significantly less than 2 for some system geometries and in Subsection D we seek to explain the failure of the mean-field argument. In our view, this failure casts doubt on the argument’s applicability to real granular systems. Useful supplementary information is covered in Subsections B and E.

### B Analytic results for limiting cases

It is instructive to consider two cases for which the force curves can be explicitly calculated.

Case 1: $`N_y=1`$. This system consists of only a single layer of random-length springs, as shown in Fig. 7a. Assuming the equilibrium spring lengths to be uniformly distributed in a finite interval, one immediately obtains $`F\mathrm{\Delta }^2`$. In this case the mean-field argument (Appendix A) is exact, as all gaps do close at the same rate.

Case 2: $`N_x=1`$. Because of the periodic boundary conditions, this case is equivalent to a single column of grains without periodic boundary conditions, as shown in Fig. 7b. Again we take the distribution of equilibrium lengths to be a square distribution of width $`2ϵ`$. In this case it is more convenient to fix the force $`F`$ and compute the displacement $`\mathrm{\Delta }(F)`$ since the compressive stress on each grain must be the same. Let us decompose the displacement

$$\mathrm{\Delta }=\underset{n=1}{\overset{N_y}{}}\delta _n,$$
(5)   
where $`\delta _n`$ is the displacement of the $`n`$th spring. When $`F=0`$, we have $`\delta _n=0`$ for each $`n`$, the longer spring at each level being just at the threshold of compression. The growth of $`\delta _n`$ with $`F`$ depends only on $`a_n`$, the difference between the equilibrium lengths of the two springs in the $`n`$th layer. Specifically,

$$\frac{d\delta _n}{dF}=\{\begin{array}{cc}\frac{1}{k}\hfill & \text{if }0<F<ka_n\hfill \\ \frac{1}{2k}\hfill & \text{if }ka_n<F\text{.}\hfill \end{array}$$
(6)   
Note that the $`a_n`$’s are independent and all have the probability density $`P(a)=(2ϵa)/2ϵ^2`$ for $`0<a<2ϵ`$.

The expected value of $`d\mathrm{\Delta }/dF`$ is defined as

$$\frac{d\mathrm{\Delta }}{dF}=\underset{n=1}{\overset{N_y}{}}_0^{2ϵ}\frac{d\delta _n}{dF}P(a_n)𝑑a_n.$$
(7)   
Using the fact that all the terms in the summation are equal and performing the integration, we obtain $`d\mathrm{\Delta }/dF`$, which in turn may be integrated to yield

$$\delta =\{\begin{array}{cc}f\frac{1}{4}f^2+\frac{1}{24}f^3\hfill & \text{if }0<f<2\hfill \\ \frac{1}{2}f+\frac{1}{3}\hfill & \text{if }2<f\text{,}\hfill \end{array}$$
(8)   
where $`\delta \mathrm{\Delta }/N_yϵ`$ and $`fF/kϵ`$ are the nondimensionalized displacement and force per layer, respectively. Eq. (8) shows that, in this case, there is no simple power-law behavior in the large $`N_y`$ limit. Note that Eq. (8) is approximately linear near $`f=0`$ (i.e., $`lim_{f0^+}\left[d(\mathrm{ln}\delta )/d(\mathrm{ln}f)\right]=1`$), and it is linear near and beyond $`f=2`$. If we attempt to identify a single power law for intermediate values of $`f`$, the natural choice is the derivative $`d(\mathrm{ln}\delta )/d(\mathrm{ln}f)`$, evaluated at the point where this quantity is most slowly varying; i.e., where $`d^2(\mathrm{ln}\delta )/d(\mathrm{ln}f)^2`$ vanishes. This method yields $`f\delta ^\nu `$ with $`\nu =1.27\mathrm{}`$, where we have returned to displacement as the independent variable.

The analysis of Case 2 shows that it is possible for the force curve to exhibit a region corresponding to a power less than 2. It also shows that the emergence of a true power law should not be taken for granted in these systems. We find, however, that in sufficiently wide systems, a power law does arise (cf. discussion below).

### C Numerical results

Short, wide systems: Fig.6a shows two typical force curves for an $`80\times 20`$ system. It appears here that there is a regime in which $`\nu 2`$, followed by the expected crossover to $`\nu =1`$ at large $`\mathrm{\Delta }`$. This observation lends some support to the mean field argument and is consistent with the claim of Gilabert et al. , who studied the electrical analogue of our model. (See Section IV E below.)

Tall, narrow systems: In a $`20\times 150`$ system, the exponent $`\nu `$ appears to shift noticeably. The best fit to the curve shown in Fig. 6b is $`\nu =1.7`$ in the intermediate regime of interest. The error in this measurement is estimated to be $`\pm 0.1`$ on the basis of fits made with different choices for which points to exclude from the intermediate regime. The data clearly rule out $`\nu =2`$.

Roughly square system: An estimate of $`\nu =1.9`$ is obtained from data on $`40\times 40`$ systems. Fig. 6c shows data from two $`40\times 40`$ systems. The apparent lack of consistent behavior is presumably due to finite size effects. Nevertheless, using data from 25 runs, an estimate of $`\nu =1.9`$ can be obtained as shown in Fig. 6d. To obtain the dotted line in the figure, the data in the intermediate regime were fit to a power-law; data from the initial linear regime and data from the two largest forces in the figure (where we expect a crossover to linear growth) were excluded in making this estimate for the exponent $`\nu `$.

We suggest that the variation of exponent with aspect ratio can be traced to the tree structure in the system. Specifically, we conjecture that if most of the force is transmitted within a single tree, then a smaller exponent will be observed, whereas if the force is spread over many trees, a larger exponent. We interpret the following numerical experiment as support for this conjecture. Starting with a $`60\times 60`$ grid, we removed all sites lying outside the $`90^{}`$ cone emanating from the center of the bottom boundary. We then measured $`F(\mathrm{\Delta })`$ for this reduced system. In this geometry all the active sites at every stage of compression are connected to the floor at the same point, so all the active bonds are contained within a single tree. As shown in Fig. 6e the exponent $`\nu `$ is close to $`1.5`$, indicating that the behavior within a single tree has a different character than in a system with many trees.

To explore this issue further, we measured the rate at which trees expand as a function of height in the initial configuration. (Roux et al. have investigated a similar, but not identical, model .) While building the initial configuration, we keep track of the root of the tree associated to each site, thereby counting the number $`T(z)`$ of trees that survive up to layer $`z`$. Since only initial configurations are involved — not their subsequent evolution — rather large systems can be simulated. For a system of width 50,000 with data collected for up to $`10^7`$ layers, we obtained an excellent power law fit $`T(z)Az^\gamma `$, with $`\gamma =0.66`$ and $`A=6.5\times 10^4`$, with only a slight deviation for very small heights ($`z<10`$) and very large heights (corresponding to $`T(z)<5`$). Note that $`A=1.3N_x`$. (Incidentally, since $`T(1)=N_x`$, it is clear that the power law fit must be inaccurate for $`z`$ of order $`1`$.) Generalizing from this system to one of arbitrary width, we estimate that in the initial configuration of an $`N_x\times N_y`$ system, approximately $`1.3N_xN_y^\gamma `$ trees will reach the top layer.

Let us speculate on possible consequences of this estimate. The rates of advance of all sites within a single tree are correlated, since compression of the bottom bond of the tree affects the rates of all the branches above it. We therefore expect the sequence of gap closings to depend upon the initial tree structure. Assuming that the system obeys a simple anisotropic scaling law, we conjecture that in large systems the exponent $`\nu `$ is a function of $`N_xN_y^\gamma `$ alone. Limited support for this conjecture comes from two additional runs on $`20\times 80`$ and $`40\times 225`$ systems. In both cases, $`N_xN_y^{0.66}1.1`$, and in both cases we observe power-law behavior with $`\nu 1.8`$, as shown in Fig. 6f. A more rigorous test would require runs on significantly larger systems, for which a more efficient code is needed. Note that the $`40\times 225`$ system was already too large for us to run to completion; the curve in Fig. 6f ends well before the crossover to the linear regime occurs. Note also that in the $`40\times 225`$ system, a clean power-law regime extends over at least two decades in $`\mathrm{\Delta }`$ and three decades in $`F`$, giving us some confidence that a true power-law regime does exist in large systems.

A consequence of our scaling conjecture is that increasing the system size while keeping the aspect ratio fixed should result in exponents that eventually approach that of a wide, short system. As we have shown that $`\nu =2`$ in the $`\mathrm{}\times 1`$ limit (Case 1, Section IVB), we expect to observe the mean-field result $`\nu =2`$ in very large, approximately square systems as well.

### D Why the mean field argument fails

The key assumption in the mean-field argument is that nearest neighbor sites approach each other at the mean compression rate. In terms of the probability density for inter-site distances, the assumption is expressed quantitatively in Eq. (A5) of the appendix. which asserts that the shape and width of the distribution of inter-site distances are independent of the compression $`\mathrm{\Delta }=N_y\delta `$. As discussed in the appendix, in our model the probability density $`P(x,\delta )`$ is well defined only for $`x>0`$; i.e., for bonds whose springs are compressed. However, even restricting our attention to the active bonds, we find that our data are inconsistent with the above assumption. Specifically, in Fig. 8a we show the probability densities for the (nonzero) forces $`f`$, where $`f=\mathrm{max}\{kx,0\}`$, at various stages of compression of a $`40\times 40`$ system, and it is clear that the distribution broadens as $`\delta `$ increases. More quantitatively, in Fig. 8b we plot the widths of the best fits to the data by Gaussian distributions restricted to $`\{x>0\}`$, for several values of $`\delta `$.

The broadening of $`P(x,\delta )`$, which occurs during the early and intermediate regime of compression, is a simple consequence of force balance at branching points of chains of active bonds. Consider a site at which three active bonds meet, forming a “Y” (either right-side-up or upside-down). Because there are no tensile forces, the force in the unpaired branch of the Y is greater than either of those in the paired branches. Moreover, as the Y is further compressed, the rate at which the force in the unpaired branch increases must equal the sum of the rates of increase of the forces in the paired branches. Since the larger force evidently increases more rapidly than the two smaller forces (the forces almost always being increasing functions of $`\mathrm{\Delta }`$), the distribution broadens.

During the late regime of compression, the addition of active bonds will convert Y’s to X’s, for which any correlation between force and bond compression rate must be more subtle. Indeed, when all of the bonds are activated, all bonds compress at exactly the same rate and the force distribution simply shifts uniformly to higher forces – the mean field assumption becomes exact. Thus we expect to observe a broadening of the force distribution during early and intermediate stages of compression, with a rate that decreases as the density of Y’s becomes smaller. This is precisely the behavior displayed in Fig. 8, where the width of the positive $`x`$ portion of the distribution is seen to grow roughly as a power of $`\mathrm{\Delta }`$ less than unity for small $`\mathrm{\Delta }`$ and level off at high compression. A quantitative calculation of the rate of broadening would require a detailed understanding of the statistics of branching in the active network and is beyond the scope of this work.

Let us argue that the broadening of the force distribution generically leads to an exponent $`\nu `$ smaller than the mean field value. Suppose, for definiteness, that $`P(x,\delta )`$ has a sharp leading edge at $`x_1(\delta )`$; i.e., that $`P(x,\delta )=0`$ for $`x>x_1(\delta )`$ and $`lim_{xx_1^{}}P(x,\delta )`$ is bounded away from zero. The existence of Y’s in the active bond network leads to an advance of the leading edge $`x_1(\delta )`$ that is more rapid than the advance of $`\delta `$ itself. Taking into account the formation of new contacts and branch points during compression, and assuming simple asymptotic behavior, we expect that $`x_1(\delta )\delta ^\beta `$ with $`\beta <1`$.  Thus, to lowest order in $`\delta `$, we find $`F=_0^{x1}P(x,\delta )kx𝑑x\delta ^{2\beta }`$, and therefore $`\nu 2`$. Corrections to this exponent should become significant for $`\delta `$ of order unity, which is also the order of the width of $`P(x)`$. This expectation is consistent with Fig. 6, where it can be seen that the crossover to the linear regime for large $`\mathrm{\Delta }`$ occurs for $`\mathrm{\Delta }N_y/5`$, which is the last half decade in the plots.

These qualitative results do not depend on the linearity (under compression) of the springs in our model, nor do they depend on the dimensionality of the model. They should also apply to the vertical forces in a system with horizontal degrees of freedom, providing only that the creation of new branch points is sufficiently common. Thus the mean-field argument seems problematic for physical bead packs. The analysis shows that new contact formation can lead to exponents smaller than the mean field value, an important point in light of other mechanisms proposed to explain experimental observations of this exponent. 

### E Relation to random resistor networks

The analogy between the elastic properties of a network of linear springs and the electrical properties of the same network of resistors has been exploited in numerous studies . In order to clarify the relation between our model and others, and particularly with random networks near the directed percolation threshold, let us consider the analogy in some detail. We will see that there are good reasons to be skeptical of the applicability of percolation model results to our model, but there is also an intriguing numerical coincidence.

There is a formal identity between our equations for mechanical equilibrium

$`f_{ij}^\pm `$ $`=`$ $`k(u_{i,j}u_{i\pm 1,j}s_{ij}^\pm )`$ (9)   
and the electrical equations

$`I_{ij}^\pm `$ $`=`$ $`{\displaystyle \frac{1}{R}}(V_{i,j}V_{i\pm 1,j}v_{ij}^\pm ).`$ (10)   
In the latter equation, $`I_{ij}^\pm `$ specifies the current between sites that are joined by a resistance $`R`$ in series with a battery generating a potential difference $`v_{ij}^\pm `$. Note that the resistors’ conductance plays the role of the springs’ stiffness. The property that the springs function only under compression can be modeled in the electrical system by the insertion of perfect diodes, all directed “downward”, in series with each resistor. Roux et al. studied this very model , but did not report results for varying aspect ratios or single trees and did not study its relation to the $`q`$-model, which had not yet been introduced.

A system with batteries is substantially more complicated than a simple, randomly diluted resistor-diode network. To relate the randomly diluted network at the directed percolation threshold to our model just beyond the initial linear regime, one must first assume that the network of active springs in our model has the same structure as the current carrying paths in a directed network of resistors placed at random on the lattice. One must also assume a relation between $`\mathrm{\Delta }`$ of our model and the probability $`p`$ that a bond exists in the percolation system. The natural assumption here would be that the probability $`pp_c`$ is proportional to $`\mathrm{\Delta }`$, where $`p_c`$ is the critical value for percolation, since $`\mathrm{\Delta }=0`$ is the point where a single force chain (or current carrying path) first forms and on average the bonds in our model are compressed by an amount proportional to $`\mathrm{\Delta }`$ . In our system, however, the batteries play two roles. First, they generate potentials that affect the current distribution even when all the diodes are forward biased, which would correspond to a trivial, homogeneous state of the simple resistor network. Perhaps more importantly, however, they determine which diodes will be forward biased for a given applied potential difference across the whole network. This dynamical process of selecting the current carrying paths may or may not yield structures well-modeled by the random addition of resistive links in a percolation model.

In spite of the difficulties in establishing a connection, it is interesting to compare our results for the power law obeyed by the stiffness $`dF/d\mathrm{\Delta }`$ to the conductivity exponent obtained from the theory of directed percolation in randomly diluted resistor networks . The conductivity exponent for directed percolation has been calculated both numerically and using renormalization group methods and appears to be approximately $`0.7\pm 0.1`$,  which corresponds to a value of $`1.7\pm 0.1`$ for the exponent $`\nu `$. The fact that this agrees with our measurements on tall, narrow systems deserves further study. Other authors have investigated additional details of the statistics near threshold in our model and argued that the system is closely related to the percolation one. 

## V Statistics of forces on individual grains

We now consider the statistics of forces transmitted by individual springs, exploring, in particular, the relation of our results to the $`q`$-model of Coppersmith et al. . Our model is simpler than the granular packings the $`q`$-model was intended to describe. However, our model would appear to be a better candidate for description by the $`q`$-model than the original granular packings, since we have removed the tensorial stresses from our system but still have contacts whose formation is governed by quenched randomness.

In the $`q`$-model in two dimensions, $`q_{i,j}`$ refers to the fraction of the force on site $`(i,j)`$ from above that is transmitted to its neighbor at $`(i1,j1)`$, the complementary fraction $`1q_{i,j}`$ being transfered to site $`(i+1,j1)`$. One assumes that each $`q_{i,j}`$ is a random variable that (1) is independent of the other $`q`$’s and (2) has the same distribution $`\eta (q)`$ at every site, independent of the force supported by that site.

Analytical studies of the $`q`$-model show that the distribution of forces supported by a single grain has an exponential tail at large forces whenever $`\eta (q)`$ is non-vanishing at $`q=1`$.  By contrast, the force distributions in our model, plotted in Fig. 8, show no evidence of exponential tails. At no stage during the compression does it appear that the $`q`$-model distributions are a good match for the distributions we observe. In particular, consider the fourth curve from the left in Fig. 8, which is made at a compression for which almost all sites are active but there remains an appreciable fraction of inactive bonds. The distribution of $`q_{i,j}`$’s directly measured from our data in this regime (see Fig. 9a) is reasonably represented by

$$\eta (q)=0.1\left[\delta (q)+\delta (1q)\right]+0.80.2\mathrm{cos}(2\pi q),$$
(11)   
for which 20% of the $`q`$’s are either 0 or 1. The $`q`$-model would predict an exponential tail in the single grain force distribution.  Even for the small systems we study here, the exponential tail would be clearly distinguishable from the rapid decay we observe. Fig. 9b compares numerical results from our model and from a simulation of the $`q`$-model with $`\eta (q)`$ given by Eq. (11). For the $`q`$-model simulation, forces on the top row were chosen randomly from a uniform distribution on the interval $`(0,2)`$.

To understand the discrepancy between our results and the predictions of the $`q`$-model, we examine our data as regards assumption (2) above: i.e., that the distribution of $`q`$ at each site is independent of the force supported by that site. Fig. 9c shows the distribution of $`q`$’s obtained from equilibrium configurations of our model corresponding to the same conditions as in Fig. 9a, but separated according to the force supported by the site. Different symbols in the plot indicate different levels of force as described in the figure caption. It is clear that the contribution to $`\eta (q)`$ from larger forces is peaked more strongly about 1/2 and has little weight near $`0`$ or $`1`$, which explains why the $`q`$-model prediction fails for this system. As a quantitative measure of the correlation, we have computed the covariance $`C_{qw}(q_{ij}\frac{1}{2})^2(w_{ij}w_{ij})`$ for our model, where $`w_{ij}f_{ij}^{}+f_{ij}^+`$ and the averages are performed over space for one realization of a $`40\times 40`$ system. (The square in this definition is necessary because left-right symmetry guarantees that a correlation function linear in $`q1/2`$ would vanish.) $`C_{qw}`$ vanishes identically in the $`q`$-model since $`q_{i,j}`$ and $`w_{i,j}`$ are independent in that model. As shown in Fig. 9d, $`C_{qw}`$ develops a significant negative value when the system is under compression, indicating that larger forces are associated with $`q`$’s closer to 1/2. The heavy dot in the figure indicates the point corresponding to the data in parts (b) and (c).

## VI Conclusions

We briefly summarize the results of our study of the scalar model, discuss its generalization to three dimensions, and finally draw two conclusions concerning the implications for real systems or more realistic models.

In the context of a toy model, we have tested arguments that have been applied to stresses in static, non-cohesive granular materials. Further study of the model is needed, especially simulations of larger systems, but already two important facts have been established. (1) Correlations in the stress configuration are responsible for substantial effects, both at the level of single-grain forces and that of macroscopic stresses. (2) Depending on the (properly scaled) aspect ratio, new contact formation may play a decisive role in determining the macroscopic stress-strain relationship, with sufficiently tall systems showing power-law behavior with a nontrivial exponent. Regarding (1), we have identified two important effects: (i) the distribution of forces broadens under compression because of force balance constraints at sites where three active bonds meet; and (ii) larger forces tend to divide more evenly between supporting grains as a result of the dynamical process that determines the structure of the active bond network. Theories that neglect these correlations fail when applied to our model.

One may wonder whether qualitatively new features might appear in a 3D generalization of our model in which the vertical direction is taken to be the 111 direction of a simple cubic lattice. As a preliminary check, we have measured the statistics of trees in the initial configuration and observed their general morphology. The number of surviving trees decays as $`z^{1.2}`$ and the trees remain relatively compact. In other words, the diameter of surviving trees grows as $`z^{0.60}`$ (compared to $`z^{0.66}`$ in two dimensions) and the branches of separate trees do not become heavily entangled with each other. On this basis, we conjecture that the essential physics of the stress-strain relation will not be qualitatively different in three dimensions. In particular, we expect a variation in the exponent $`\nu `$ with the scaled aspect ratio.

We draw two general conclusions. First, it appears that the mechanism for the broadening of the vertical force distribution should be present in the full tensorial problem as well, and therefore should be re-examined as a possible explanation for the observation of smaller exponents than those derived from the mean-field argument. Experiments measuring the dependence of the exponent $`\nu `$ on aspect ratio would be especially interesting. Second, the success of the $`q`$-model in predicting the distribution of single-grain forces in bead packs stands in need of explanation. The $`q`$-model has already been criticized for its neglect of tensorial stresses . One may have expected it to provide a better description of the stresses in a scalar model of the type we have studied, but this is not borne out by our results.

We thank Sue Coppersmith and Long Nguyen for helpful conversations, and also Scott Zoldi for his suggestions. This work was supported by NSF through grants DMS-98-9803305 (DGS and MGS) and DMR-94-12416 (JESS).

## A The mean-field argument predicting $`\nu =2`$

In this appendix, we summarize a mean-field argument that was developed to describe forces during the compression of a bead pack between parallel plates (See, e.g., Refs. .)

We define the overlap between two adjacent beads as $`(d_0d)`$ where $`d`$ is the distance between their centers and $`d_0`$ is the nominal distance at which the beads touch but exert no force. Note that because of the minus sign, the overlap is positive for beads in contact, and negative for beads not in contact.

We introduce the random variable, $`X`$ to be the overlap between two randomly chosen adjacent beads. Thus the sample space for $`X`$ includes both the choice of configuration (spring lengths in our model), and the choice of a pair of adjacent beads. Of course $`X`$ also depends on the displacement $`\mathrm{\Delta }`$ of the floor plate. (As in Section I, we normalize so that nonzero forces start at $`\mathrm{\Delta }=0`$). We rescale the independent variable, defining $`\delta =\mathrm{\Delta }/N_y`$ where the bead pack is $`N_y`$ layers thick. Let $`P(x,\delta )`$ be the probability distribution for $`X`$. The average force per spring, $`\sigma `$, is given by;

$$\sigma (\delta )=_{\mathrm{}}^{\mathrm{}}f(x)P(x,\delta )𝑑x,$$
(A1)   
where $`f(x)=\mathrm{max}\{kx,0\}`$. Since

$$f(x)=0\mathrm{for}x<0,$$
(A2)   
the integral may be restricted to $`(0,\mathrm{})`$. Of course, by our normalizations,

$$P(x,0)=0\mathrm{for}x>0,$$
(A3)   
so $`\sigma (0)=0`$. The rescaled variable $`\delta =\mathrm{\Delta }/N_y`$ equals the average change in the overlap $`X`$ resulting from motion of the floor; thus,

$$_{\mathrm{}}^{\mathrm{}}xP(x,\delta )𝑑x=\delta +\delta _0,$$
(A4)   
where $`\delta _0`$ is the average overlap when $`\delta =0`$.

The first main assumption of the mean-field argument is a considerable strengthening of (A4): One assumes that the shape and width of the probability distribution are independent of $`\delta `$, i.e.,

$$P(x,\delta )=P(x\delta ,0).$$
(A5)   
Substituting (A5) into (A1) and rescaling (A2,A3), we deduce that

$$\sigma =k_\delta ^0(x+\delta )P(x,0)𝑑x,$$
(A6)   
which is equivalent to Eq. (1.2) of Ref. .

The second main assumption of the mean-field argument is that $`P`$ is continuous in both arguments and

$$\underset{x0^{}}{lim}P(x,0)=C,$$
(A7)   
where $`C`$ is a finite positive constant. Combining the two assumptions, we find that

$$\sigma =k_\delta ^0(x+\delta )\left[C+𝒪(x)\right]𝑑x,$$
(A8)   
which immediately yields $`\nu =2`$. Under assumption (A5), the only way to accommodate a value of $`\nu `$ less than 2 is to posit that the limit in Eq. (A7) is infinite.

Note that in our model, the distribution for $`X`$ is not uniquely defined for $`X<0`$ since inactive grains are free to relocate slightly. Despite this ambiguity, in section IVD we are able to test hypothesis (A5) by focusing on positive values of $`X`$.

Figures