A mean field description of jamming in non--cohesive frictionless particulate systems
D. A. Head
A theory for kinetic arrest in isotropic systems of repulsive, radially--interacting particles is presented that predicts exponents for the scaling of various macroscopic quantities near the rigidity transition that are in agreement with simulations, including the non--trivial shear exponent. Both statics and dynamics are treated in a simplified, one--particle level description, and coupled {\em via} the assumption that kinetic arrest occurs on the boundary between mechanically stable and unstable regions of the static parameter diagram. This suggests the arrested states observed in simulations are at (or near) an elastic buckling transition. Some additional numerical evidence to confirm the scaling of microscopic quantities is also provided.
Point_J_Long_v2.tex
%%
%% MMA (Mean-Mode-Approximation) applied to "Point J", i.e. the
%% simple (simplest?) jamming model: T=0, spherical particles
%% interacting via strictly repulsive, short ranged, central force
%% potentials
%%
%% Long version, started DAH 2004/11/03
%%

%% This is version 2, 2005/06/09, altered in response to the referee's comments
%% and a query from Craig Maloney regarding the exponents. Actually pretty slight -
%% the referral was basically positive.

%% Document class : can choose between <preprint>/<twocolumn>,
%% <groupedaddress>/<superscriptaddress> (the latter saves space).
%% <draft> highlights overfull boxes
%% <showpacs> shows the PACS numbers
%% <showkeys> shows keywords.
\documentclass[aps,pre,twocolumn]{revtex4}
\usepackage{graphicx}

\begin{document}

\newcommand{\PE}{{\mathcal V}}
\newcommand{\fexp}{\alpha}
\newcommand{\Kzexp}{c}
\newcommand{\Pzexp}{e}
\newcommand{\Vzexp}{g}
\newcommand{\Sm}{\scriptstyle}
\newcommand{\polyexp}{\gamma}
\newcommand{\polyexptwo}{\nu}


\title{A mean field description of jamming in non--cohesive frictionless particulate systems}

\author{D. A. Head}
\affiliation{Division of Physics and Astronomy, Vrije Universiteit, Amsterdam, The Netherlands.\\
Department of Applied Physics, School of Engineering, University of Tokyo, Japan.}

\date{\today}

\begin{abstract}
A theory for kinetic arrest in isotropic systems of repulsive, radially--interacting particles is presented that predicts exponents for the scaling of various macroscopic quantities near the rigidity transition that are in agreement with simulations, including the non--trivial shear exponent. Both statics and dynamics are treated in a simplified, one--particle level description, and coupled {\em via} the assumption that kinetic arrest occurs on the boundary between mechanically stable and unstable regions of the static parameter diagram. This suggests the arrested states observed in simulations are at (or near) an elastic buckling transition. Some additional numerical evidence to confirm the scaling of microscopic quantities is also provided.
\end{abstract}

\pacs{PACS Numbers: 81.05.Rm  83.80.Iz  64.70.Pf  82.70.-y}

\maketitle

%%
%% BODY OF PAPER
%%
\section{Introduction}
\label{s:intro}


Determining when and how a self--assembled particulate system is capable of supporting macroscopic loads appears to be becoming established as a core paradigm in non--equilibrium condensed matter. For many dissipative athermal systems, such as granular media~\cite{Cates:98,Geng:03,Moukarzel:04,Kasahara:04a,Kasahara:04b,Barker:93,Head:00}, foams~\cite{Rafi:DryFoam,Durian:Foam,Bolton:90} or emulsions~\cite{Mason:95,Zhang:05}, the removal of a heat bath or other energy source results in a monotonic decrease of particle motion until a state of kinetic arrest is reached. The morphology, and hence mechanical response, of this static phase is determined by the dynamics of its preparation~\cite{Cates:98}, demonstrating the central importance of the initial dynamic phase and highlighting the non--equilibrium nature of the problem. It is this aspect that delineates self--assembled systems to those with a predefined geometry, such as the tensegrity structures of great importance to engineering~\cite{Calladine:78,Connelly:96,Guest:03}, or the disordered discrete and continuous models of rigidity percolation~\cite{Sahimi:Book,Jacobs:95,Astrom:00,LatvaKokko:01,Head:03,Bergman:84,Feng:84,Arbabi:88,Tang:88,Zhou:03}.

A generic feature of these distinct but related problems is the existence of a {\em rigidity} transition between states that have non--zero elastic moduli, and those that do not~\cite{footnote,Donev:04,OHern:04,OHern:03,Aharonov:99,Makse:EMAFail,Makse:05,Donev_condmat}. For bond diluted lattices, the rigidity transition corresponds to a critical bond dilution, or equivalently a critical mean coordination number~$z_{\rm c}$, that is well approximated by the Maxwell constraint counting method~\cite{Calladine:78}. This transition appears to bear some of the hallmarks of a continuous phase transition in equilibrium systems, including a diverging length scale, although other properties such as universality remain unproven~\cite{Head:03}. For particulate systems, the picture is less clear. Central force models, such as frictionless elastic spheres in the limit of small deformation, appear to reach a similar transition~\cite{Kasahara:04b,OHern:03}, also with a diverging length scale~\cite{Silbert:05}, in the limit of infinite interaction stiffness under controlled pressure. If the volume is controlled, a critical volume fraction must be fine tuned. However, the introduction of friction, at least, appears to overshoot the critical $z_{\rm c}$ by at least a few percent when in gravity~\cite{Kasahara:04a,Kasahara:04b}, the exact amount apparently depending on the damping.



% Schematic figure energy/enthalpy minimisation
\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Min_PV.ps}
\caption{Combined dynamic flow and static stability diagrams in $(z,\omega)$ space, with $z$ the mean coordination number and $\omega>0$ a quantity proportional to the mean particle overlap. Solid (dashed) lines correspond to lines of constant volume $V$ (pressure $P$), with the direction of minimising internal energy $U$ (enthalpy $H$) given by the arrows. Lines to the lower left corner correspond to smaller pressure or larger volume. The system becomes kinetically arrested at the boundary of the (shaded) stable region. The transition point at $(z_{\rm c},\omega)$ is shown as a solid disc (see text for discussion of $z_{\rm c}$).}
\label{f:Min_PV}
\end{figure}


Even when the transition can be approached arbitrarily closely, there remain crucial differences between disordered lattices and particulate systems. In particulate systems, the morphology of the solid state depends on the dynamics of the immediate precursor to arrest, and may therefore include non--trivial, correlated structures. Most studies of lattice models assume morhpologies with no structure beyond the one--bond level. Furthermore particles may become arrested in a state with a finite pressure $P>0$, quite unlike the stressless configurations typically adopted in lattice models.  Attempts have been made to circumvent some of the deficiencies of lattice models by postulating initial configuration generating algorithms that are hoped to mimic the dynamic self--assembly of particulate or atomic systems. Indeed one such scheme, known as bootstrap percolation, has already been applied to this problem~\cite{Schwarz:05}, and by infinite dimensional calculations claim to determine the exponents observed in simulations of central force systems.

The purpose of this paper is to describe an analytical scheme that, when applied to isotropic systems of particles interacting {\em via} central forces, derives the exponents relating pressure, shear modulus and volume fraction to $z-z_{\rm c}$ observed in simulations. There is no mapping to a known model, or assumption of any statistical mechanical analogy, as recently suggested by the elegant theory of  Henkes {\em et al.}~\cite{Henkes:05}. Instead, the properties of the static system are approximated according to a minimum assumption philosophy that proposes a series of intuitive approximations to locally close the equations. For the dynamics, a one--contact level closure is employed, whereas for the statics the approximation takes place somewhere between the one--contact and one--coordination shell level, in that it considers all contacts acting on a particle but ignores correlations between them. The static aspect of this scheme  has already been applied to mixed tensile--compressive systems~\cite{Head:04}.

The central result of this paper is schematically represented in Fig.~\ref{f:Min_PV}, in which the statics and dynamics of the theory have been overlayed into a single plot. The axes are $z$ and $\omega$, where $\omega$ is proportional to the mean particle overlap, so higher $\omega$ corresponds to increasingly compressive contacts. The rigidity transition lies at $z=z_{\rm c}$, $\omega=0$, which corresponds to the same unstressed transition observed in {\em e.g.} lattice models. The anomalous low $z_{\rm c}=3$ (which should be closer to $2d=6$ in this $d=3$ dimension example) is explained later, but is essentially due to the reduction of degrees of freedom inherent in the approximation scheme. The shaded region corresponds to systems that are mechanically stable in that they generate a positive restoring force in response to a localised, linear perturbation. The upper boundary of this stable region, which represents some form of buckling, is locally quadratic, and it is this shape that ultimately determines the exponents.

This static picture is augmented by a global description of the system dynamics from excited to arrested states. Again approximations are required, in this instance the kinetic energy is ignored and the system is assumed to evolve towards state with a lower internal energy (in the case of fixed volume) or enthalpy (in the case of fixed pressure), not unlike thermal systems~\cite{Weiner:Book} but without temperature/entropy terms. Some choice is required in defining volume in terms of $z$ and~$\omega$, but for rapid quenches from highly excited states, a single--particle volume function should be valid, and will robustly give the direction given in the figure. This drives the system towards the stable region. Kinetic arrest is assumed to take place when the boundary of the stable region has been reached; that is, the configurations observed in simulations lie on a line of buckling transitions. By suitably controlling pressure or volume, it is possible to bring the system to the same rigidity transition as in stressless systems, but along a different line to the $\omega\equiv0$ systems considered in disordered lattices. 



%%
%% STATICS: MMA
%%

% Initial E and assumptions in deriving it
\section{Statics: The mean mode approximation}
\label{s:mma}


A variety of approximate analytical treatments for predicting the mechanical response of athermal materials have been devised; just two will be mentioned here by way of comparison. Perhaps the most straightforward is to assume the induced deformation field is affine, so the interparticle displacements are just scaled--down versions of the macroscopic field. This has been applied to frictional granular media~\cite{Walton:87}, but is incapable of predicting a rigidity transition at finite density or volume fraction. Even side--stepping this deficiency by directly comparing elastic moduli to pressure does not resolve its failings~\cite{Makse:EMAFail}, indicating that modes near the transition are inherently non--affine. Enhanced versions have been proposed in which additional degrees of freedom are introduced and determined by variational principles~\cite{Kruyt:98,Trentadue:01}, but only at the cost of significant additional complexity and still not suitable for studying the transition.

Another approach, extensively applied to unstressed systems of predefined morphology, is known as EMA or the {\em effective medium approximation}  (although this term is sometimes also applied to the affine approximation~\cite{Walton:87}). Here the disordered system is replaced by a homogenous analogue with effective interaction parameters. The procedure can be crudely summarised as follows: a single disordered element is inserted into the homogenous bulk and the response to this isolated defect determined using the Green's function. The effective parameters are determined by demanding that the mean response averaged over all possible states of the inserted element is zero. Although this has been successfully applied to a range of continuous and lattice systems~\cite{Sahimi:Book,Feng:85,Schwartz:85,Astrom:00}, it has a fatal shortcoming when applied to particulate systems: the calculation cannot proceed without first identifying an analogous homogeneous material with a known Green's function. It is difficult to postulate such an analogue for finite--size particles, except in certain special cases such as quasi--ordered systems~\cite{Velicky:02}. Furthermore we have already argued that prestress is likely to be important, and even for lattice spring networks the Green's function with prestresses is not known exactly~\cite{Tang:88}. For granular media the problem is particularly acute, as the Green's response has been the subject of substantial debate in recent years (see {\em e.g.}~\cite{Cates:98,Geng:03,Moukarzel:04} and refs therein). The approximation scheme detailed below circumvents these issues, as we describe below.



%The materials under consideration here can all be described as multiscale in that they exhibit qualitatively different structural paradigms depending on what length scale they are viewed, from atomic lengths to the dimensions of the entire macrosopic system. It is therefore important to decide on a suitable level of description before embarking on any analytical procedures. In this work we shall take as our fundamental building block discrete mesoscopic entities of size $\sim10\mu{\rm m}-100\mu{\rm m}$ or larger, namely gas bubbles in wet foams, solid beads in granular media, and immiscible droplets in emulsions. These entities will be hereafter simply referred to as `particles.' The significant mass of these particles means that, to good approximation, they are not susceptible to thermal fluctuations, and so such systems can be treated as athermal.

%At some suitable level of description, all of the above--mentioned particles interact {\em via} short ranged interactions that essentially extend no further than the physical extent of the particles themselves. Furthermore, these interactions are dissipative: for non--zero relative velocities, the total kinetic energy of an interacting subsystem is reduced by viscous damping or other inelastic effects. This energy is transferred to heat which, according the the arguments given above, has no effect on the dynamics of the particles and is therefore effectively lost from the system. Thus without any energy input to the system, such as some form of mechanical noise or `shaking,' the total energy will monotonically decrease until the kinetic energy vanishes simultaneously with the system configuration falling into a local minimum of its potential energy landscape.

For this first exposition of the theory, the non--damping part of the interparticle interaction is assumed to be a purely radial pair potential, generating central forces acting along lines connecting particle centres. This choice is made to reduce the parameter space to a manageable size and to allow a systematic and lucid unfolding of the resulting phase structure of the system. Nonetheless this simplification should hold approximately true for emulsions or wet foams near the transition, when all particles are only slightly deformed, although it will clearly not apply to strongly deformed particles (or dry foams~\cite{Rafi:DryFoam}). For real granular media there will also arise transverse surface forces mediated by non--zero friction at the bead--bead interface; purely central forces correspond to frictionless spherical beads, which only truly exist in numerical simulations.
%Nonetheless we shall continue to refer to simulations of frictionless granular models, in the understanding that a comparison to real grains must await an enhanced version of the theory.

\subsection{Determining the mechanical stability}

Given central force interactions, there is no need to track particle orientations and the static system is fully specified by the $d$--dimensional position vectors ${\bf x}^{\beta}$ of all particles $\beta$. The force on $\beta$ due to $\gamma$ is denoted by~${\bf f}^{\gamma\beta}$,

\begin{equation}
{\bf f}^{\gamma\beta}
%=
%{\bf f}({\bf x}^{\beta}-{\bf x}^{\alpha})
=
f(r^{\gamma\beta}){\hat {\bf n}}^{\gamma\beta}
\end{equation}

\noindent{}where $r^{\gamma\beta}=\mid{\bf x}^{\beta}-{\bf x}^{\gamma}\mid$ is the distance between particle centres and ${\hat {\bf n}}^{\gamma\beta}=({\bf x}^{\beta}-{\bf x}^{\gamma})/r^{\gamma\beta}$ is the unit vector from $\gamma$ to $\beta$. In this context the scalar central force $f(r)$ is usually taken to be of the form

\begin{equation}
f(r)
=
\left\{
\begin{array}{lll}
\mu\left(1-\frac{r}{r_{0}}\right)^{\fexp}
& : & r<r_{0} \\
0 & : & {\rm otherwise}
\end{array}
\right.
\label{e:f_r}
\end{equation}

\noindent{}with $\fexp=1$ (truncated Hookean) or $\fexp=3/2$ (Hertzian), and $r_{0}$ is the sum of the two particle radii. The prefactor $\mu>0$ is typically treated as a particle--independent parameter, although strictly speaking it is a function of the radii of curvature at the point of contact for Hertzian interactions~\cite{LandauLifshitz,Makse:05}. Note that with this sign convention, positive $f$ corresponds to compressive forces and negative $f$ to tensile ones.

Suppose we are given an initial configuration $\{{\bf x}^{\beta}\}$ that is static, {\em i.e.} the vector sum of all contact forces on each particle vanishes. To determine its linear stability, apply an arbitrarily small external force $\delta{\bf f}^{\rm ext}$ onto the particle lying nearest some arbitrary point in space --- call this particle~$\alpha$ (not to be confused with the force law exponent). If the system is stable, it will move to a nearby static configuration in which all particles $\beta$ have been displaced to ${\bf x}^{\beta}+\delta{\bf x}^{\beta}$ with $|\delta{\bf x}|\ll r_{0}$. Force balance must again be obeyed, {\em i.e.} the changes in contact forces on $\beta$ sum to zero for $\beta\neq\alpha$, and to $-\delta{\bf f}^{\rm ext}$ for particle~$\alpha$.

The response $\delta{\bf x}^{\alpha}$ for a particular configuration $\{{\bf x}^{\beta}\}$, even if tractable, would be of no practical interest and we must instead ensemble average, keeping fixed a set of parameters that dominate the mechanical response of the system. A wealth of data has shown the mean coordination number $z$ to be a crucial factor in determining stability, and Alexander~\cite{AlexanderRev} has highlighted the importance of prestresses, so we also assume that both $z$ and some measure of the initial contact force distribution is kept fixed.  Once the ensemble (denoted by the angled brackets $\langle\ldots\rangle$ below) has been suitably defined, the requirement of force balance on the perturbed bead $\alpha$ can be written

\begin{equation}
\delta {\bf f}^{\rm ext}
-
\left\langle
\sum_{\beta\sim\alpha}
\delta{\bf f}^{\alpha\beta}
\right\rangle
=
0
\label{e:force_balance}
\end{equation}

\noindent{}where the sum is over all $\beta$ interacting with~$\alpha$. The change in contact force $\delta{\bf f}^{\alpha\beta}$ can be related to the particle displacements $\delta{\bf x}^{\alpha}$ and $\delta{\bf x}^{\beta}$ by

\begin{equation}
\delta f^{\alpha\beta}_{i}
=
A^{\alpha\beta}_{ij}
\left(
\delta x^{\beta}_{j}
-
\delta x^{\alpha}_{j}
\right)
\label{e:delta_f_ab}
\end{equation}

\noindent{}with summation over Roman indices only. The $d\times d$ matrix $A^{\alpha\beta}$ is defined by~\cite{Tanguy:02}

\begin{equation}
A_{ij}=
\frac{f(r)}{r}
(\delta_{ij}-\hat{n}_{i}\hat{n}_{j})
+f^{\prime}(r)
\hat{n}_{i}\hat{n}_{j}
\label{e:A_ij}
\end{equation}

\noindent{}assuming $f(r)$ is continuous with a finite first derivative $f^{\prime}(r)\neq0$ over all $r$ of interest. Here and below the suffices $\alpha$, $\beta$ are dropped whenever the meaning is clear.



% Actual derivation of the MMA
\subsection{Derivation of the MMA}

%Each $\delta{\bf f}^{\alpha\beta}$ can be written in terms of the displacements $\delta{\bf x}^{\alpha}$ and $\delta{\bf x}^{\beta}$ using the $A^{\alpha\beta}_{ij}$ matrix defined earlier~(\ref{e:A_ij}).


So far this is exact but intractable. Now we approximate. Consider that, for an isotropic system, the perturbed bead must move parallel to the external force after averaging, $\langle\delta{\bf x}^{\alpha}\rangle=\lambda\delta{\bf f}^{\rm ext}$ with an unknown compliance~$\lambda$. The philosophy of the MMA is to impose this form {\em before} averaging, {\em i.e.} inside the brackets in~(\ref{e:force_balance}). In this way the dependency of $\delta{\bf x}^{\alpha}$ on the entire initial configuration $\{ {\bf x}^{\gamma}\}$ is subsumed into the single scalar parameter~$\lambda$. This is clearly a significant saving in terms of complexity, although it disallows the transverse motion of the particle and hence reduces the degrees of freedom; the consequence of this on the location of the rigidity transition will be discussed later.

The logical continuation of this approach is to similarly replace the $\delta{\bf x}^{\beta}$ by $\langle\delta{\bf x}^{\beta}\rangle$; however, this averaged form cannot be determined by symmetry considerations alone. Instead we assume here that the change in the contact force with $\alpha$ can be treated as an external force on~$\beta$, so that $\delta{\bf x}^{\beta}=\lambda\delta{\bf f}^{\alpha\beta}$ with the same $\lambda$ as before. Intuitively, this corresponds to the statement that the displacement of $\beta$ is dominated by the change in contact force with $\alpha$, which, for a monotonically decaying force field extending outwards from~$\alpha$, should at least not be embarrassingly wrong. These two approximations taken together allows each contact force $\delta{\bf f}^{\alpha\beta}$ to be uniquely determined from $\delta{\bf f}^{\rm ext}$, as found by inserting $\delta{\bf x}^{\alpha}=\lambda\delta{\bf f}^{\rm ext}$ and $\delta{\bf x}^{\beta}=\lambda\delta{\bf f}^{\alpha\beta}$ into (\ref{e:delta_f_ab}) and (\ref{e:A_ij}) and inverting,

\begin{eqnarray}
%\begin{equation}
\delta f^{\alpha\beta}_{i}
&=&
S^{\alpha\beta}_{ij}
\delta f^{\rm ext}_{j},
\nonumber\\
%\end{equation}
%\noindent{}with
%\begin{eqnarray}
S^{\alpha\beta}_{ij}
&=&
\left[
1+(\lambda|f^{\prime}(r^{\alpha\beta})|)^{-1}
\right]^{-1}
{\hat n}^{\alpha\beta}_{i}{\hat n}^{\alpha\beta}_{j}
\nonumber\\
&&+
\left[
1-\left(
\frac{\lambda f(r^{\alpha\beta})}
{r^{\alpha\beta}}
\right)^{-1}
\right]^{-1}
\left(\delta_{ij}-{\hat n}^{\alpha\beta}_{i}{\hat n}^{\alpha\beta}_{j}\right)
\nonumber\\
\label{e:Sij}
\end{eqnarray}

\noindent{}Thus each $\delta{\bf f}^{\alpha\beta}$ is independent of the others. Note that the unphysical singularity at $\lambda f/r=1$ is avoided by the stability equation below.

It is apparent from (\ref{e:Sij}) that the MMA has reduced the global problem to a local one, in which the response of each contact $\delta{\bf f}^{\alpha\beta}$ depends only on the interparticle separation $r^{\alpha\beta}$ (through which $f(r)$ and $f^{\prime}(r)$ are found), the unit vector $\hat{\bf n}^{\alpha\beta}$, and the coordination number for this bead $z^{\alpha}$, implicit in the summation~(\ref{e:force_balance}). Before the averaging can be completed, it is necessary to specify how these quantities vary. Here we will deliberately take simple forms to facilitate transparent interpretation of the results. Firstly, $z^{\alpha}$ is taken to be independent of the contact forces and orientations, allowing the $z^{\alpha}$--averaging to be performed and the the force balance equation (\ref{e:force_balance}) rewritten as

\begin{equation}
\delta f^{\rm ext}_{j}
\left\{
\delta_{ij}
-
z
\left\langle
S^{\alpha\beta}_{ij}
\right\rangle_{\alpha\beta}
\right\}
=
0
\label{e:forcebal_sep_z}
\end{equation}

\noindent{}where $z$ is the mean coordination number, and the averaging is now over ${\bf n}^{\alpha\beta}$ and $r^{\alpha\beta}$. Since $\delta{\bf f}^{\rm ext}$ is arbitrary, the quantity inside the brackets in (\ref{e:forcebal_sep_z}) must vanish.

We further assume that the bond orientations are independent of the contact forces. The $\hat{\bf n}^{\alpha\beta}$ are taken to be independent, identically distributed random variables, uniformly distributed over the $d-1$ dimensional unit hypersphere. This neglects any correlations in the topology of the network. It also ignores excluded volume effects, as it allows the same particle to have bonds arbitrarily close together (so the contacting particles would significantly overlap); this should not be crucial near the transition, which is determined by stability requirements rather than excluded volume, but will be relevant at higher~$z$. Performing the average gives

\begin{eqnarray}
z
\Bigg\langle
\frac{d-1}{d}
\left[
1-
\left(
\frac{\lambda f}{r}
\right)^{-1}
\right]^{-1}
&
+
&
\frac{1}{d}
\left[
1-(\lambda f^{\prime})^{-1}
\right]^{-1}
\Bigg\rangle
\nonumber\\
&=&
1
\label{e:forcebal_z_n}
\end{eqnarray}

\noindent{}where the identity $\langle \hat{n}_{i}\hat{n}_{j}\rangle=\frac{1}{d}\delta_{ij}$ has been used.

It remains to specify how the contact forces are distributed. In principle only configurations consistent with force balance should be allowed, but this complication becomes redundant given the approximations leading to~(\ref{e:Sij}), which allows the response from each contact to be calculated independently. A natural choice is then to assume that each contact force $f(r)$, or equivalently each interparticle separation~$r$, are identically and independently distributed according to some given distribution. Clearly this neglects any correlations in the initial force network, but force balance is ensured on average by virtue of the uniform distribution of contact angles already employed. Some calculations for general force distributions will be described later. For now, the simplest choice possible is made, namely a delta--function distribution corresponding to a monodisperse separation $r$, force $f(r)$ and gradient $f^{\prime}(r)\neq0$ for each contact. It is then possible to insert (\ref{e:Sij}) into (\ref{e:forcebal_z_n}) and integrate; the result is finally

\begin{equation}
d\left(\frac{1}{z}-1\right)
=
(d-1)
\frac{1}{\frac{\displaystyle\lambda f(r)}{\displaystyle r} - 1}
-
\frac{1}{1+\lambda\mid f^{\prime}\mid}
\label{e:mma_lambda}
\end{equation}

\noindent{}This is a quadratic equation in~$\lambda$, which is easily solved.

Although $\lambda$ is in principle a measurable quantity, it is more useful to specify results in terms of the shear modulus $G$ or bulk modulus $K$. These can be related to the compliance $\lambda$ by specifying some suitable closed surface $S$ and applying an external force to each particle it cuts, {\em i.e.} forces parallel to $S$ for $G$ and normal for~$K$. Given the displacement of each particle is $\propto\lambda^{-1}$ according to the MMA, it is straightforward to see that

\begin{equation}
G\sim B_{\rm G}\lambda^{-1},
\quad
K\sim B_{\rm K}\lambda^{-1},
\label{e:KGDef}
\end{equation}

\noindent{}where the prefactors $B$ have dimension (length)$^{2-d}$. The relevant length scales are the characteristic length of the enclosing surface $L\sim\sqrt[d-1]{S}$ and the particle radius $r_{0}/2$, but the local closure of the MMA equations means we are unable to determine their weighting in $B_{\rm G}$ and~$B_{\rm K}$. However, if $\lambda$ diverges with~$L$, as is trivially true for a $d=1$ system with fixed boundaries, then $B\sim Lr_{0}^{1-d}$ to ensure finite moduli for arbitrarily large systems. Some form of divergence of $\lambda$ with $L$ is also expected from linear continuum elasticity~\cite{LandauLifshitz}. This subtlety is sidestepped below, where $L$ is assumed to be fixed and finite.






%Here we will use the bulk modulus~$K$, which can be related to $\lambda$ as follows: apply an external force $\delta{\bf f}^{\rm ext}$ to every particle that cuts the surface $S$ of a macroscopic region of the system with volume~$V$. Each $\delta{\bf f}^{\rm ext}$ is of equal magnitude but perpendicular to~$S$, so the change in mean isotropic stress inside $V$ is $\delta P=\delta f^{\rm ext}N_{\rm S}/S$ with $N_{\rm S}$ the number of perturbed particles. By the MMA, the surface will contract by $\lambda\delta f^{\rm ext}$ at every point, so the volume will change by $\delta V=-\lambda S\delta f^{\rm ext}$ and

%\begin{equation}
%K
%=
%-
%V\frac{\delta P}{\delta V}
%=
%B
%\frac{L}{\lambda}
%\label{e:K_def}
%\end{equation}

%\noindent{}where $L=V/S$ is the characteristic length of the region and $B^{-1}=S/N_{\rm S}$ is the portion of the surface area per bead, so $B\sim r_{0}^{1-d}$. A negative $K$ denotes an unstable system. Furthermore, a complex $\lambda$ means that the assumption that the nearby perturbed configuration obeys force balance is unachievable, also signifying instability. 

%For $K$ to take an $L$--independent form, it must be true that $\lambda\sim L$, {\em i.e.} $\lambda$ diverges with system size, even far from the transition point. This should be expected: in linear continuum elasticity, the displacement of a local region perturbed by a point force also diverges. For one dimensional systems bounded by fixed walls, it is also clear that $\lambda$ should scale with $L$ (since each bead will contract by the same amount). The local closure of the MMA equations means that it cannot be used to directly probe finite size effects to check this dependency, and here it is simply assumed that it holds true. Instead we focus on the boundaries between regions where $\lambda$ is positive, negative or complex, which remain fixed under the scaling of $\lambda$ by~$L$. We start with the special case $d=1$, and then consider $d\geq2$ separately.


% d=1 : trivial 
%\subsubsection{Dimension d=1}
%For $d=1$, (\ref{e:mma_lambda}) reduces to $\lambda | f^{\prime} |=(z-1)^{-1}$ and hence $G,K\sim|f^{\prime}|(z-1)$, suggesting a rigidity transition at $z_{\rm c}=1$. In fact rigidity should occur at $z_{\rm c}=2$, a discrepancy that is discussed below for $d\geq2$. More serious is that the MMA equations completely fail to predict the order of the transition in this dimension. In a real $d=1$ system, $K$ is zero for all $z<2$, but then jumps to a value proportional to $|f^{\prime}|$ at $z=2$ when a trivial spanning cluster emerges. By contrast, the MMA predicts a continuous transition in all dimensions. We conclude that the MMA cannot be applied to $d=1$, but since this is an anomalous case (and can in any case can be analytically solved), this should not be viewed as a significant failing and is mentioned here principally for completeness.


% Schematic figure of stable and unstable regions - no constant P/V lines
\begin{figure}
\centering
\includegraphics[width=\columnwidth]{CompTens.ps}
\caption{Stable regimes of $(z,\omega)$ space in the MMA model, for $d=3$. The black disc at $(d,0)$ is the (unstressed) rigidity percolation transition. All points $z<d$, $\omega=0$ are unstable. For $\omega<0$, corresponding to tensile bonds, the system is always stable (light grey region). For $\omega>0$, corresponding to compression, only the dark grey shaded region is stable.}
\label{f:CompTens}
\end{figure}


% d>1 : Interesting-er.
%\subsubsection{Dimensions $d\geq2$}
\subsection{Stability regimes}

Ignoring the trivial $d=1$ case, dimensionality only enters {\em via} the prefactors and so all $d\geq2$ will be discussed together. The solutions to (\ref{e:mma_lambda}) can be conveniently expressed in terms of the two scalars $z$ and $\omega$, where

\begin{eqnarray}
\omega
&=&
\frac{f(r)/r}{\mid f^{\prime}(r)\mid}
\\
&\sim&
\frac{1}{\fexp}
\left(
1-\frac{r}{r_{0}}
\right)
\label{e:omega}
\end{eqnarray}

\noindent{}is a dimensionless measure of the prestress in the system. The second form (\ref{e:omega}) holds for the particulate potentials (\ref{e:f_r}) in the limit $r\rightarrow r_{0}^{-}$, {\em i.e.} close to the rigidity transition, which is the regime of interest here. A schematic description of the predictions of the MMA are given in Fig.~\ref{f:CompTens}. Some important features are now discussed.

\vspace{\baselineskip}

\noindent{\em $\omega\equiv0$:}
This is the unstressed case in which all contact forces are initially zero, nullifying use of the particulate potentials~(\ref{e:f_r}), which have no linear response at $r=r_{0}$, but still attainable for non--truncated Hookean springs. The equation for $\lambda$ gives the single solution $\lambda|f^{\prime}|=(z/d-1)^{-1}$, or $G,K\sim(z/z_{\rm c}-1)^{f}$ with the transition point $z_{\rm c}=d$ and an exponent $f=1$. The effective medium theory for diluted spring lattices also predicts $f=1$, but at the higher transition point $z_{\rm c}=2d$ in accord with the Maxwell counting estimate~\cite{Feng:85}. The transition value found here, $z_{\rm c}=d$, seems anomalous until one recalls that the basic assumptions of the MMA restrict the motion of the particles to mean forms, thereby reducing their degrees of freedom and hence lowering~$z_{\rm c}$. Despite this, the MMA still predicts a finite transition, and can therefore can be used qualitatively. Any unease over the actual value could be lessened by referring to it as an effective coordination number~$z^{\rm eff}$ if desired.

\vspace{\baselineskip}

\noindent{\em $\omega<0$:}
When all of the bonds are tensile, there is always one real, positive solution of $\lambda$ extending from arbitrarily large $z$ down to a lower value $z_{\rm min}=1$. Again this value is too small; $z_{\rm min}=2$ is more probable, {\em i.e.} infinite chains of particles spanning the system. As $\omega\rightarrow0^{-}$ with $z>z_{\rm c}$ fixed, the single root of $\lambda$ continuously approaches the unstressed solution given above. Repeating this procedure for $z<z_{\rm c}$, however, reveals that $\lambda$ diverges as $|\omega|^{-1}$ and hence $G,K\sim|\omega|$ vanishes continuously as the unstressed axis is reached. Thus just below the $\omega=0$ line, the elastic moduli are very small and the system is inherently weak, becoming weaker as $z$ decreases. This may explain why the few attempts to survey this region in disordered lattices~\cite{Tang:88,Zhou:03} have observed a rapid but gradual crossover of the transition from $z_{\rm c}$ to $z_{\rm min}$: numerical noise and/or arithmetic precision may incorrectly attribute zero values to small but finite moduli.

\vspace{\baselineskip}

\noindent{\em $\omega>0$:}
For compressed bonds, the $(z,\omega)$ plane is partitioned into a stable region with two distinct real, positive roots, and an unstable region for which both roots are either complex or negative. The boundary between the stable and unstable regions is quadratic near $z_{\rm c}$,

\begin{equation}
\omega_{\rm bdy}
\approx
\frac
{(z-z_{\rm c})^{2}}
{4d^{2}(d-1)},
\quad\quad z>z_{\rm c}.
\label{e:zw_boundary}
\end{equation}

\noindent{}Both roots of $\lambda$ coincide on the boundary,

\begin{equation}
\lambda_{\rm bdy}
\approx
\frac{z-z_{\rm c}}{2d(z-1)}
\frac{r_{0}}{\mu}
\left(
\frac{4d^{2}(d-1)}{\fexp(z-z_{\rm c})^{2}}
\right)^{\fexp}
\label{e:lambda_boundary}
\end{equation}

\noindent{}and hence $G_{\rm bdy},K_{\rm bdy}\sim\lambda_{\rm bdy}^{-1}\sim(z-z_{\rm c})^{2\fexp-1}$. Starting from the stable regime and decreasing $\omega$ to zero, one of the roots diverges as $\omega^{-1}$ while the other continuously approaches the unstressed solution, crossing over to become the single root in the tensile regime (where the other root becomes negative).

The manner in which the compressive system becomes unstable is noteworthy. On the boundary, $r_{0}-r\sim(z-z_{\rm c})^{2}$, $f\sim(z-z_{\rm c})^{2\fexp}$ and $f^{\prime}\sim(z-z_{\rm c})^{2(\fexp-1)}$, which according the (\ref{e:Sij}) means that the force transfer is predominantly longitudinal. As already noted by Alexander~\cite{AlexanderRev}, in such cases the change in energy will be positive, from which he infers the system should be stable. However, there are other ways of buckling. An established alternative is a bifurcation to a different class of solution~\cite{Libai&Simmonds}; we might also speculate that the energy landscape may exhibit discontinuities in the limit of infinite system size, allowing some form of catastrophic buckling. In fact, the buckling as envisaged by Alexander, which corresponds to $\lambda<0$ here, {\em does} arise within the MMA, but only for $z<z_{\rm c}$ and small $\omega>0$. The upper boundary in Fig.~\ref{f:CompTens} rather corresponds to when $\lambda$ becomes complex.

%%
%% ENERGY/ ENTHALPY MINIMISATION
%%

\section{Dynamics: Energy minimisation}
\label{s:min}

A system not in a shaded region in Fig.~\ref{f:CompTens} will destabilise under any non--zero noise, evolving its contact network according to the dynamical particle interactions and hence allowing $z$ and~$\omega$ to vary. It can only come to rest in a mechanically stable region. Indeed for sufficiently damped interactions, as assumed here, kinetic arrest can be identified with the point at which the system first touches a stable region. Strong damping also means that the kinetic energy is always small, so that the system will evolve to minimise some suitably defined energy potential. For constant volume $V$, this potential is the internal energy $U$, here just the total potential energy stored in the interparticle bonds. For controlled pressure, the corresponding potential is the enthalpy $H=U+PV$~\cite{Weiner:Book}.

A crucial problem in integrating the fixed $P$ or $V$ dynamics and the $(z,\omega)$ stability diagram is writing down expressions relating $V$ to $z$ and~$\omega$. This is likely to be a subtle issue; under gentle shaking, the particles may form spatially extended structures that would necessitate many--particle variables to calculate~$V$~\cite{Edwards:05,Barker:93,Head:00}, which is clearly beyond the one--particle closure of the MMA equations. For now we ignore such potential pitfalls, and instead assume the following, one--particle description, in the expectation that it will hold in the initial dynamic phase from a highly excited initial state. Simply assume that $V$ is a decreasing function of both $z$ and~$\omega$, as might be expected for uniform, global changes of these variables. The precise choice of $V(z,\omega)$ should incorporate the large changes in $z$ that are possible for small changes in $r$ when the particles are barely touching, {\em i.e.} when $\omega\ll1$. This can be written as

\begin{equation}
\frac{V,_{z}}{V,_{\omega}}
=
D\omega^{b}
\quad
{\rm as}
\quad
\omega\rightarrow 0
\label{e:exponent_b}
\end{equation}


\noindent{}where the unknown exponent $b$ is assumed here to obey $b\geq1$ ($b<1$ alters the scaling behaviour of $V$ with $z$ described later, but not of $P$, $G$ or $K$). The dimensionless constant $D>0$ is some material--dependent parameter.
%Then lines of constant $V$ can be found by solving ${\rm d}V={\rm d}\omega+D\omega^{b}{\rm d}z=0$, which reaches the unstressed $\omega=0$ axis only if $b<1$. For $b\geq1$, the lines asymptote to $\omega=0$ but never reach it within finite~$z$.
 
%First, we consider the simple case of systems with a given, immutable list of contacts, so only $\omega$ may evolve with $z$ fixed. This naturally arises in the case of lattice spring models, and so we shall focus attention on Hookean interactions for this case, allowing consideration of compressed, tensed and unstressed systems. The repulsive--only particulate potentials (\ref{e:f_r}) also allow $z$ to evolve and will be considered separately.

% Not allowing the contact number to change
%\subsection{Fixed $z$ and variable $\omega$ for Hookean interactions}

%For $\fexp=1$, the internal energy is just $U=\sum_{\rm bonds}\frac{1}{2}r_{0}\mu\omega^{2}$ which is minimised when all the $\omega\equiv0$, {\em i.e.} every bond is unstressed. This is the usual starting condition for lattice models. However, it may not always be reached: if the system is initially in a state of tension, or a state of compression with $z>z_{\rm c}$, it will have to first pass through a stable region before reaching the unstressed axis $\omega=0$. It will therefore become arrested in some stable state with $\omega\neq0$, where it will remain until non--linearly perturbed or somehow "shaken," when it may again start to flow to a state with a lower~$U$. Once reaching the minimum of~$U$, any perturbations will merely cause the system to fluctuate around $U=0$.

%The constant--$NP$ case is here somewhat trivial. Since $P$ should be a monotonic increasing function of both $z$ and~$\omega$, fixing $P$ with $z$ constant means a fixed $\omega$, {\em i.e.} the choice of pressure determines~$\omega$. Any $P<0$ will give a stable, tensile system, whereas $P\geq 0$ may only be stable if $z>z_{\rm c}$ (and $P$ sufficiently small). If an unstable combination of $z$ and $P$ is chosen, the system is unable to reach a static configuration and will flow as a fluid. For $P=0$ this situation can continue indefinitely, but for $P>0$ any system flow will be accompanied by a reduction in volume, which will bring particles closer together and hence increase~$z$, at which point the assumption of a fixed contact network must be relaxed.


% Allowing the contact number to change, applied to the repulsive-only case
%\subsection{Variable $z$, $\omega$ for repulsive particles}

It is straightforward to derive the internal energy $U$ by employing the same approximations as used to perform the integration of the MMA equations, namely a constant $\omega=\frac{1}{\fexp}(1-r/r_{0})$ and isotropic bond orientations $\hat{\bf n}$,

\begin{equation}
U
=
\frac{Nz}{2}
\frac{r_{0}\mu}{\fexp+1}
(\fexp\omega)^{\fexp+1}
\label{e:U}
\end{equation}

\noindent{}which is the total number of contacts $Nz/2$ multiplied by the bond potential. Similarly, the isotropic pressure $P\delta_{ij}$ is the sum of the $r_{i}f_{j}$ for each bond, divided by the volume~$V$, or (for $r\approx r_{0}$)

\begin{equation}
PV
=
\frac{Nz}{2}
\frac{r_{0}\mu}{d}
(\fexp\omega)^{\fexp}
\label{e:PV}
\end{equation}

\noindent{}where the identity $\langle {\hat n}_{i}{\hat n}_{j}\rangle=\frac{1}{d}\delta_{ij}$ has been used.

Performing the minimisation for both cases reveals broadly the same behaviour; the system will evolve in the direction of increasing $z$ and decreasing $\omega$, as schematised in Fig.~\ref{f:Min_PV} and already discussed in the introduction. For example, for fixed~$V$, the extremum of $U$ (which we assume is the minimum) is found by solving $dU=0$ simultaneously with $dV=0$, where the latter gives the constraint of constant volume. This can be rearranged to give $U,_{\omega}V,_{z}=U,_{z}V,_{\omega}$ and hence from (\ref{e:U}),

\begin{equation}
z\omega^{a}V,_{z}
=
\frac{\omega^{\fexp+1}}{\fexp+1}
V,_{\omega}
\label{e:min_constV}
\end{equation}

\noindent{}This admits the single solution $\omega=0$ in the small--$\omega$ regime of interest here, so $U$ is minimised when all particles are at the limits of their interaction potentials (or at their natural lengths for Hookean springs). Given that the system is constrained to move on lines of constant~$V$~(\ref{e:exponent_b}), the minimum corresponds to a divergent~$z$, clearly unobtainable in a real system. Excluded volume and ordering effects must be incorporated into the theory before any large--$z$ treatment can be attempted. Repeating the enthalpy minimisation at fixed $P$ gives essentially the same behaviour.

%Repeating these calculations with constant $P$ is straightforward. The condition for the enthalpy $H=U+PV$ to be minimised is $H,_{z}P,_{\omega}=H,_{\omega}P,_{z}$, which for small $\omega$ leads to\begin{equation}\omega^{\fexp}P,_{\omega}=za\omega^{\fexp-1}P,_{z}\label{e:min_constP}\end{equation}\noindent{}The derivatives $P,_{\omega}$ and $P,_{z}$ near the line $\omega=0$ can be found from the explicit expression for PV~(\ref{e:PV}) and the assumed asymptotic form for $V$~(\ref{e:exponent_b}). It then becomes clear that there are no solutions for small $\omega$, and the solution for $\omega=0$ when $a>1$ is never reached since lines of constant $P\neq0$ never touch the $\omega=0$ axis. The enthalpy minimum therefore corresponds to the limit $z\rightarrow\infty$, $\omega\rightarrow0^{+}$, which again takes the system into the high--$z$ regime.



\subsection{Kinetic arrest and scaling behaviour}

Given minimisation drives the system in the direction of small $\omega$ and large $z$, they will enter the mechanically stable region of the $(z,\omega)$ diagram Fig.~\ref{f:Min_PV} somewhere along the stability boundary. For overdamped dynamics, we assume it also stops there, allowing the scaling behaviour of various quantities with $z-z_{\rm c}$ to be determined. According to (\ref{e:zw_boundary}), the microscopic variable $\omega\sim(z-z_{\rm c})^{2}$, which is confirmed by $d=3$ numerical simulation of Hertzian spheres shown in Fig.~\ref{f:Z_Omega}. The relation between the elastic modulus $G$ and $z-z_{\rm c}$ on the boundary has already been given~(\ref{e:lambda_boundary}),

\begin{equation}
G
\sim
(z-z_{\rm c})^{\Kzexp}
\label{e:Kzexp}
\end{equation}

\noindent{}with $\Kzexp=2\fexp-1$. Note that this diverges as $\fexp\rightarrow\frac{1}{2}^{+}$, signifying the breakdown of linear response in this admittedly atypical class of pair potentials.
%Interestingly for the Hookean case $\fexp=1$, this predicts the same linear dependence of $K$ as found in the unstressed $\omega=0$ case, although closer inspection reveals that the prefactors are different. 
At the transition point $(z,\omega)=(z_{\rm c},0)$, the pressure $P$ is zero with a finite volume $V_{0}$, so $PV\approx PV_{0}$ to first order and scaling relations involving $P$ can be found independently of the choice of~$V$~(\ref{e:exponent_b}). From (\ref{e:PV}),

\begin{equation}
P
\sim
(z-z_{\rm c})^{\Pzexp}
\label{e:Pzexp}
\end{equation}

\noindent{}with $\Pzexp=2\fexp$. Finally, for relations involving $V$, and hence the volume fraction $\phi\sim V^{-1}$, we find

\begin{equation}
V_{0}-V
\sim
(z-z_{\rm c})^{\Vzexp}
\end{equation}

\noindent{}with $\Vzexp=2$ when $b\geq1$ ($\Vzexp=2b$ for $b<1$). These results are summarised in Table~\ref{t:exponents}, where they are compared to simulation results on various central force systems.

% OCTA data
\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Z_Omega.eps}
\caption{Confirmation of the predicted scaling $\langle\omega\rangle\sim(z-z_{\rm c})^{2}$ for $d=3$ Hertzian spheres relaxed at constant volume. The straight line is a fit to $\omega\propto(z-z_{\rm c})^{h}$ with $h=1.92(6)$ and $z_{\rm c}=5.86(4)$, where figures in brackets denote single standard errors on the last digit and $z_{\rm c}<2d=6$ since (rigidly) disconnected rattler particles were included in the average. The simulations were performed on an $N=1000$ monodisperse particle system using the conjugate gradient minimisation routine in the COGNAC module of the open source numerical suite OCTA~\cite{OCTA}.
}
\label{f:Z_Omega}
\end{figure}

It should be noted that the postulated forms for $V_{0}-V$ and $P$ already allow the bulk modulus~$K$ to be determined without recourse to calculating the mechanical response, and as expected (given the one--bond expressions for $P$ and $V$), the predicted exponent is that of an affine deformation, $K\sim P/(V_{0}-V)\sim(z-z_{\rm c})^{2\alpha-2}$, different to that of $G$ given above. These different exponents have also been seen in simulations~\cite{OHern:03}, although it counters the consensus of lattice models, where all elastic moduli obey the same scaling behaviour~\cite{Head:03,Bergman:84,Feng:84,Arbabi:88}. It also disagrees with the MMA exponents presented earlier~(\ref{e:KGDef}). It is not yet clear if this is a real physical phenomenon or a coincidence of anomalous simulation results and naive choices for $P$ and $V$ in this theory. In any case, we are forced to conclude that the MMA predictions only apply to non--affine deformations.

The similarity of the MMA exponent for $G$ and that of simulations is noteworthy, as the theory relied on a local closure of equations and hence cannot incorporate the long wavelength modes observed in simulations~\cite{Silbert:05}. Persisting with the continuous phase transition analogy, it is possible that the scaling regime for non--trivial exponents is incredibly narrow; this is not without precedent, as ionic fluids exhibit mean field exponents except very near the critical temperature~\cite{Fisher:96,Luijten:02}. Alternatively the upper critical dimension may simply be~2. It should be noted that most (but not all~\cite{OHern:03}) simulations maintain a fixed system size, which will also give mean field exponents when the correlation length exceeds the system size.

A brief survey of other simulations not shown in the table indicate what new features will alter the scaling behaviour and hence represent relevant perturbations. Friction is an obvious candidate, and indeed recent simulations of Zhang {\em et al.}~\cite{Zhang:05} demonstrate that infinite friction will alter the exponents. They also show that finite friction introduces history dependency, suggesting a more advanced description of the dynamical phase will be required for a full theory. Also, the molecular dynamics simulations of Kasahara and Nakanishi~\cite{Kasahara:04a,Kasahara:04b} appear to find exponents consistent with simple rationals that are nonetheless different to those predicted here. This may be because their system has gravity, and hence the contact network is anisotropic and the force balance equations couple directly to an external field, either of which may be a relevant perturbation.


% Distribution of overlaps
\subsection{Distributed contact forces}

Of all the enhancements to the MMA theory that could be incorporated, perhaps the most pressing is to relax the assumption that every contact force~$f$, or equivalently every particle overlap $\delta=r_{0}-r$, is the same. Simulations have demonstrated that these distributions are in fact continuously distributed right down to zero forces or overlaps~\cite{Silbert:02,Ouaguenouni:97,OHern:01}. Below we present some calculations that probe the effect of polydisperse~$\delta$ within the MMA framework. Our conclusion is that it in fact makes very little difference to the overall behaviour of the system, and does not change the exponents already quoted. We are also able to confirm the scaling of the overlap distribution as observed in simulations~\cite{OHern:03}.

The simplest way to incorporate distributed overlaps is to retain the approximations leading to (\ref{e:forcebal_z_n}), namely independent isotropic bond orientations $\hat{\bf n}$ and a local coordination number independent of the contact forces, and perform the integral over a known, fixed distribution~$P(\delta)$. This can be performed explicitly for a convenient choice of parameters, for instance Hookean interactions $\fexp=1$ with a uniform overlap distribution  $P(\delta)=1/\delta_{0}$ for $0\leq\delta\leq\delta_{0}$. This then produces an equation for $\lambda$ that collapses into the form already studied (\ref{e:mma_lambda}) in the limit  $\delta_{0}\rightarrow0$, with the separation $r$ replaced by a mean overlap $r_{0}-\delta_{0}/2$. Given the volume function (\ref{e:exponent_b}) still holds, we see that $\delta$ being distributed uniformly down to $\delta=0$ does not significantly alter the statics or dynamics of the system in this case, but merely modifies the prefactors.

More general distributions can be considered by assuming that the behaviour observed for the monodisperse case still broadly applies. Specifically, we assume that a unique stressless rigidity transition exists at some point $z=z_{\rm c}$, and a boundary between stable and unstable compressive regimes extends continuously from this point into the region $z>z_{\rm c}$, similar to Fig.~\ref{f:CompTens}. Furthermore, we assume that the system becomes kinetically arrested on this boundary under energy minimisation. Therefore the averaged compliance $\lambda$ is expected to scale with the distance from the transition $\varepsilon=(z-z_{\rm c})/z_{\rm c}>0$ as

\begin{equation}
\lambda
=
\lambda_{0}
\varepsilon^{-\polyexptwo}
\label{e:lambda_ansatz}
\end{equation}

\noindent{}for small $\varepsilon$, with $\polyexptwo$ an unknown positive exponent. Then we make the following scaling {\em ansatz} for the distribution of overlaps $P(\delta)$,

\begin{equation}
P(\delta)
=
\varepsilon^{-\polyexp}
q(\varepsilon^{-\polyexp}\delta)
\label{e:polyexp}
\end{equation}

\noindent{}where $q(x)$ is a fixed distribution. (\ref{e:polyexp}) states that the distribution of overlaps will uniformly contract as the transition is approached, with a width $s$ that vanishes as $s\sim\varepsilon^{\gamma}$ with $\gamma>0$.

Inserting (\ref{e:lambda_ansatz}) and (\ref{e:polyexp}) into (\ref{e:forcebal_z_n}) allows the integration to be performed. Different results occur for different combinations of exponents, but only the form of solution for $\fexp\polyexp>\polyexptwo$ and $\polyexp(\fexp-1)<\polyexptwo$ admits a stressless rigidity transition. The equation for $\lambda_{0}$ in this case is

\begin{eqnarray}
1-\frac{d}{z}
&=&
(d-1)
\frac{\lambda_{0}\mu}{r_{0}}
\varepsilon^{-\polyexptwo+\fexp\polyexp}
\left\langle x^{\fexp}\right\rangle_{q(x)}
\nonumber\\
&+&
\frac{1}{\lambda_{0}\mu\fexp}
\varepsilon^{\polyexptwo-\polyexp(\fexp-1)}
\left\langle
x^{1-\fexp}
\right\rangle_{q(x)}
\label{e:poly_lambda}
\end{eqnarray}

\noindent{}where the angled brackets here denote averaging over the fixed distribution~$q(x)$. (\ref{e:poly_lambda}) admits a solution $\varepsilon=0$ at $z=z_{\rm c}=d$ given the exponent inequalities just quoted. It is similar to the monodisperse expression (\ref{e:mma_lambda}) with $f$ replaced by $\left\langle x^{\fexp}\right\rangle_{q(x)}$, and $\left\langle x^{1-\fexp}\right\rangle_{q(x)}$ the quantity related to $f^{\prime}$.
 

Stable solutions of (\ref{e:poly_lambda}) correspond to $\lambda_{0}$ real and positive, and so the boundary between the stable and unstable regions can be found using the familiar quadratic equation formulae. For consistency with the earlier assumption that the compressed stability boundary for $\left\langle x^{\fexp}\right\rangle_{q(x)}>0$
is continuously connected to the point $z=z_{\rm c}$, $\varepsilon=0$, we find it is necessary to impose $\gamma=2$. This means that the width of the overlap distribution scales as $s\sim(z-z_{\rm c})^{2}$ near the transition. The simulations of O'Hern {\em et al.} have shown that $s\sim(\phi-\phi_{\rm c})^{\Delta}$ with $\Delta$ close to 1 for Hookean interactions $\fexp=1$ in $d=3$~\cite{OHern:03}. According to Table~\ref{t:exponents}, this corresponds to $s\sim(z-z_{\rm c})^{2\Delta}$, in agreement with the prediction $\gamma=2\Delta=2$ found here. A second consistency check is that $\lambda_{0}={\mathcal O}(1)$, which allows the second exponent to be fixed, $\nu=2\fexp-1$. Note that both of these exponents are equal to their counterparts in the monodisperse--$\delta$ case, {\em i.e.} (\ref{e:zw_boundary}) for $\gamma$ and (\ref{e:lambda_boundary}) for $\nu$. Although this analysis in non--rigourous, it strongly suggests that polydisperse contacts does not alter the scaling picture presented earlier. 


% Exponents table v2 : vertical
 \begin{table}
\caption{\label{t:exponents}Table of the scaling relations between $\Delta z=z-z_{\rm c}$, $G$, $P$ and $\Delta V=V_{0}-V\sim(\phi-\phi_{0})$ (with $\phi_{0}$ the critical volume fraction) as predicted by the MMA theory, where $\alpha$ is the force law exponent~(\ref{e:f_r}), $d$ dimension. For comparison, results from simulations of central force systems and the trivial (affine) predictions are also shown.}
\begin{ruledtabular}
\begin{tabular}{c@{\quad}|@{\quad}ccc}
Model &
$G\sim\Delta z^{\Kzexp}$ &
$P\sim\Delta z^{\Pzexp}$ &
$\Delta V\sim\Delta z^{\Vzexp}$
\\
\hline
%% MMA approx (here)
MMA\\
$\alpha>0$, $d\geq2$ & $2\alpha-1$ & $2\alpha$ & 2 \\
%$\left\{\begin{array}{c@{\quad:\quad}l}2b&b<1\\2&b\geq1\end{array}\right.$
\hline
Affine (see~\cite{OHern:03}) & $2\alpha-2$ & $2\alpha$ & -\\ 
\hline
%% Doug Durian's wet foam model
Wet foam~\cite{Durian:Foam} \\
$\alpha=1$, $d=2$ & $\approx1$\footnote{Result for shear modulus shown.}  & $2\pm0.4$ & $2\pm0.4$ \\
\hline
%% Corey's min. model
O'Hern {\em et al.}~\cite{OHern:03} \\
$\alpha=1$, $d=2,3$ & $1.01\pm0.1^{a}$ & $2.1\pm0.2$ & $2.04\pm0.1$ \\
$\alpha=3/2$, $d=2,3$ & $2.08\pm0.1^{a}$ & $3.15\pm0.3$ & $2.08\pm0.1$ \\
\hline
%% Zhang/Makse emulsion/frictional granular model
Zhang {\em et al.}~\cite{Zhang:05} \\
$\alpha=1.28$, $d=3$ & - & $\approx2.45$ & $\approx1.96$ \\
%$\alpha=3/2$\footnote{Results for infinite friction case shown.}, $d=3$
%& - & $\approx3.69$ & - & $\approx2.22$ \\
%$\alpha=3/2$\footnote{Results for finite friction case shown.}, $d=3$
%& - & $\approx3.25$ & - & $\approx2.17$ \\
\hline
% Makse EMA paper
Makse {\em et al.}~\cite{Makse:05}\footnote{Only frictionless data shown.} \\
%$\alpha=3/2$, $d=2$ & - & $3.6\pm0.9$ & $1.43\pm0.7$ \\
$\alpha=3/2$, $d=3$ & - & $3.3\pm0.5$ & $2.1\pm0.6$ \\
\end{tabular}
\end{ruledtabular}
\end{table}






%%
%% Discussion
%%
\section{Discussion}
\label{s:discussion}

A central feature of this class of problem is the intrinsic interweaving between the dissipative dynamics as the system cools, with the mechanical response of the arrested state. In this paper, a minimal coupling has been presented, namely that the dynamics proceeds in an independent--particle manner until a mechanically stable region is reached, Fig.~\ref{f:Min_PV}. This is highly simplified, and a more elaborate theory is desirable, perhaps along the lines of the bootstrap percolation model approach~\cite{Schwarz:05}. Intuitively, we expect that, during the dynamic phase, transient overconstrained, rigid clusters will become stressed by interparticle collisions, and relieve this stress by becoming non--rigid, {\em i.e.} expanding into neighbouring, underconstrained regions. Kinetic arrest occurs when a spanning rigid cluster forms. Note that this argument suggests a dynamic homogenisation process, perhaps explaining the apparent appearance of mean field exponents in $d=2$ and 3 simulations in table~\ref{t:exponents}.

This first application of the mean mode approximation has been applied to arguably the simplest particulate problem, namely repulsive central forces in an isotropic system. The simplicity of its results suggests that additional features could be included while remaining tractable. For instance, friction and gravity would be needed before any sensible comparison with real granular media could be made. A problem that may emerge is closing the equations; the averaged force balance equation (\ref{e:force_balance}) only gives one scalar equation, which is why the proposed displacement modes were parameterised by single scalar~$\lambda$. If a future application had too few equations, one possible approach would be to assume the response will minimise the increase in elastic energy, converting it to a minimisation problem with any known equations as constraints. In principle, displacement modes with any number of unknown parameters could be introduced by this approach.




% If you have acknowledgments, this puts in the proper section head.
\begin{acknowledgments}

%The author gratefully acknowledges discussions with M. Cates,
DAH was jointly funded by the European Union Marie Curie program and the JSPS (Japanese Society for the Promotion of Science) program. The author would also like to thank David Dean for suggesting references~\cite{Fisher:96,Luijten:02}.

\end{acknowledgments}




% Specify following sections are appendices. Use \appendix* if there
% only one appendix.
%\appendix*
%\section{Energetic description of stability with prestress}
%A static system must be at a turning point on its potential energy landscape $\PE(\{{\bf x}^{\alpha}\})$, {\em i.e.} $\frac{\partial\PE}{\partial x^{\alpha}_{i}}=0$ $\forall \alpha=1\dots N$, $i=1\dots d$. Linear stability requires the stronger constraint that it is also at a local minimum, so that the second derivatives are strictly positive. For an arbitrary infinitesimal change in particle positions ${\bf x}^{\alpha}\rightarrow{\bf x}^{\alpha}+\gamma {\bf u}^{\alpha}$ with $\gamma\ll1$, and ignoring contact formation or scission (which should only make ${\mathcal O}(\gamma^{2})$ contribution to the total energy, since each such bond is by definition weak), the total change in potential energy is
%\begin{equation}
%\Delta\PE
%=
%-
%\frac{1}{2}
%\gamma^{2}
%\sum_{{\rm all\,bonds\,} \alpha\beta}
%\Delta{\bf u}^{\alpha\beta}
%A^{\alpha\beta}
%\Delta{\bf u}^{\alpha\beta}
%\label{e:work}
%\end{equation}
%\noindent{}in terms of the relative displacements between contacting particles $\Delta{\bf u}^{\alpha\beta}={\bf u}^{\alpha}-{\bf u}^{\beta}$ and the same matrix $A$ as in~(\ref{e:A_ij}). For $f^{\prime}(r)<0$ this can be rewritten as
%\begin{eqnarray}
%\Delta\PE
%&=&
%\frac{\gamma^{2}}{2}
%\left\{
%\sum_{(\alpha\beta)}
%|f^{\prime}(r^{\alpha\beta})|
%\left(\Delta u^{[\alpha\beta]}_{\parallel}\right)^{2}
%\right.
%\nonumber\\
%&&-
%\left.
%\sum_{(\alpha\beta)}
%\frac{f(r^{\alpha\beta})}{r^{\alpha\beta}}
%\left(\Delta u^{[\alpha\beta]}_{\perp}\right)^{2}
%\right\}
%\label{e:E_prestressed}
%\end{eqnarray}
%\noindent{}where $\Delta u_{\perp,\parallel}$ are the components of $\Delta {\bf u}$ perpendicular and parallel to $\hat{\bf n}$, respectively. The three cases of interest are:
%\vskip \baselineskip
%\noindent{\em Stressless:} For all $f\equiv0$, the truncated potentials of (\ref{e:f_r}) have $f^{\prime}$ either zero or undefined, neither of which is suitable for studying the linear response. For non--truncated Hookean springs, however, (\ref{e:E_prestressed}) suggests a rigidity transition when not all $\Delta u_{\parallel}$ vanish, which can be estimated {\em via} the usual constraint counting methods. 
%\vskip \baselineskip
%\noindent{\em Tension:} For all $f<0$, (\ref{e:E_prestressed}) is positive unless both $\Delta u_{\parallel}$ and $\Delta u_{\perp}$ simultaneously vanish for each bond, which only occurs in the trivial case of a disconnected system. Hence tensile systems are stable.
%\vskip \baselineskip
%\noindent{\em Compression:} When all $f>0$, for which the particulate potentials (\ref{e:f_r}) are once more valid, the two terms in (\ref{e:E_prestressed}) have opposite sign, signifying competing influences. In contrast to tension, compression acts to destabilise the system, although it is not possible to calculate the stability boundary from this general analysis.


%%
%%
%%

%Some insight into the nature of the transition can be gleamed by observing that the Maxwell counting method described above is {\em never} exact for the systems of interest here. For particulate matter, isolated non--bonded particles or {\em rattlers} must first be removed before it comes within error bars~\cite{Donev_condmat,OHern:03,Zhang:05}. For disordered lattices the estimate generally overestimates the $z_{\rm c}$ determined in simulations by a few percent or more, significantly beyond errors, a situation that persists when transverse forces are included~\cite{Schwartz:85,LatvaKokko:01,Head:03}. This is {\em not} due to colinearity or coplanarity of bonds in the lattice, as efficient numerical schemes are used (at least for $d=2$) which effectively perturb each node by a small amount in a random direction, ensuring linear independence~\cite{Jacobs:95}. It is rather due to the implicit assumption that the set of equations imposed by the constraints are non--decomposable. That is, they cannot be separated into two or more independent subsets of equations, corresponding to regions of the system that are not rigidly connected. In fact this is not true: for particulate matter, the rattler particles mentioned above each represent an isolated, independent subset of equations. That $z_{\rm c}=2d$ appears to be recovered after their removal suggests only a single, compact, rigidly--connected cluster remains. For disordered lattices, their are a large number of such subsystems; indeed, even the largest, system--spanning rigid cluster has a fractal dimension at the rigidity transition and thus occupies a vanishingly small fraction of the total system volume in the asymptotic limit. This also means that the system--spanning rigid cluster first appears when (almost) all of the system is non--rigid, clearly demonstrating why estimates based on just the mean coordination number will always overestimate the true location of the transition.



%: References
%%
%% REFERENCES
%%
\begin{thebibliography}{99}

% RoySoc granular
\bibitem{Cates:98} M. E. Cates, J. P. Wittmer, J. --P. Bouchaud and P. Claudin, Phil. Trans. R. Soc. Land. A {\bf 356}, 2535 (1998); also {\em cond-mat/9803266}.

% Greens exps (long)
\bibitem{Geng:03} J. Geng, G. Reydellet, E. Cl\'ement and R. P. Behringer, Physica D {\bf 182}, 274 (2003).

% Moukarzel greens
\bibitem{Moukarzel:04} C. F. Moukarzel, H. Pacheco--Mart\'inez, J. C. Ruiz--Suarez and A. M. Vidales, Gran. Mat. {\bf 6}, 61 (2004).

% Kasahara & Nakanishi sims
\bibitem{Kasahara:04a} A. Kasahara and H. Nakanishi, Phys. Rev. E {\bf 70}, 051309 (2004).

\bibitem{Kasahara:04b} A. Kasahara and H. Nakanishi, J. Phys. Soc. Japan {\bf 73}, 789 (2004).

% R. Stinchcombe et al.'s "spatially heterogeneous dynamics" preprint
%\bibitem{Lefevre:04} A. Lef\`evre, L. Berthier and R. Stinchcombe, {\em cond-mat/0421741}.

% Compaction papers
\bibitem{Barker:93} G. C. Barker and A. Mehta, Phys. Rev. E {\bf 47}, 184 (1993).

\bibitem{Head:00} D. A. Head, Phys. Rev. E {\bf 62}, 2439 (2000).

\bibitem{Rafi:DryFoam} R. Blumenfeld, J. Phys. A {\bf 36}, 2399 (2003).

\bibitem{Durian:Foam} D. J. Durian, Phys. Rev. Lett {\bf 75}, 4780 (1995).

% Foam : early transition paper
\bibitem{Bolton:90} F. Bolton and D. Weaire, Phys. Rev. Lett. {\bf 65}, 3449 (1990).

% Emulsion: rheology exps
\bibitem{Mason:95} T. G. Mason, J. Bibette and D. A. Weitz, Phys. Rev. Lett. {\bf 75}, 2051 (1995).

% Makse/Zhang recent emulsion/granular paper (some with finite friction)
\bibitem{Zhang:05} H. P. Zhang and H. A. Makse, {\em cond-mat/0501370}.

% Calladine - Maxwell counting
\bibitem{Calladine:78} C. R. Calladine, Int. J. Solids Structures {\bf 14}, 161 (1978).

% Connelly/Whitely: definitions of (2nd order) rigidity
\bibitem{Connelly:96} R. Connelly and W. Whitely, SIAM J. Disc. Math. {\bf 9}, 453 (1996).

% Guest/Hutchinson, repetitive structures
\bibitem{Guest:03} S. D. Guest and J. W. Hutchinson, J. Mech. Phys. Sol. {\bf 51}, 383 (2003).

\bibitem{Sahimi:Book} M. Sahimi, {\em Heterogeneous Materials I: Linear transport and optical properties} (Springer--Verlag, New York, 2003).

% Original pebble game paper
\bibitem{Jacobs:95} D. J. Jacobs and M. F. Thorpe, Phys. Rev. Lett. {\bf 75}, 4051 (1995).

% Finite element of rigidly bonded fibre nets
\bibitem{Astrom:00} J. A. \AA strom, J. P. M\"akinen, M. J. Alava and J. Timonen, Phys. Rev. E {\bf 61}, 5550 (2000).

% Pebble game for filaments
\bibitem{LatvaKokko:01} M. Latva--Kokko and J. Timonen, Phys. Rev. E {\bf 64}, 066117 (2001).

% Earlier filament rigidity paper
\bibitem{Head:03} D. A. Head, F. C. MacKintosh and A. J. Levine, Phys. Rev. E 68, 025101(R) (2003).

% Universal Poisson at transition paper
\bibitem{Bergman:84} D. J. Bergman and Y. Kantor, Phys. Rev. Lett. {\bf 53}, 511 (1984).

% Feng's BB model
\bibitem{Feng:84} S. Feng, P. N. Sen, B. I. Halperin and C. J. Lobb, Phys. Rev. B {\bf 30}, 5386 (1984).

% Arbabi/Sahaimi P(f) distribution modes
\bibitem{Arbabi:88} S. Arbabi and M. Sahimi, Phys. Rev. B {\bf 38}, 7173 (1988).

% Tensed lattice 1
\bibitem{Tang:88} W. Tang and M. F. Thorpe, Phys. Rev. B {\bf 37}, 5539 (1988).

% Tensed lattice 2
\bibitem{Zhou:03} Z. Zhou, B. Jo\'os and P. Lai, Phys. Rev. E {\bf 68}, 055101(R) (2003).

\bibitem{footnote} This definition of rigidity was chosen for its experimental practicality, although other definitions can be entertained; see {\em e.g.}~\cite{Connelly:96}.

% Comment/reply about relationship to packing models
\bibitem{Donev:04} A. Donev, S. Torquato, F. H. Stillinger and R. Connelly, Phys. Rev. E {\bf 70}, 043301 (2004).

\bibitem{OHern:04} C. O'Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E {\bf 70}, 043302 (2004).


% Point J - style papers
% Point-J long paper
\bibitem{OHern:03} C. O'Hern, L. E. Silbert, A. J. Liu and S. R Nagel, Phys. Rev. E {\bf 68}, 011306 (2003).

% Diverging length scale in point J
\bibitem{Silbert:05} L. E. Silbert, A. J. Liu and S. R. Nagel, {\em cond-mat/0501616}.

% Rigidity phase transition - constant V & P
\bibitem{Aharonov:99} E. Aharonov and D. Sparks, Phys. Rev. E {\bf 60}, 6890 (1999).

% Makse et al - EMA Fail
\bibitem{Makse:EMAFail} H. A. Makse, N. Gland, D. L. Johnson and L. S. Schwartz, Phys. Rev. Lett. {\bf 83}, 5070 (1999).

\bibitem{Makse:05} H. A. Makse, N. Gland, D. L. Johnson and L. Schwarz, Phys. Rev. E {\bf 70}, 061302 (2004).

% Hard sphere
\bibitem{Donev_condmat} A. Donev, S. Torquato and F. H. Stillinger, {\em cond-mat/0408550}.

% A. Liu et als. recent bootstrap percolation paper
\bibitem{Schwarz:05} J. M. Schwarz, A. Liu and L. Q. Chayes, {\em cond-mat/0410595}.

% V. recent field theory paper
\bibitem{Henkes:05} S. Henkes and B. Chakraborty, {\em cond-mat/0504371}.

% itshhhmeee....
\bibitem{Head:04} D. A. Head, J. Phys. A {\bf 37}, 10771 (2004).

% Weiner book
\bibitem{Weiner:Book} J. H. Weiner, {\em Statistical Mechanics of Elasticity}, 2nd ed. (Dover Publications, New York, 2002).

% Granular "EMA" (friction; no transition)
\bibitem{Walton:87} K. Walton, J. Mech. Phys. Solids {\bf 35}, 213 (1987).

% Ring-removal / self-organisation thy paper
%\bibitem{Thorpe:00} M. F. Thorpe, D. J. Jacobs, M. V. Chubynsky and J. C. Phillips, J. Non--Cryst. Sol. {\bf 266}, 859 (2000).

% Affine : energy principles
\bibitem{Kruyt:98} N. P. Kruyt and L. Rothenburg, Int. J. Eng. Sci. {\bf 36}, 1127 (1998).

% Affine: + non-affine modes
\bibitem{Trentadue:01} F. Trentadue, Int. J. Sol. Struct. {\bf 38}, 7319 (2001).

% EMA : central forces
\bibitem{Feng:85} S. Feng, M. F. Thorpe and E. Garboczi, Phys. Rev. B {\bf 31}, 276 (1985).
 
% EMA: transverse
\bibitem{Schwartz:85} L. M. Schwartz, S. Feng, M. F. Thorpe and P. N. Sen, Phys. Rev. B {\bf 32}, 4607 (1985).

% EMA: Sound propagation ; HCP
\bibitem{Velicky:02} B. Velick\'y and C. Caroli, Phys. Rev. E {\bf 65}, 021307 (2002).

% Moukarzel's isostaticity paper
%\bibitem{Moukarzel:98} C. F. Moukarzel, Phys. Rev. Lett. {\bf 81}, 1634 (1998).

\bibitem{LandauLifshitz} L. D. Landau and E. M. Lifshitz, {\em Theory of elasticity}, 3rd ed. (Butterworth--Heinemann, Oxford, 1986).

% Granular green's
\bibitem{Tanguy:02} A. Tanguy, J. P. Wittmer, F. Leonforte and J.--L. Barrat, Phys. Rev. B {\bf 66}, 174205 (2002).

\bibitem{AlexanderRev} S. Alexander, Phys. Rep. {\bf 296}, 65 (1998).

% Shells book
\bibitem{Libai&Simmonds} A. Libai and J. G. Simmonds, {\em The nonlinear theory of elastic shells}, 2nd ed. (Cambridge University Press, Cambridge, 1998).

% Recent preprint on gran stat mech - mentions V functions
\bibitem{Edwards:05} S. F. Edwards, J. Brujic/'c and H. A. Makse, {\em cond-mat/0503057}.

%  Silbert al. sims : P(f)
\bibitem{Silbert:02} L. E. Silbert, G. S. Grest and J. W. Landry, Phys. Rev. E {\bf 66}, 061303 (2002).

% JN Roux - P(f)
\bibitem{Ouaguenouni:97} S. Ouaguenouni and J.--N. Roux, Euro. Phys. Lett. {\bf 39}, 117 (1997)

% P(f) paper
\bibitem{OHern:01} C. O'Hern, S. A. Langer, A. J. Liu and S. R. Nagel, Phys. Rev. Lett. {\bf 86}, 111 (2001).

% Makse "compressible granular media" paper
%\bibitem{Makse:2000} H. A. Makse, D. L. Johnson and L. M. Schwartz, Phys. Rev. Lett. {\bf 84}, 4160 (2000).

% Generic OCTA reference
\bibitem{OCTA} Available at \verb"http://octa.jp"

% Narrow scaling regime example
\bibitem{Fisher:96} M. E. Fisher, J. Phys. Cond. Mat. {\bf 8}, 9103 (1996).

\bibitem{Luijten:02} E. Luijten, M. E. Fisher and A. Z. Panagiotopoulos, Phys. Rev. Lett. {\bf 88}, 185701 (2002).

\end{thebibliography}

\end{document}


%% Official blurb from APS:
%
% If in two-column mode, this environment will change to single-column
% format so that long equations can be displayed.
%\begin{widetext}
% (long equation)
%\end{widetext}

%% Figures: requires "graphics" or "graphicx" packages
% Use the figure* environment if the figure should span across the
% entire page. There is no need to do explicit centering.
% \begin{figure}
% \includegraphics{}%
% \caption{Caption\label{}
% \end{figure}




 
cond-mat0504418.extracted_bib
\begin{thebibliography}{99}

% RoySoc granular
\bibitem{Cates:98} M. E. Cates, J. P. Wittmer, J. --P. Bouchaud and P. Claudin, Phil. Trans. R. Soc. Land. A {\bf 356}, 2535 (1998); also {\em cond-mat/9803266}.

% Greens exps (long)
\bibitem{Geng:03} J. Geng, G. Reydellet, E. Cl\'ement and R. P. Behringer, Physica D {\bf 182}, 274 (2003).

% Moukarzel greens
\bibitem{Moukarzel:04} C. F. Moukarzel, H. Pacheco--Mart\'inez, J. C. Ruiz--Suarez and A. M. Vidales, Gran. Mat. {\bf 6}, 61 (2004).

% Kasahara & Nakanishi sims
\bibitem{Kasahara:04a} A. Kasahara and H. Nakanishi, Phys. Rev. E {\bf 70}, 051309 (2004).

\bibitem{Kasahara:04b} A. Kasahara and H. Nakanishi, J. Phys. Soc. Japan {\bf 73}, 789 (2004).

% R. Stinchcombe et al.'s "spatially heterogeneous dynamics" preprint
%\bibitem{Lefevre:04} A. Lef\`evre, L. Berthier and R. Stinchcombe, {\em cond-mat/0421741}.

% Compaction papers
\bibitem{Barker:93} G. C. Barker and A. Mehta, Phys. Rev. E {\bf 47}, 184 (1993).

\bibitem{Head:00} D. A. Head, Phys. Rev. E {\bf 62}, 2439 (2000).

\bibitem{Rafi:DryFoam} R. Blumenfeld, J. Phys. A {\bf 36}, 2399 (2003).

\bibitem{Durian:Foam} D. J. Durian, Phys. Rev. Lett {\bf 75}, 4780 (1995).

% Foam : early transition paper
\bibitem{Bolton:90} F. Bolton and D. Weaire, Phys. Rev. Lett. {\bf 65}, 3449 (1990).

% Emulsion: rheology exps
\bibitem{Mason:95} T. G. Mason, J. Bibette and D. A. Weitz, Phys. Rev. Lett. {\bf 75}, 2051 (1995).

% Makse/Zhang recent emulsion/granular paper (some with finite friction)
\bibitem{Zhang:05} H. P. Zhang and H. A. Makse, {\em cond-mat/0501370}.

% Calladine - Maxwell counting
\bibitem{Calladine:78} C. R. Calladine, Int. J. Solids Structures {\bf 14}, 161 (1978).

% Connelly/Whitely: definitions of (2nd order) rigidity
\bibitem{Connelly:96} R. Connelly and W. Whitely, SIAM J. Disc. Math. {\bf 9}, 453 (1996).

% Guest/Hutchinson, repetitive structures
\bibitem{Guest:03} S. D. Guest and J. W. Hutchinson, J. Mech. Phys. Sol. {\bf 51}, 383 (2003).

\bibitem{Sahimi:Book} M. Sahimi, {\em Heterogeneous Materials I: Linear transport and optical properties} (Springer--Verlag, New York, 2003).

% Original pebble game paper
\bibitem{Jacobs:95} D. J. Jacobs and M. F. Thorpe, Phys. Rev. Lett. {\bf 75}, 4051 (1995).

% Finite element of rigidly bonded fibre nets
\bibitem{Astrom:00} J. A. \AA strom, J. P. M\"akinen, M. J. Alava and J. Timonen, Phys. Rev. E {\bf 61}, 5550 (2000).

% Pebble game for filaments
\bibitem{LatvaKokko:01} M. Latva--Kokko and J. Timonen, Phys. Rev. E {\bf 64}, 066117 (2001).

% Earlier filament rigidity paper
\bibitem{Head:03} D. A. Head, F. C. MacKintosh and A. J. Levine, Phys. Rev. E 68, 025101(R) (2003).

% Universal Poisson at transition paper
\bibitem{Bergman:84} D. J. Bergman and Y. Kantor, Phys. Rev. Lett. {\bf 53}, 511 (1984).

% Feng's BB model
\bibitem{Feng:84} S. Feng, P. N. Sen, B. I. Halperin and C. J. Lobb, Phys. Rev. B {\bf 30}, 5386 (1984).

% Arbabi/Sahaimi P(f) distribution modes
\bibitem{Arbabi:88} S. Arbabi and M. Sahimi, Phys. Rev. B {\bf 38}, 7173 (1988).

% Tensed lattice 1
\bibitem{Tang:88} W. Tang and M. F. Thorpe, Phys. Rev. B {\bf 37}, 5539 (1988).

% Tensed lattice 2
\bibitem{Zhou:03} Z. Zhou, B. Jo\'os and P. Lai, Phys. Rev. E {\bf 68}, 055101(R) (2003).

\bibitem{footnote} This definition of rigidity was chosen for its experimental practicality, although other definitions can be entertained; see {\em e.g.}~\cite{Connelly:96}.

% Comment/reply about relationship to packing models
\bibitem{Donev:04} A. Donev, S. Torquato, F. H. Stillinger and R. Connelly, Phys. Rev. E {\bf 70}, 043301 (2004).

\bibitem{OHern:04} C. O'Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E {\bf 70}, 043302 (2004).


% Point J - style papers
% Point-J long paper
\bibitem{OHern:03} C. O'Hern, L. E. Silbert, A. J. Liu and S. R Nagel, Phys. Rev. E {\bf 68}, 011306 (2003).

% Diverging length scale in point J
\bibitem{Silbert:05} L. E. Silbert, A. J. Liu and S. R. Nagel, {\em cond-mat/0501616}.

% Rigidity phase transition - constant V & P
\bibitem{Aharonov:99} E. Aharonov and D. Sparks, Phys. Rev. E {\bf 60}, 6890 (1999).

% Makse et al - EMA Fail
\bibitem{Makse:EMAFail} H. A. Makse, N. Gland, D. L. Johnson and L. S. Schwartz, Phys. Rev. Lett. {\bf 83}, 5070 (1999).

\bibitem{Makse:05} H. A. Makse, N. Gland, D. L. Johnson and L. Schwarz, Phys. Rev. E {\bf 70}, 061302 (2004).

% Hard sphere
\bibitem{Donev_condmat} A. Donev, S. Torquato and F. H. Stillinger, {\em cond-mat/0408550}.

% A. Liu et als. recent bootstrap percolation paper
\bibitem{Schwarz:05} J. M. Schwarz, A. Liu and L. Q. Chayes, {\em cond-mat/0410595}.

% V. recent field theory paper
\bibitem{Henkes:05} S. Henkes and B. Chakraborty, {\em cond-mat/0504371}.

% itshhhmeee....
\bibitem{Head:04} D. A. Head, J. Phys. A {\bf 37}, 10771 (2004).

% Weiner book
\bibitem{Weiner:Book} J. H. Weiner, {\em Statistical Mechanics of Elasticity}, 2nd ed. (Dover Publications, New York, 2002).

% Granular "EMA" (friction; no transition)
\bibitem{Walton:87} K. Walton, J. Mech. Phys. Solids {\bf 35}, 213 (1987).

% Ring-removal / self-organisation thy paper
%\bibitem{Thorpe:00} M. F. Thorpe, D. J. Jacobs, M. V. Chubynsky and J. C. Phillips, J. Non--Cryst. Sol. {\bf 266}, 859 (2000).

% Affine : energy principles
\bibitem{Kruyt:98} N. P. Kruyt and L. Rothenburg, Int. J. Eng. Sci. {\bf 36}, 1127 (1998).

% Affine: + non-affine modes
\bibitem{Trentadue:01} F. Trentadue, Int. J. Sol. Struct. {\bf 38}, 7319 (2001).

% EMA : central forces
\bibitem{Feng:85} S. Feng, M. F. Thorpe and E. Garboczi, Phys. Rev. B {\bf 31}, 276 (1985).
 
% EMA: transverse
\bibitem{Schwartz:85} L. M. Schwartz, S. Feng, M. F. Thorpe and P. N. Sen, Phys. Rev. B {\bf 32}, 4607 (1985).

% EMA: Sound propagation ; HCP
\bibitem{Velicky:02} B. Velick\'y and C. Caroli, Phys. Rev. E {\bf 65}, 021307 (2002).

% Moukarzel's isostaticity paper
%\bibitem{Moukarzel:98} C. F. Moukarzel, Phys. Rev. Lett. {\bf 81}, 1634 (1998).

\bibitem{LandauLifshitz} L. D. Landau and E. M. Lifshitz, {\em Theory of elasticity}, 3rd ed. (Butterworth--Heinemann, Oxford, 1986).

% Granular green's
\bibitem{Tanguy:02} A. Tanguy, J. P. Wittmer, F. Leonforte and J.--L. Barrat, Phys. Rev. B {\bf 66}, 174205 (2002).

\bibitem{AlexanderRev} S. Alexander, Phys. Rep. {\bf 296}, 65 (1998).

% Shells book
\bibitem{Libai&Simmonds} A. Libai and J. G. Simmonds, {\em The nonlinear theory of elastic shells}, 2nd ed. (Cambridge University Press, Cambridge, 1998).

% Recent preprint on gran stat mech - mentions V functions
\bibitem{Edwards:05} S. F. Edwards, J. Brujic/'c and H. A. Makse, {\em cond-mat/0503057}.

%  Silbert al. sims : P(f)
\bibitem{Silbert:02} L. E. Silbert, G. S. Grest and J. W. Landry, Phys. Rev. E {\bf 66}, 061303 (2002).

% JN Roux - P(f)
\bibitem{Ouaguenouni:97} S. Ouaguenouni and J.--N. Roux, Euro. Phys. Lett. {\bf 39}, 117 (1997)

% P(f) paper
\bibitem{OHern:01} C. O'Hern, S. A. Langer, A. J. Liu and S. R. Nagel, Phys. Rev. Lett. {\bf 86}, 111 (2001).

% Makse "compressible granular media" paper
%\bibitem{Makse:2000} H. A. Makse, D. L. Johnson and L. M. Schwartz, Phys. Rev. Lett. {\bf 84}, 4160 (2000).

% Generic OCTA reference
\bibitem{OCTA} Available at \verb"http://octa.jp"

% Narrow scaling regime example
\bibitem{Fisher:96} M. E. Fisher, J. Phys. Cond. Mat. {\bf 8}, 9103 (1996).

\bibitem{Luijten:02} E. Luijten, M. E. Fisher and A. Z. Panagiotopoulos, Phys. Rev. Lett. {\bf 88}, 185701 (2002).

\end{thebibliography}