Jamming II: Edwards' statistical mechanics of random packings of hard spheres
aa r X i v : . [ c ond - m a t . s o f t ] N ov Jamming II: Edwards’ statistical mechanics of randompackings of hard spheres
Ping Wang a , Chaoming Song a , Yuliang Jin a , Hern´an A. Makse a, ∗ a Levich Institute and Physics Department, City College of New York, New York, NY10031, US
AbstractThe problem of finding the most efficient way to pack spheres hasan illustrious history, dating back to the crystalline arrays conjec-tured by Kepler and the random geometries explored by Bernal inthe 60’s. This problem finds applications spanning from the mathe-matician’s pencil, the processing of granular materials, the jammingand glass transitions, all the way to fruit packing in every grocery.There are presently numerous experiments showing that the loosestway to pack spheres gives a density of ∼ % (named random loosepacking, RLP) while filling all the loose voids results in a maximumdensity of ∼ − % (named random close packing, RCP). Whilethose values seem robustly true, to this date there is no well-acceptedphysical explanation or theoretical prediction for them. Here we de-velop a common framework for understanding the random packingsof monodisperse hard spheres whose limits can be interpreted as theexperimentally observed RLP and RCP. The reason for these limitsarises from a statistical picture of jammed states in which the RCPcan be interpreted as the ground state of the ensemble of jammedmatter with zero compactivity, while the RLP arises in the infinitecompactivity limit. We combine an extended statistical mechanicsapproach ’a la Edwards’ (where the role traditionally played by theenergy and temperature in thermal systems is substituted by thevolume and compactivity) with a constraint on mechanical stabilityimposed by the isostatic condition. We show how such approachescan bring results that can be compared to experiments and allowfor an exploitation of the statistical mechanics framework. The keyresult is the use of a relation between the local Voronoi volumes ofthe constituent grains (denoted the volume function) and the num-ber of neighbors in contact that permits to simple combine the twoapproaches to develop a theory of volume fluctuations in jammed ∗ Corresponding author
Email address: [email protected] (Hern´an A. Makse)
Preprint submitted to Elsevier October 27, 2018 atter. Ultimately, our results lead to a phase diagram that pro-vides a unifying view of the disordered hard sphere packing problemand further shedding light on a diverse spectrum of data, includingthe RLP state. Theoretical results are well reproduced by numericalsimulations that confirm the essential role played by friction in deter-mining both the RLP and RCP limits. The RLP values depend onfriction, explaining why varied experimental results can be obtained.1. Introduction
Filling containers with balls is one of the oldest mathematical puzzles knownto scientists. The study of disordered sphere packings raises an interestingproblem: How efficient and uniform will spheres pack if assembled randomly?This problem has an important application in the jamming phenomenon, whichtakes place in particulate systems where all particles are in close contact withone another. The study of jammed granular media offers unexpected chal-lenges in physics since the equilibrium statistical mechanics fail for these out-of-equilibrium systems. The goal of the present study is to develop an ensembletheory of volume fluctuations to describe the statistical mechanics of jammedmatter with the aim of shedding light to the long-standing problem of charac-terizing random sphere packings.
The study of sphere packing problem started four centuries ago when Jo-hannes Kepler conjectured that the most efficient arrangement of spheres is theFCC lattice (an important part of the 18th problem proposed by Hilbert in1900). Even though this is a tool used for centuries in fruit markets around theglobe, nearly 400 years passed before this conjecture was considered a mathe-matical proof, which has been developed only recently by Hales [1] in a seriesof articles covering 250 pages supplemented by 3Gb of computer code to deter-mine the best ordered packing through linear programming. The difficulty arisessince in 3d it is not enough to look at the packing of one cell, but it is necessaryto consider several Voronoi cells at once. That is, the packing that minimizethe volume locally (the dodecahedron) does not tile the system globally. Sucha situation does not arise in 2d, where the hexagonal packing minimizes thevolume locally and globally; the equivalent of Kepler conjecture in 2d, which isalso known as Thue’s theorem, was proved long ago [2, 3].The analogous problem for disordered packings has also an illustrious butunfinished history. This problem was initiated by the pioneering work of Bernalin the 1960’s [4] (although earlier attempts can be found). The traditionalview states that [5] “packings of spherical particles have been shaken, settled indifferent fluids and kneaded inside rubber balloons and all with no better resultsthan a maximum density of 63%”. This is the so-called random close packingRCP limit [4, 5, 6, 7, 8]. On the other hand, other experiments have shown thatdensities as low as 55% can be obtained in random loose packings, RLP [6, 8, 9].2hile the two limits seem reproducible in various experiments and computersimulations, a mathematically-rigorous definition is still unavailable. Indeedit was conjectured [10] that the RCP conception is mathematically ill-definedand should be replaced by the maximally random jammed (MRJ) conceptionin terms of an ensemble of order parameters.To this date there is no well-accepted physical explanation of this phe-nomenon, no well-accepted theoretical prediction of such density values andheated debates are still found in the literature regarding the existence of rigor-ous definitions of the RCP and RLP, the uniqueness of the RCP state and thenature of their state of randomness.
One of the most important physical applications of the sphere packing prob-lem is to understand the intrinsic geometrical structures of jammed systems.A great deal of research has been devoted by physicists to the properties ofpackings of granular materials and other jammed systems such as compressedemulsions, dense colloids and glasses [12, 11], due to their “out-of-equilibrium”nature. Generally speaking, there are three interrelated approaches to under-stand the nature of the jammed state of matter and its ensuing jamming tran-sition : (a) Statistical mechanics of jammed matter : The jammed phase isdescribed with the ensemble of volume and force fluctuations proposed by Ed-wards [13, 14, 16, 15, 21, 17, 18, 20, 19]. Ideas from the physics of glassesprovide interesting cross-fertilization between the jamming transition and theglass transition [22, 23, 24, 25]. (b) Critical phenomena of deformable particle : jamming is seen asa critical point at which observables such as pressure, coordination numberand elastic moduli behave as power-law near the critical volume fraction, φ c ,[26, 27, 29, 31, 32, 30, 28, 33, 34]. The jammed state is modelled as granularmedia composed of soft particles, typically Hertz-Mindlin [32, 35, 36] grainsunder external pressure and emulsions compressed under osmotic pressure [37],as they approach the jamming transition from the solid phase above the criticaldensity: φ → φ + c . (c) Hard-sphere glasses : In parallel to studies in the field of jammingthere exists a community attempting to understand the packing problem ap-proaching jamming from the liquid phase [38, 39, 40], that is, φ → φ − c . Here,amorphous jammed packings are seen as infinite pressure glassy states [38, 41].Therefore, the properties of the jamming transition are intimately related tothose of the glass transition [38].A great deal of effort has been devoted to the study of jamming apply-ing these approaches. A large body of experiments and simulations have fullycharacterized the jammed state in the critical phenomena framework of soft-particles. The hard-sphere field studies frictionless amorphous packings inter-acting with hard-core normal forces only. Their results are a particular case ofthe more general problem of jammed matter composed of frictional and friction-less hard spheres as treated in the present work.3tatistical mechanics studies of jammed matter were initiated in 1989 withthe work of Sir Sam Edwards. It was first recognized that the main theoreti-cal difficulty to develop a statistical formulation results from the lack of well-defined conservation laws on which an ensemble description of jammed mattercould be based. Unlike equilibrium statistical mechanics where energy conser-vation holds, granular matter dissipates energy through frictional inter-particleforces. Therefore, it is doubtful that energy in granular matter could describethe microstates of the system and a new ensemble needs to be considered.Stemming from the fact that it is possible to explore different jammed con-figuration at a given volume in systematic experiments [18, 42, 20, 19], Edwardsproposed the volume ensemble (V-ensemble) as a replacement of the energy en-semble in equilibrium systems. A simple experiment is merely pouring grainsin a fixed volume and applying perturbations, sound waves or gentle tapping,to explore the configurational changes. Concomitant with a given network ge-ometry there is a distribution of contact forces or stresses in the particulatemedium. This means that the V-ensemble must be supplemented by the forceor stress ensemble (F-ensemble) determined by the contact forces [16, 43, 44]for a full characterization of jammed matter.Following this theoretical framework, a large body of experiments, theoryand simulations [13, 14, 16, 15, 21, 17, 18, 20, 19, 46, 47, 45, 37] have focused onthe study of volume fluctuations in granular media. Statistical studies have beenconcerned with testing for the existence of thermodynamics quantities such aseffective temperatures and compactivity as well as challenging the foundationsbased on the ergodic hypothesis or equal probability of the jammed states. It isrecognized that thermodynamic analogies may illuminate methods for attempt-ing to solve certain granular problems, but will inevitably fail at some point intheir application. The mode of this failure is an interesting phenomenon, as therange of phenomena for which ergodicity of some kind or other will apply or notis an interesting question. This has been illustrated by the compaction experi-ments of the groups of Chicago, Texas, Rennes and Schlumberger [18, 42, 20, 19].They have shown that reversible states exist along a branch of compaction curvewhere statistical mechanics is more likely to work. Conversely, experiments alsoshowed a branch of irreversibility where the statistical framework is not expectedto work. Poorly consolidated formations, such as a sandpile, are irreversible anda new “out-of-equilibrium” theory is required to describe them.In this paper we will describe the reversible states which are amenable tothe ergodic hypothesis while the unconsolidated irreversible states will be ad-dressed in a future work. Unfortunately, there is no first principle derivation ofthe granular statistical mechanics analogous to the Liouville theorem for equi-librium systems [48]. Therefore, advancing the statistical mechanics of granularmatter requires well-defined theoretical predictions that can be tested exper-imentally or numerically. While the possibility of a thermodynamic principledescribing jammed matter is recognized as a sensible line of research, the prob-lem with the statistical approach is that after almost 20 years there are nopractical applications yet. We require predictions of practical importance suchas equations of state relating the observables: pressure, volume, coordination4umber, entropy, etc, that may lead to new phenomena to be discovered andtested experimentally that will allow for a concrete exploitation of the thermo-dynamic framework. Here, we explore this problem by developing a systematic study of the V-ensemble of jammed matter. Our systems of interest are primarily packingsof granular materials, frictional colloids, infinitely rough grains and frictionlessdroplets mimicking concentrated emulsions. We attempt to answer the followingquestions using statistical mechanics methods:— What is a jammed state?— What are the proper state variables to describe the granular system atthe ensemble level?— What are the ensuing equations of state to relate the different observables?For instance, the simplest equilibrium thermal system, the ideal gas, is describedby pressure, volume and temperature ( p, V, T ) through pV = nRT . Is it possibleto identify the state variables for jammed matter which can describe the systemin a specified region of the phase space?— Which are the states with the maximum and minimum entropy (mostand least disordered)?— Can we define a ”compactimeter” to measure compactivity?— Can we provide a statistical interpretation of RLP and RCP, and use theensemble formalism to predict their density values?To answer these questions, we illuminate a diverse spectrum of data onsphere packings through a statistical theory of disordered jammed systems.Our results ultimately lead to a phase diagram providing a unifying view ofthe disordered sphere packing problem. An extensive program of computer sim-ulations of frictional and frictionless granular media tests the predictions andassumptions of the theory, finding good general agreement. The phase diagramintroduced here allows the understanding of how random packings fill space in3d for any friction, and eventually can be extended as a mathematical modelfor packings in 2d and in higher dimensions, anisotropic particles and other sys-tems. While physical application may be limited to 2d or 3d, higher dimensionsfind application in the science of information theory, where data transmissionand encryption find use for packing optimization. This paper is the second installment of a series of papers dealing with dif-ferent aspects of jammed matter from a statistics mechanics point of view. Ina companion paper, Jamming I [49], we describe the definition of the volumefunction based on the Voronoi volume of a particle and its relation to the coor-dination number.In the present paper, Jamming II, we describe calculations leading to thephase diagram of jammed matter and numerical simulations to test the theoret-ical results; a short version of the present paper has been recently published in550]. Jamming III [51] describes the entropy calculation to characterize random-ness in jammed matter and a discussion of the V-F ensemble. Jamming IV [52]studies the distribution of volumes from the mesoscopic ensemble point of viewand a generalized z-ensemble to understand the distribution of coordinationnumber. Jamming V [53] details the calculations for jamming in two dimen-sions at the mesoscopic level. The theory could be also extended to bidisperse[54] and high-dimensional systems [55].The outline of the present paper is as follows: Section 2 explain the basisof the Edwards statistical mechanics. Section 3 summarizes the results of Jam-ming I regarding the calculation of the volume function. Section 4 describesthe isostatic condition that defines the ensemble of jammed matter through theΘ jam function, and explains the difference between the geometrical, z , and me-chanical, Z , coordination number, which is important to define the canonicalpartition function. Section 5 defines the density of states in the partition func-tion. Section 6 explains the main results of this paper in the equations of state(21) and (23), and the phase diagram of Fig. 4. Finally, Section 7 explains thenumerical studies to test the theoretical predictions.
2. Statistical mechanics of jammed matter
Conventional statistical mechanics uses the ergodic hypothesis to derive themicrocanonical and canonical ensembles, based on the quantities conserved,typically the energy E [48]. Thus the entropy in the microcanonical ensembleis S ( E ) = k B log R δ ( E − H ( p, q )) dpdq , where H ( p, q ) is the Hamiltonian. Thisbecomes the canonical ensemble with exp[ −H ( ∂S/∂E )].The analogous development of a statistical mechanics of granular and otherjammed materials presents many difficulties. Firstly, the macroscopic size ofthe constitutive particles forbids the equilibrium thermalization of the system.Secondly, the fact that energy is constantly dissipated via frictional interparticleforces further renders the problem outside the realm of equilibrium statisticalmechanics due to the lack of energy conservation. In the absence of energyconservation laws, a new ensemble is needed in order to describe the systemproperties.Following this theoretical perspective, Edwards proposed the statistical me-chanics of jammed matter and interest in the problem of volume fluctuationshas flourished [14]. The central concept is that of a volume function W replacingthe role of the Hamiltonian in describing the microstates of the system in theV-ensemble and the stress boundaryΠ ij = Z σ ij dV with the stress σ ij = 1 / (2 V ) P c f ci r cj describing the F-ensemble [16, 43, 44],where f ci , r ci are the force and position at contact c . For simplicity only theisotropic case is described. Thus, only the pressure σ = σ ii / W i (for instance with a Voronoi tessellation as it will be done bellow),then the total volume, W , of a system of N particles is given by: W = N X i =1 W i . (1)The ensemble average of the volume function W provides the volume of thesystem, V = hWi , in an analogous way to the average of the Hamiltonian is theenergy in the canonical ensemble of equilibrium statistical mechanics.The full canonical partition function in the V-F ensemble is the startingpoint of the statistical analysis [16]: Z ( X, A ) = Z g ( W , Π) exp h − Π A − W X i Θ jam d W d Π , (2)where g ( W , Π) is the density of states for a given volume and boundary stress.Here Θ jam formally imposes the jamming restriction and therefore defines theensemble of jammed matter. This crucial function will be discussed at lengthbelow. As a minimum requirement it should ensure touching grains, and obe-dience to Newton’s force laws.Just as ∂E/∂S = T is the temperature in equilibrium systems, the temperature-like variables in granular systems are the compactivity [14]: X = ∂V∂S , (3)and the angoricity [from the Greek “ ´ αγχoζ ” (ankhos) = stress] [16] : A = ∂ Π ∂S . (4)The compactivity measures the power to compactify while the angoricitymeasures the power to stress. These quantities have remained quite abstractso far, perhaps, this fact being the primary reason for the great deal of contro-versy surrounding the statistical mechanics of jammed matter. In the presentpaper we attempt to provide a meaning and interpretation for the compactivity(the angoricity will be treated in subsequent papers) by developing equationsof states as well interpretations in terms of a “compactimeter” to measure the“temperature” of jammed matter. In Eq. (2), the analogue of the Boltzmannconstant are set to unity for simplicity, implying that the compactivity is mea-sured in units of volume and the angoricity in units of boundary stress (stresstimes volume). Note that in [16], the angoricity is denoted Z . Here we use A for the isotropic angoricitysince Z is used to denote the mechanical coordination number. In general the angoricity is atensor, A ij = ∂ Π ij ∂S , here we simplify to the isotropic case.
7n the limit of vanishing angoricity, A →
0, the system is described by theV-ensemble alone. This is the hard sphere limit which will be the focus of thepresent work where the following partition function describes the statistics: Z hardsph ( X ) = Z g ( W ) e −W /X Θ jam d W , (5)where g ( W ) is the density of states for a given volume W .
3. Volume Function
While it is always possible to measure the total volume of the system, it isunclear how to partition the space to associate a volume to each grain. Thus, thefirst step to study the V-ensemble is to find the volume function W i associatedto each particle that successfully tiles the system. This is analogous to theadditive property of energy in equilibrium statistical mechanics. This derivationis explained in Jamming I [49] and is the main theoretical result that leads to thesolution of the partition function for granular matter. We refer to this paper fordetails. Here we just state the main two results, Eqs. (6) and (7), and explaintheir significance for the solution of the statistical problem.Initial attempts at modelling W included, (a) a model volume function undermean-field approximation [14], (b) the work of Ball and Blumenfeld [56, 15] in2d and in 3d [57], and (c) simpler versions in terms of the first coordinationshell by Edwards [45]. These definitions are problematic: (a) is not given interms of the contact network, (c) is not additive, (b) and (c) are proportionalto the coordination number of the balls contrary to expectation (see below). InJamming I [49] we have found an analytical form of the volume function in anydimensions and demonstrated that it is the Voronoi volume of the particle.The definition of a Voronoi cell is a convex polygon whose interior consistsof all points closer to a given particle than to any other (see Fig. 1). Its formulain terms of particle positions for monodisperse spherical packings in 3d is [49]: W vor i = 13 *(cid:18) R min j r ij cos θ ij (cid:19) + s ≡ hW si i s , (6)where ~r ij is the vector from the position of particle i to that of particle j , theaverage is over all the directions ˆ s forming an angle θ ij with ~r ij as in Fig. 1,and R is the radius of the grain. W si defines the orientational Voronoi volumewhich is obtained without the integration over ˆ s . This formula has a simpleinterpretation depicted in Fig. 1.The Voronoi construction is additive and successfully tiles the total volume.Prior to Eq. (6), there was no analytical formula to calculate the Voronoi volumein terms of the contact network r ij . Equation (6) provides such a formula whichallows theoretical analysis in the V-ensemble. However this microscopic versionis difficult to incorporate into the partition function since it would necessitate afield theory. The next step is to develop a theory of volume fluctuations to coarse8rain W vori over a mesoscopic length scale and calculate an average volumefunction. The coarsening reduces the degrees of freedom to one variable, thecoordination number of each grain, and defines a mesoscopic volume functionwhich is more amenable to statistical calculations than Eq. (6). We find a(reduced) free volume function [49]: w ( z ) ≡ hW si i i − V g V g ≈ √ z , (7)approximately valid for monodisperse hard spheres with volume V g , where z isthe geometrical coordination number (which is different from the mechanicalcoordination number, Z , see below). The average is over all the balls.The inverse relation with z in Eq. (7) is in general agreement with ex-periments [46]. Many experiments have focused on the analysis of the systemvolume and single particle volume fluctuations in granular media [19, 46, 47, 18].Of particular importance to the present theory are the recent advances in X-ray tomography [46] and confocal microscopy [37, 58] which have revealed thedetailed internal structure of jammed matter allowing for the study of the freevolume per particle. By partitioning the space with Voronoi diagrams, thesestudies show that W ivor is distributed with wide tails [46, 47, 45, 37]. More im-portantly, the X-ray tomography experiments with 300,000 monodisperse hardspheres performed by the Aste group [46] find that the Voronoi volumes areinversely proportional, on average, to the coordination number of the particle, z , in reasonable agreement with the prediction of the volume function Eq. (7).Such data is displayed in Fig. 6 in [46] where the local volume fraction definedas φ − i − W ivor /V g is plotted against the coordination number. From Eq.(7) we find: φ i = 1 w + 1 = zz + 2 √ , (8)which agrees relatively well with the shape of the curves displayed in Fig. 6 of[46] for different packing preparations. It should be noted that for a more pre-cise comparison a coarse grained volume fraction and geometrical coordinationnumber should be considered in Eq. (8). A detailed discussion is provided inSection 7.6.The free volume function decreases with z as expected since the more con-tacts per grain, the more jammed the particle is and the smaller the free volumeassociated with the grain. The coordination number z in Eq. (7) can be con-sidered as a coarse-grained average associated with “quasiparticles” with meso-scopic free volume w . By analogy to the theory of quantum energy spectra,we regard each state described by Eq. (7) as an assembly of “elementary exci-tations” which behave as independent quasiparticles. As a grain is jammed inthe packing, it interacts with other grains. The role of this interaction betweengrains is assumed in the calculation of the volume function (7) and it is implicitin the coarse-graining procedure explained above. The quasiparticles can beconsidered as particles in a self-consistent field of surrounding jammed matter.9n the presence of this field, the volume of the quasiparticles depends on thesurrounding particles as expressed in Eq. (6). The assembly of quasiparticlescan be regarded as a set of non-interacting particles (when the number of ele-mentary excitations is sufficiently low) and a single particle approximation canbe used to solve the partition function as we will formulate in Section 6.While Eq. (6) is difficult to treat analytically, the advantage of the meso-scopic Eq. (7) is that the partition function can be solved analytically since w depends on z only, instead of r ij . The key result is the relation between theVoronoi volume and the coordination number which allow us to incorporate thevolume function into a statistical mechanics approach of jammed hard spheres,by using the constraint of mechanical stability as we show below.
4. Definition of jamming: isostatic conjecture
The definition of the constraint function Θ jam is intimately related to theproper definition of a jammed state and should contain the minimum require-ment of mechanical equilibrium.Distinguishing between metastable and mechanically stable packings thatdefine the jammed state through the Θ jam function remains a problem underdebate, related to the more fundamental question of whether or not a jammedpacking is well-defined [10]. In practice, it is widely believed that the isostaticcondition is necessary for a jammed disordered packing following the Alexanderconjecture [59, 60, 61] which was tested in several works [27, 29, 30, 33]It is well known that mechanical equilibrium imposes an average coordinationnumber larger or equal than a minimum coordination where the number of forcevariables equals the number of force and torque balance equations [59, 61, 60].The so-called isostatic condition.In the case of frictionless spherical particles the isostatic condition is Z =2 d = 6 (in 3d), while the coordination in the case of infinitely rough particles(with interparticle friction coefficient µ → ∞ ) is Z = d + 1 = 4, where d is the dimension of the system. Numerical simulations and theoretical worksuggest that at the jamming transition the system becomes exactly isostatic[29, 31, 27, 34, 61, 33, 60]. But no rigorous proof of this statement exists. In thefollowing, we will use this isostatic conjecture to define the ensemble of jammedmatter.While isostaticity holds for perfectly smooth and infinitely rough grains, themain problem is to extend it to finite frictional systems. For finite friction,the Coulomb condition takes the form of an inequality between the normalinterparticle contact force F n and the tangential one F t : F t ≤ µF n and thereforeno trivial solution to the minimum number of contacts can be obtained.The problem can be understood as an optimization of an outcome basedon a set of constraints, i.e., minimizing a Hamiltonian of a system over a con-vex polyhedron specified by linear and non-negativity constraints. The isostaticcondition can be augmented to indicate the number of extra equations for con-tacts satisfying the Coulomb condition, which is analogous to the number of10 able 1: Number of constraints and variables determining the isostatic condition for differentsystems of spherical particles. Friction N n N t E f E t µ = 0 N Z dN µ finite N Z ( d − N Zf ( µ ) dN d ( d − N f ( µ ) µ = ∞ N Z ( d − N Z dN d ( d − N redundant constraints in Maxwell constraint counting of rigidity percolation.This suggests that, for a finite value of µ , the original nonlinear problem canbe mapped to a linear equation problem if we know how many extra equationsshould be added. This problem is treated with more details in a future paper.In the following, we present a rough approximation suggesting that the co-ordination number could be related to friction. We consider the following argu-ment. Consider a set of spherical particles interacting via normal and tangentialcontact forces. These can be the standard Hertz and Mindlin/Coulomb forces ofcontact mechanics [35, 36, 32], respectively, see Section 7.1. We set N : numberof particles, N n : number of unknown normal forces, N t : number of unknowntangential forces, E f : number of force balance equations, E t : number of torquebalance equations, Z = 2 M/N : average coordination number of the packing,where M is the total number of contacts, f ( µ ): undetermined function of thefriction coefficient µ such that 1 − f ( µ ) is the fraction of spheres that can ro-tate freely ( f (0) = 0 and f ( ∞ ) = 1), and f ( µ ): undetermined function of µ indicating the ratio of contacts satisfying F t < µF n , which satisfies f (0) = 0and f ( ∞ ) = 1.A packing is isostatic when N n + N t = E f + E t . (9)The average coordination number at the isostatic point is then (see Table 1): Z ( µ ) = 2 d / d − f ( µ )1 + ( d − f ( µ ) , (10)reducing to the known Z = 2 d for frictionless particles and Z = d + 1 forinfinitely rough particles.In what follows, we use the numerical fact that interpolating between the twoisostatic limits, there exist packings of finite µ with the coordination numbersmoothly varying between Z ( µ = 0) = 6 and Z ( µ → ∞ ) → Z vs µ dependence is independent of the preparationprotocol as obtained in our simulations (see Fig. 2). This result then generalizesthe isostatic conditions from µ = 0 and µ = ∞ to finite µ .It should be noted that we cannot rule out that other preparation protocolscould give rise to other dependence of Z on µ . However, we will see that the11btained phase diagram is given in terms of Z ; the main prediction of the theorywould still be valid irrespective of the particular dependence Z ( µ ). That is, thetheory does not assume anything about the relation between the interparticlefriction and Z .It is worth noting that other attempts to define the jammed state have beendeveloped. A rigorous attempt is that of Torquato et al. who propose threecategories of jamming [62]: locally jammed, collectively jammed and strictlyjammed based on geometrical constraints. However, the definition of [62] isbased purely on geometrical considerations and therefore only valid for friction-less particles. Thus it is not suitable for granular materials with inter-particlefrictional tangential forces; their configurational space is influenced by the me-chanics of normal and shear forces. Other approaches to define a jammed statebased on the potential energy landscape [29] fail for granular materials too sincesuch a potential does not exist for frictional grains due to their inherent path-dependency. Thus, the definition of jammed state for granular materials mustconsider interparticle normal and tangential contact forces beyond geometry. Inthe companion paper [51] we elaborate on this problem. We have acknowledged a difference between the geometrical coordinationnumber z in Eq. (7) and the mechanical coordination number Z which countsonly the contacts with non-zero forces. Below we discuss the bounds of z .Since some geometrical contacts may carry no force, then we have: Z ≤ z. (11)To show this, imagine a packing of infinitely rough ( µ → ∞ ) spheres withvolume fraction close to 0 .
64. There must be z = 6 nearest neighbors aroundeach particle on the average. However, the mechanical balance law requires only Z = 4 contacts per particle on average, implying that 2 contacts have zero forceand do not contribute to the contact force network.Such a situation is possible as shown in Fig. 3: starting with the contactnetwork of an isostatic packing of frictionless spheres having Z = 6 and all con-tacts carrying forces (then z = 6 also as shown in Fig. 3a), we simply allow theexistence of tangential forces between the particles and switch the friction coef-ficient to infinity. Subsequently, we solve the force and torque balance equationsagain for this modified packing of infinitely rough spheres but same geometricalnetwork, as shown in Fig. 3b [Notice that the shear force is composed of anelastic Mindlin component plus the Coulomb condition determined by µ , seeSection 7.1 for details. Thus when µ → ∞ , the elastic Mindlin component stillremains].The resulting packing is mechanically stable and is obtained by setting tozero the forces of two contacts per ball, on average, to satisfy the new force andtorque balance condition for the additional tangential force at the contact. Sucha solution is guaranteed to exist due to the isostatic condition: at Z = 4 thenumber of equations equals the number of force variables. Despite mechanical12quilibrium, giving Z = 4, there are still z = 6 geometrical contacts contributingto the volume function.Therefore, we identify two types of coordination number: the geometricalcoordination number, z , contributing to the volume function and the mechanicalcoordination number, Z , measuring the contacts that carry forces only. Thisdistinction is crucial to understand the sum over the states and the bounds inthe partition function as explained below.These ideas are corroborated by numerical simulations in Section 7.6, below.The packings along the vertical RCP line found in the simulations (see Fig. 11)have approximately the same geometrical coordination number, z ≈
6. However,they differ in mechanical coordination number, going from the frictionless point Z = 6 to Z ≈ N d positions of the particles areconstrained by the
N z/ | ~r ij | = 2 R , of rigidity. Thus,the number of contacts satisfies N z/ ≤ N d , and z is bounded by: z ≤ d. (12)Notice that this upper bound applies to the geometrical coordination, z andnot to the mechanical one, Z , and it is valid for any system irrespective of thefriction coefficient, from µ = 0 → ∞ .Furthermore, in relation with the discussion of frictionless isostaticity, it isbelieved that above 2 d , the system is partially crystallized. To increase thecoordination number above 6, it is necessary to create partial crystallization inthe packing, up to the point of full order of the FCC lattice with coordinationnumber 12. Thus, by defining the upper bound at the frictionless isostatic limitwe also exclude from the ensemble the partially crystalline packings. This is animportant point, akin to mathematical tricks employed in replica approaches toglasses [38].In conclusion, the mechanical coordination number, Z , ranges from 4 to 6as a function of µ , and provides a lower bound to the geometrical coordinationnumber, while the upper bound is 2 d . A granular system is specified by theinterparticle friction which determines the average mechanical coordination atwhich the system is equilibrated, Z ( µ ). The possible microstates in the ensembleavailable for this system follow a Boltzmann distribution Eq. (5) for statessatisfying the following bounds: Z ( µ ) ≤ z ≤ d = 6 . (13)
5. Density of states
According to the statistical mechanics of jammed matter, the volume par-tition function Z is defined by Eq. (5). In the quasiparticle approximation wecan write: 13 hardsph ( X ) = Z . . . Z g (cid:16) X w ( z i ) (cid:17) e − P w ( z i ) /X ×× Θ jam N Y i dw ( z i ) . (14)Considering N non-interacting quasiparticles with free volume w ( z ), thepartition function can be written as: Z hardsph = (cid:16) Z g ( w ) e − w/X Θ jam dw (cid:17) N . (15)Here, g ( w ) is the density of states for a given quasiparticle free volume.Since the mesoscopic w is directly related to z through Eq. (7), we changevariables to the geometrical coordination number in the partition function. Thedensity of states for a single quasiparticle, g ( w ), is: g ( w ) = Z Z P ( w | z ) g ( z ) dz, (16)where P ( w | z ) is the conditional probability of a free volume w for a given z ,and g ( z ) is the density of states for a given z . Here, we have used the boundsin Eq. (13).The next step in the derivation is the calculation of g ( z ) which is developedin three steps:First, we consider that the packing of hard spheres is in a jammed configura-tion where there can be no collective motion of any contacting subset of particlesleading to unjamming when including the normal and tangential forces betweenthe particles. This definition is an extended version of the collectively jammedcategory proposed by Torquato et al. [62] that goes beyond the merely locallyjammed configuration of packings, unstable to the motion of a single particle.While the degrees of freedom are continuous, the fact that the packing is col-lectively jammed implies that the jammed configurations in the volume spaceare not continuous. Otherwise there would be a continuous transformation inthe position space that would unjam the system contradicting the fact that thepacking is collectively jammed.Thus, we consider that the configuration space of jammed matter is discretesince we cannot change one configuration to another in a continuous way. Similarconsideration of discreteness has been studies in [62]. Notice that the volumelandscape could be continuous since we could change the volume as well, or,in the case of soft particle, we can deform them. Additionally, in the case offrictionless packings of soft particle, the energy of deformation is well-defined,and we can define the collectively jammed configuration as a minimum in theenergy with definitive positive Hessian (or a zero order saddle) like in [29]. Inthis case, there could be no continuous transformation of the particle coordinates14hat brings one jammed state to the next, unless we deform the particles. Thus,the space is discrete in this case as well.Second, we call the dimension per particle of the configuration space as D and consider that the distance between two jammed configurations is notbroadly distributed (meaning that the average distance is well-defined). We callthe typical (average) distance between configurations in the volume space as h z , and therefore the number of configurations per particle is proportional to1 / ( h z ) D . The constant h z plays the role of the Planck’s constant, h , in quantummechanics.Third, we add z constraints per particle due to the fact that the particle isjammed by z contacts. Thus, there are N z position constrains ( | ~r ij | = 2 R )for a jammed state of hard spheres as compared to the unjammed “gas” state.Therefore, the number of degrees of freedom is reduced to D− z , and the numberof configurations is then 1 / ( h z ) D− z . Since the term 1 / ( h z ) D is a constant, it willnot influence the average of the observables in the partition function (althoughit changes the value of the entropy, see Jamming III [51]). Therefore, the densityof states g ( z ) is assumed to have an exponential form: g ( z ) = ( h z ) z = e − z/z c , (17)with z − c ≡ ln(1 /h z ). Physically, we expect h z ≪
1, then g ( z ) is an exponen-tially rapid decreasing function with z . The exact value of h z can be determinedfrom a fitting of the theoretical values to the simulation data. The most popu-lated state is the highest volume at z = 4 while the least populated state is theground state at z = 6. A constant could be added in (17) but it has no conse-quence for the average of the observables in the partition function. However, itaffect the value of the entropy. For the calculation of the entropy we considerthat there is a single mesoscopic ground state and use g ( z ) = e − ( z − d ) /z c [51].We see that the negative of the geometrical coordination number, − z , playsthe role of a number of degrees of freedom for a packing, due to the extraposition constrains of the contacting particles. The coordination z can be thenconsidered as the number of degrees of “frozen” per particle. Another way tounderstand Eq. (17) is the following. In the case of a continuous phase space ofconfigurations we would obtain g ( z +1) g ( z ) = 0. However, since the space of volumeconfigurations is discrete as discussed above, the ratio g ( z +1) g ( z ) ∼ h z . This impliesagain Eq. (17).The conditional probability P ( w | z ) depends on the w function, w = √ z .The average is taken over a certain mesoscopic length scale since the volume ofa particle depends on the positions of the particles surrounding it. Practically,such length scale is approximately of several particle diameters. w is a coarse-grained volume and independent of the microscopic partition of the particles,implying: P ( w | z ) = δ ( w − √ /z ) . (18)The meaning of Eq. (18) is that we neglect the fluctuations in the coordinationnumber due to the coarse graining procedure. A more general ensemble can be15onsidered where the fluctuations in the geometrical coordination number aretaken into account. We have extended our calculations to consider fluctuationsin z and find a similar phase diagram as predicted by the present partition func-tion. This generalized z-ensemble will be treated in a future paper, where theboundaries of the phase diagram are regarded as a second-order phase transi-tion. The generalized z-ensemble allows for the calculation of the probabilitydistribution of the coordination number, beyond the assumption of the deltafunction distribution in Eq. (18). The resulting prediction of P ( z ) is in generalagreement with simulations (see [52]).Substituting Eq. (18) and Eq. (17) into Eq. (15), we find the isostaticpartition function which is used in the remaining of this study: Z iso ( X, Z ) = Z Z ( h z ) z exp − √ zX ! dz. (19)
6. Phase diagram
Next, we obtain the equations of state to define the phase diagram of jammedmatter by solving the partition function. From Eq. (19), we calculate theensemble average volume fraction φ = ( w + 1) − = z/ ( z + 2 √
3) as: φ ( X, Z ) =1 Z iso ( X, Z ) Z Z zz + 2 √ − √ zX + z ln h z ! dz. (20)We start by investigating the limiting cases of zero and infinite compactivity.(a) In the limit of vanishing compactivity ( X → z = 6 contributes to the partition function. Then weobtain the ground state of jammed matter with a density: φ RCP = φ ( X = 0 , Z ) = 66 + 2 √ ≈ . , Z ( µ ) ∈ [4 , . (21)The meaning of the subscript RCP in (21) will become clear below.(b) In the limit of infinite compactivity ( X → ∞ ), the Boltzmann factorexp[ − √ / ( zX )] →
1, and the average in (20) is taken over all the states withequal probability. We obtain: φ RLP ( Z ) = φ ( X → ∞ , Z ) == 1 Z iso ( ∞ , Z ) Z Z zz + 2 √ z ln h z ) dz. (22)The constant h z determines the minimum volume in the phase space. Weexpect h z ≪
1, such that the exponential in Eq. (22) decays rapidly. Then16he leading contribution to Eq. (22) is from the highest volume at z = Z andtherefore: φ RLP ( Z ) ≈ ZZ + 2 √ , Z ( µ ) ∈ [4 , . (23)This dependence of the volume fraction on Z suggests using the ( φ, Z ) planeto define the phase diagram of jammed matter as plotted in Fig. 4. Theequations of state (21) and (23) are plotted in the ( φ, Z ) plane in Fig. 4 providingtwo limits of the phase diagram. Since the mechanical coordination number islimited by 4 ≤ Z ≤ Z = 4 for infinitely roughgrains, denoted the granular-line or G-line in Fig. 4.All mechanically stable disordered jammed packings lie within the confininglimits of the phase diagram (indicated by the yellow zone in Fig. 4), while thegrey shaded area in Fig. 4 indicates the forbidden zone. For example, a packingof frictional hard spheres with Z = 5 (corresponding to a granular materialwith interparticle friction µ ≈ . φ < φ RLP ( Z = 5) = 5 / (5 + 2 √
3) = 0 .
591 or above φ > φ
RCP = 0 . (i) The RCP limit.–
Stemming from the statistical mechanics approach,the RCP limit arises as the result of the relation (21), which gives the maximumvolume fraction of disordered packings under the mesoscopic framework. To theright of the RCP-line, packings exist only with some degree of order (for instancewith crystalline regions). The prediction, φ RCP = 66 + 2 √ ≈ . , (24)is valid for all friction coefficients and approximates the experimental and nu-merical estimations [5, 6, 7, 8] which find a close packing limit independent offriction in a narrow range around 0.64.Beyond the fact that 63-64% is commonly quoted as RCP for monodispersehard spheres, we present a physical interpretation of that value as the groundstate of frictional hard spheres characterized by a given interparticle frictioncoefficient. In this representation, as µ varies from 0 to ∞ and Z decreases from6 to 4, the state of RCP changes accordingly while its volume fraction remainsthe same, given by Eq. (24). The present approach has led to an unexpectednumber of states that all lie in the RCP line from Z = 6 to Z = 4 as depictedin Fig. 4, suggesting that RCP is not a unique point in the phase diagram.An important prediction is that for frictionless systems there is only onepossible state at Z = 6. It is important to note that there is one state only atthe mesoscopic level used in the theory. However, for a single mesoscopic state,we expect many microstates, which are averaged out in the mesoscopic theory17f the volume function. Thus, there could be more jammed states surroundingthe frictionless point in the phase diagram. However, we expect that thesestates are clustered in a narrow region around the frictionless point. To accessthese microstates it would require different preparation protocols, analogous tothe dependence of the glass transition temperature on cooling rates in glasses[40, 38].It is interesting to note that replica approaches to the jammed states [38] pre-dict many jammed isostatic states with different volume fractions for frictionlesshard spheres. The first question is whether these states in the force/energy en-semble are the same as the volume ensemble states that we treat in our approach.It should be noted that jamming in [40, 38] is obtained when the particles arerattling infinitely fast in their cages, in the limit of infinite number of collisionsper unit time. That is when the dynamic pressure related to the momentumof the particles diverges. Second, from our point of view, an investigation ofthe fluctuations of the microstates as well as more general ensembles allowingfor fluctuations in the coordination number could be considered. These stud-ies, may reveal whether the frictionless point is unique or not. This point isinvestigated further in Section 7.5. (ii) The RLP limit.—
Equation of state (23) provides the lowest volumefraction for a given Z and represents a statistical interpretation of the RLPlimit depicted by the RLP-line in Fig. 4. We predict that to the left of this line,packings are not mechanically stable or they are experimentally irreversible asdiscussed in [18, 19, 20].A review of the literature indicates that there is no general consensus on thevalue of RLP as different estimations have been reported ranging from 0.55 to0.60 [6, 9, 8], proposing that there is no clear definition of RLP limit. The phasediagram proposes a solution to this problem. Following the infinite compactivityRLP-line, the volume fraction of the RLP decreases with increasing frictionfrom the frictionless point ( φ, Z ) = (0 . , Z →
4. Indeed, experiments [6] indicate that lower volumefractions are achieved for larger coefficient of friction. We predict the lowestvolume fraction in the limit: µ → ∞ , X → ∞ and Z → h z →
0) at φ minRLP = 44 + 2 √ ≈ . . (25)Even though this is a theoretical limit, our results indicate that for µ > . ± . φ, Z ( µ )) plane: 18a) The frictionless point µ = 0, denoted J-point in [29], at J ≡ ( φ RCP , Z (0)) = (0 . , , corresponds to a system of compressed emulsions in the limit of small osmoticpressure as measured by J. Bruji´c [63].(b) The lowest coordination number Z = 4 plotted as the G-line defines twoassociated points from the lowest volume fraction of loose packings at infinitecompactivity, L-point, L ≡ ( φ minRLP , Z ( ∞ )) = (0 . , , to the zero compactivity state of close packing, C-point, C ≡ ( φ RCP , Z ( ∞ )) = (0 . , . The full JCL triangle defines the isostatic plane where the frictional hard spherepackings reside. (iii)
Intermediate isocompactivity states.—
For finite X , Eq. (20) canbe solved numerically. For each X , the function φ ( X, Z ) can be obtained and isplotted as each isocompactivity color line in Fig. 4. Between the two limits Eqs.(21) and (23), there are packings inside the yellow zone in Fig. 4 with finitecompactivity, 0 < X < ∞ . Since X controls the probability of each state, like incondensed matter through a Boltzmann-like factor in Eq. (5), it characterizesthe number of possible ways to rearrange a packing having a given volume andentropy, S . Thus, the limit of the most compact and least compact stablearrangements correspond to X → X → ∞ , respectively. Between theselimits, the compactivity determines the volume fraction from RCP to RLP. Dependence on h z and negative compactivity.– Of interest is the dependence of our prediction on the “Planck constant” ofjammed matter, h z , that determines the minimum volume in the phase space.The equation of state (23) and the prediction of the minimum RLP, Eq. (25)have been obtained by considering h z → X →∞ , the only state contributing to the volume partition function is the mostpopulated at z = Z ( µ ).The approximation h z ≪ h z = e − in order to fit the theoretical values with the numerical ones for finitecompactivity. This extremely small constant shifts the value of the minimumRLP a little bit to the right of the phase diagram from the prediction of Eq.(25). Indeed, when we plot the phase diagram of Fig. 4 with a h z = e − aslightly larger value of φ minRLP is obtained as seen in Fig. 4. In the unphysicallimit of h z → φ minRLP + φ RCP ) / h z .It is interesting to note that we can extend the compactivity to negativevalues and study the range X : 0 − → −∞ . Indeed, for any value of h z , inthe limit X → − we obtain the lowest volume fraction of the prediction of Eq.(23). Thus, properly speaking, the minimum value of RLP is defined in the limitof negative zero compactivity, and this limit is independent of the value of h z .The entropy has an interesting behavior in this regime which will be discussedlater.It is important to realize that the region from X : 0 − → ∞ disappears when h z →
0, thus the only meaningful limit is that of X → + ∞ (which equals boththe limits X → −∞ and X → − when h z → h z .Since we always expect h z ≪
1, it may not be necessary to use the negativecompactivity states to describe RLP. We leave this interesting observation forfurther investigations where the entropy of jamming is treated in more detailsin [51].
Equation (7) plays the role of the Hamiltonian of the jammed system and thejammed configurations can be considered as the minima of (7). Inspired by thephysics of glasses and supercooled liquids, we imagine a ”volume landscape”,analogous to the energy landscape in glasses [64], as a pictorial representationto describe the states of jammed matter. Each mesoscopic jammed state (de-termined by the positions of the particles, denoted ~r i , and its correspondingvolume w ) is depicted as a point in Fig. 5. At the mesoscopic level, the vol-ume landscape has different levels of constant z , analogous to energy levels inHamiltonian systems.The lowest volume corresponds to the FCC/HCP structure (with kissingnumber z = 12), as conjectured by Kepler [1]. Other lattice packings, such asthe cubic lattice and tetrahedron lattice, have higher volume levels in this rep-resentation. Beyond these ordered states, the ensemble of disordered packingsis identified within the yellow area in Fig. 5a, corresponding to a system withinfinite friction. In this case the partition function is integrated from z = 4 to6, thus all the states are sampled in the configuration average as indicated bythe arrow in Fig. 5a. When X → ∞ , all the states are sampled with equalprobability, and, when X →
0, the ground state is the most probable. As thecompactivity is varied, the states along the G-line in the phase diagram result.For µ → ∞ , the maximum volume, w ( z = 4) = √ / z = 4 when µ → ∞ in analogy with the high energy states in a classical system.As we set the friction coefficient to a finite value, the available states in thevolume are less, since the integration in the partition function is in the region20 ( µ ) ≤ z ≤
6. This is indicated in Fig. 5b for a generic Z ( µ ). Finally when thefriction vanishes, we obtain only the ground state z = 6 as indicated in Fig. 5c.For any value of µ , the lowest state is always at z = 6. This correspondsto the states exemplified by the RCP-line in the phase diagram, all of themwith a geometrical coordination z = 6. Equation (21) indicates that the RCPcorresponds to the ground state of disordered jammed matter for a given frictionwhich determines Z , while the RLP states are achieved for higher volume levelsas indicated in Fig. 5. An important conclusion is the following: The statesalong the RCP line all have z = 6, independent of Z (and h z ) which ranges from4 to 6 as a function of the friction coefficient µ . The states along the RLP lineall have z = Z (when h z → ~r i ,parameterized by a common value of z with a density of states g ( z ). The basinsrepresent single states only at the mesoscopic level providing a mesoscopic viewof the landscape of jammed states. This is an important distinction arising fromthe fact that the states defined by w ( z ) = 2 √ /z , Eq. (7), are coarse grainedfrom the microscopic states defined by the microscopic Voronoi volume Eq. (6)in the mesoscopic calculations leading to (7) as discussed in Jamming I [49].This fact has important implications for the present predictions which will bediscussed in Section 7.5. The advantage of the volume landscape picture is thatit allows visualization of the corresponding average over configurations that giverise to the macroscopic observables of the jammed states. Further statistical characterization of the jammed structures can be obtainedthrough the calculation of the equations of state in the three-dimensional space(
X, φ, S ), with S the entropy, as seen in Fig. 6.21he entropy density, s = SN , is obtained as: s ( X, Z ) = h w i /X + ln Z iso ( X, Z ) (26)This equation is obtained in analogy with equilibrium statistical mechanicsand it is analogous to the definition of free energy: F = E − T S where F = − T ln Z is the free energy. We replace T → X , E → h w i . Therefore, F = E − T S or S = ( E − F ) /T = E/T + ln Z is now s ( X, Z ) = h w i /X + ln Z iso ( X, Z ), whichis plotted as the equation of state in Fig. 6.Each curve in the figure corresponds to a system with a different Z ( µ ). Theprojections S ( X ) and S ( φ ) in Fig. 6 characterize the nature of randomnessin the packings. When comparing all the packings, the maximum entropy isat φ minRLP and X → ∞ while the entropy is minimum for φ RCP at X → < Z ( µ ) < J -point at Z = 6 the entropy vanishes. Moreprecisely, it vanishes for a slightly smaller φ than φ RCP of the order h z . Strictlyspeaking, the entropy diverges to −∞ at φ RCP as S → ln X for any value of Z , in analogy with the classical equation of state, when we approach RCP todistances smaller than h z . However, this is an unphysical limit, as it would belike considering distances in phase space smaller than the Planck constant.It is commonly believed that the RCP limit corresponds to a state with thehighest number of configurations and therefore the highest entropy. However,here we show that the states with a higher compactivity have a higher entropy,corresponding to looser packings. Within a statistical mechanics frameworkof jammed matter, this result is a natural consequence and gives support tosuch an underlying statistical picture. A more detailed study of the entropy isperformed in Jamming III [51].The interpretation of the RCP as the ground state, X →
0, with vanishingentropy ( S → φ FCC = 0 . X , being maximum for the RLP limit. We notethat the microscopic states contribute still to the entropy of RCP giving rise tomore states than predicted by the present mesoscopic approach. This case isdiscussed in more detail in [51].The equation of state φ ( X ) for different values of Z can be seen in the pro-22ection of Fig. 6. The volume fraction diminishes with increasing compactivityaccording to the theoretical picture of the phase diagram. The curves φ ( X )qualitatively resemble the reversible branch of compaction curves in the experi-ments of [18] for shaken granular materials and oscillatory compression of grains[20] suggesting a correspondence between X and shaking amplitude. The in-tention is that, different control parameters in experiments could be related toa state variable, and therefore might help experimentalists to describe resultsobtained under different protocols.For any value of Z , there is a common limit φ → φ RCP as X →
0, indicatingthe constant volume fraction for all the RCP states. The singular nature of thefrictionless J -point is revealed as the volume fraction remains constant for anyvalue of X , explaining why this point is the confluence of the isocompactivitylines, including RCP and RLP. We conclude that at the frictionless J-point thecompactivity does not play a role, at least at the mesoscopic level. A reanalysis of the available experimental data tends to agree with the abovetheoretical predictions. However, it would be desirable to perform more con-trolled experiments in light of the present results. As in other out of equilibriumsystems, such as glasses, the inherent path-dependency of jammed matter mate-rializes in the fact that different packing structures can be realized with differentpreparation protocols [17, 18, 19, 20, 42] involving tapping, fluidized beds, set-tling particles at different speeds, acoustic perturbations or pressure waves [9].Due to this reason, the value of φ RCP has not been determined yet for monodis-perse frictionless systems. The more extensive evidence for a frictionless RCPappears from simulations, which has been extensively performed for hard andsoft sphere systems. They find a common value of RCP for many differentpreparation protocols [27, 29, 40]. For frictional materials, experiments of Scottand Kilgour [6], and others display a nearly universal value of volume fractionfor RCP consistent with the theoretical estimation.Previous experiments [9] and simulations [32] find that lower volume frac-tions can be achieved for smaller settling speeds of the grains or slower com-pression (or quenching) rates during packing preparation. Colloidal glasses finda similar scenario [38, 41], with their glass transition temperature dependenceon the quench rate of cooling, although due to different reasons. This raises thequestion of whether the jamming point is unique or determined by the prepa-ration protocol [41].The experiments performed by Onoda and Liniger with glass spheres in liq-uid of varying density to adjust conditions of buoyancy in the limit of vanishinggravity show that for larger values of gravity, the volume fraction decreases.This indicates that settling speed of the particles can determine the final vol-ume fraction. Indeed, it was found numerically that the compression rate duringpreparation of the packings is a systematic way to obtain lower packing fractions,such that lower volume fractions can be achieved with quasi-static compressionrates during preparation of the packings [32].23n the other hand, measuring the mechanical coordination number in ex-periments seems to be an even more difficult task. In the phase diagram, theRCP can be found at the frictionless point with the highest possible coordina-tion number Z = 6. Such a value of Z has been observed in the experiments ofBruji´c et al. on concentrated emulsions near the jamming transition [63] andsimulations of droplets.The experiments of Bernal [4] simultaneously measured the coordinationnumber and the volume fraction. Indeed, Mason, a postgraduate student ofBernal, took on the task of shaking glass balls in a sack and freezing the packingby pouring wax over the whole system. He would then carefully take the packingapart, ball by ball, painstakingly recording the positions of contacts for each ofover 8000 particles [4]. Their result is a coordination number Z ≈ φ ≈ .
64 which corresponds to the prediction of the frictionlesspoint in the phase diagram. This may indicate that the particles used in theexperiment of Bernal have a low friction coefficient.Another explanation, more plausible, is that the coordination number mea-sured by Bernal is not the mechanical one, but the geometrical one. Indeed, wemay argue that the only coordination that can be measured in experiments ofcounting balls ’a la Bernal’ is the geometrical one, since, one is never sure if thecontacting balls were carrying a force or not. Bernal could only measure geome-tries, not forces. This is in addition to the uncertainty in the determination ofthe contacts with such a rudimentary method as pouring wax and then countingone by one the area of contact. The theory predicts that along the RCP-line,all the RCP states have geometrical z = 6, while Z ranges from 4 to 6. Thus,it is quite plausible that the coordination number in Bernal experiments is thegeometrical one, z = 6, in agreement with the theory.It is worth noting that since the labor-intensive method patented half acentury ago by Bernal, other groups have extracted data at the level of theconstituent particles using X-ray tomography [46]. Such experiments may givea clue to the relationship between coordination number and volume fraction.The experimental data, to date, seems in good agreement with the presenttheory. However, this method does not directly determine the contacts to verifywhether the particles are touching or just very close together. Furthermore,the method cannot measure forces so the distinction between geometrical andmechanical coordination is not possible to achieve in this experiment.An answer might be obtained using methods from biochemistry being devel-oped at the moment [68]. These methods promise to provide the high resolutionto determine the contact network with accuracy to develop an experimental un-derstanding of the problem.An alternate way might be generating the packings in the phase spacethrough numerical simulations, where both volume fraction and coordinationnumber could be easily determined. Since Z is directly determined by µ , andthe compactivity determines what value of the volume fraction a packing hasbetween the limits of the phase diagram, the main question is how to generatepackings with different φ for a fixed µ to allow the exploration of the phasediagram. 24 . Numerical tests In this section we perform two different numerical tests: on the predictionsof the theory and the assumptions of the theory. The former are explained inSection 7.3 while the latter are elaborated in Section 7.6. It is worth mentioningthat while the predictions can be tested with packings prepared numerically andexperimentally, the test of the assumptions of the theory is not so trivial. Thisis because the theory is based on the existence of quasiparticles which carry theinformation at a mesoscopic scale. Thus, in principle we cannot use real pack-ings, such as computer generated packings or experimental ones, to measure thequasiparticles. The information obtained from those packings already containsensemble averages through the Boltzmann distribution and density of states.That is, it is in principle not possible to isolate the behavior of quasiparticlesfrom the real measures rendering difficult to properly test some of the assump-tions of the theory at the more basic mesoscopic scale. In Section 7.6 we attemptto perform such a test, especially since we can easily obtain the packings at theRLP line or the infinite compactivity limit. These packings are fully randomobtained as flat averages in the ensemble without the corresponding Boltzmannfactor and may contain direct information on the mesoscopic fluctuations.The existence of the theoretically inferred jammed states opens such predic-tions to experimental and computational investigation. We numerically test thepredictions of the phase diagram by preparing monodisperse packings of Hertz-Mindlin [35, 36] spheres with friction coefficient µ at the jamming transitionusing methods previously developed [27, 32, 65, 66].We test the theoretical predictions and show how to dynamically generate allthe packings in the phase space of configurations through different preparationprotocols. Although diverse states are predicted by the theory, they may not beeasily accessible by experimentation due to their low probability of occurrence.For instance, we will find that the packings close to the C-point (high volumefraction, high friction, low coordination) are the most difficult to obtain.The advantage of our theoretical framework is that it systematically classifiesthe different packings into a coherent picture of the phase diagram. In thefollowing we use different preparation protocols to generate all the phase spaceof jamming. In particular we provide a scheme to reproduce the RCP andRLP-lines amenable to experimental tests.We achieve different packing states by compressing a system from an initialvolume fraction φ i with a compression rate Γ in a medium of viscosity (damp-ing) η where the particles are dispersed. The system is defined by the frictioncoefficient µ which sets Z ( µ ) according to Fig. 2. While the simulations arenot realistic (no gravity, boundaries, or realistic protocol is employed), theyprovide a way to test the main predictions of the theory. The final state ( φ, Z )is achieved by the system for every ( φ i , Γ , η, µ ) at the jamming transition ofvanishing stress with a method explained next.It should be noted that other experimental protocols could be also adaptedto the exploration of the phase diagram, including (a) Gentle tapping withservo mechanisms that adjust the system at a specified pressure [66], (b) Gentle25apping with external oscillatory perturbations, (c) Settling of grains undergravity in a variety of liquids with viscosity η , (d) fluidized beds. Relation with hard sphere simulations.—
The present algorithm findsanalogies with recent attempts to describe jamming using ideas coming fromthe theory of mean-field spin glasses and optimization problems [38, 41]. This isan interesting situation since it has been shown that in the case of hard spheresthe system always crystallizes unless an infinite quench is applied [10].However, our soft sphere simulations do not produce appreciable crystalliza-tion in the packing. This situation can be understood as following: In order toavoid the (partial) crystallization in practice, a fast compression is needed. Asimple dimensionality analysis tells that the quenching rate Γ ∝ G / , where G is the shear modulus of the grains (see Fig. 7). This analysis allows for compar-ison of our soft ball simulations with the hard sphere simulations used in otherstudies [10]. One could imagine that there exists a typical quench rate, depend-ing on the stiffness of the grains, as a division of the ordered and disordered(jammed) phases. It is interesting to note that at the hard-sphere limit, G → ∞ ,all the packings are in the ordered phase except when Γ → ∞ , where φ ≈ . G . However in prac-tice, this limit is almost impossible to achieve experimentally. For sufficientlyfast quench rates, crystallization can be avoided and a granular system is stuckin a random jammed configuration since its configurational relaxation time ex-ceeds the laboratory timescale. The contention is that, for a well defined set ofquenches, the system does not crystallize and the study of the disordered phaseensues. Indeed, numerical and experimental evidence points to the validity ofthis assertion. The present theory neglects the existence of the ordered phaseand assumes an ensemble of the disordered states. A more general theory willaccount for the ordered states as well. But this is a very difficult problem injamming: a full theory that accounts for RCP and crystallization. A recentstudy has attempted to answer the question using Edwards’ thermodynamicviewpoint of granular matter [69]. From the simulation point of view, we cansay that our soft sphere approach does not produce crystallization while thehard sphere approach does not produce randomness for any finite compressionrate. We prepare static packings of spherical grains interacting via elastic forcesand Coulomb friction (see [27, 32] for more details). The system size rangesfrom N = 1 ,
024 to N = 10 ,
000 particles. Two spherical grains in contact at26ositions ~x and ~x and with radius R interact with a Hertz normal repulsiveforce [35] F n = 23 k n R / δ / , (27)and an incremental Mindlin tangential force [36]∆ F t = k t ( Rδ ) / ∆ s, (28)Here the normal overlap is δ = (1 / R − | ~x − ~x | ] >
0. The normal force actsonly in compression, F n = 0 when δ <
0. The variable s is defined such that therelative shear displacement between the two grain centers is 2 s . The prefactors k n = 4 G/ (1 − ν ) and k t = 8 G/ (2 − ν ) are defined in terms of the shear modulus G and the Poisson’s ratio ν of the material from which the grains are made. Weuse G = 29 GPa and ν = 0 . R = 5 × − m and the density of the particles, ρ = 2 × kg/m . Viscousdissipative forces are added at the global level affecting the total velocity of eachparticle through a term − γ ˙ ~x in the equation of motion, where γ is the dampingcoefficient related to the viscosity of the medium η = γ/ (6 πR ). Sliding frictionis also considered: F t ≤ µF n . (29)That is, when F t exceeds the Coulomb threshold, µF n , the grains slide and F t = µF n , where µ is the static friction coefficient between the spheres. Wemeasure the time in units of t = R p ρ/G , the compression rate in units ofΓ = 5 . t − and the viscosity in units of η = 8 . R ρ/t . We choose thetime step to be a fraction of the time unit t , which is the time that it takes forsound waves to propagate on one grain (see [66] for more details). The dynamicsfollows integration of Newton’s equation for translations and rotations. The critical volume fraction at the jamming transition, φ c will be identifiedby the “split” algorithm as explained in [32], allowing one to obtain packingsat the critical density of jamming with arbitrary precision. Preliminary resultsindicate that lower values of φ are obtained for slow compression [32]. Thus weexpect that by varying Γ or η the phase space will be explored. Also, the initial φ i plays an important role in achieving packings close to the RCP line.The preparation protocol consists of first preparing a gas of non-interactingparticles at an initial volume fraction φ i in a periodically repeated cubic box.The particles do not interact and therefore the stress in the system is σ = 0 and Z = 0.To generate a random configuration with friction at volume fraction ∼ . φ ≈ . ∼ .
36. Then we apply an extremely slow27sotropic compression without friction on this dilute configuration until the sys-tem reaches another unjammed configuration at higher density, φ i (see Fig. 8a).The reason to add an extremely slow isotropic compression to the dilute config-uration is to avoid involving any kinetic energy to keep the system as random aspossible. After obtaining this unjammed state with initial volume fraction φ i ,we restore friction and a compression is applied with a compression rate Γ untila given volume fraction φ . Then the compression is stopped and the systemis allowed to relax to mechanical equilibrium by following Newton’s equationswithout further compression.After the compression, two things can occur (see Fig. 9):(a) The system jams: If the volume fraction φ is above the jamming point, φ > φ c , then the stress will decrease and ultimately stabilize to a finite nonzerovalue, meaning that the pressure of the system remains unchanged (usually∆ σ < − Pa, red line, Fig 9) over a large period of time (usually ∼ MDsteps). The coordination number usually has a first initial decrease, but if thesystem is jammed it will also stabilize at a constant value above the isostaticminimal number (inset Fig. 9).(b) The system is not jammed: If the volume fraction φ is below the jam-ming point, φ < φ c , then the stress and the coordination number will relax tozero. This fact is illustrated in Fig. 9. If the packing has φ > φ c , it stabilizesat a non-zero pressure above the jamming transition, but the pressure decreasesvery quickly to zero (the system is not jammed) if φ < φ c , even though φ and φ differ only by 2 × − .We identify the exact volume fraction of the jamming transition φ c , as follows[32]: A search procedure consisting of several cycles is applied such that in eachcycle we fix the lower and upper boundaries of φ c , see Fig. 10. The differencebetween the boundaries gets smaller as the cycles proceed, meaning that φ c isfixed with higher and higher precision. We start from the packing of high volumefraction φ > φ c and generate a series of packings with step-decreasing volumefractions until the first packing with zero pressure is observed, which has avolume fraction φ = φ − ∆ φ , where ∆ φ is the difference between the jammmedpacking fraction φ and the unjammed packing fraction φ . Thus, φ c is boundedbetween φ − ∆ φ and φ . Then we test φ = φ − ∆ φ/
2. If φ = φ − ∆ φ/ φ c is between φ − ∆ φ and φ − ∆ φ/
2. If φ = φ − ∆ φ/ φ c is between φ − ∆ φ/ φ . Therefore, in this cycle, we reduce the regionwhere φ c possibly lies in from ∆ φ to ∆ φ/
2. If we carry out this cycle for n times, we improve the precision to ∆ φ n = ∆ φ/ n . In our simulations cyclingceases when ∆ φ n gets below 2 × − and n = 12. A similar algorithm wasemployed in [32] to study the approach to the jamming transition by preparingpackings at a finite pressure. In the present work we are interested in jammedpackings at vanishing pressure, right at the jamming transition in the isostaticplane defined for all friction coefficients.It is important to determine whether the packings are jammed in the sensethat they are not only mechanically stable but also they are stable under per-turbations. Our numerical protocols assure that the system is at least locallyjammed since each particle is in mechanical equilibrium [62]. To test if the28ystem is collectively jammed is more involved. For frictionless systems, wheretangential forces are removed, we use the Hertz energy U hertz = k n R / δ / to test whether the Hessian of the jammed configurations is positive [29]. Wefind that the frictionless configurations have positive Hessian indicating thatthey are collectively jammed. However, this method is not useful when consid-ering frictional systems. The energy of deformation depends on the path takento deform the system and cannot be defined uniquely. In this case, a numericaltest of the stability of the packings applies a small random velocity to eachjammed particle. We find that the packings are stable to small perturbationsconsisting of external forces of the order of 0.1 times the value of the averageforce, indicating that the packings may be collectively jammed.To test for more strict jammed conditions involves studying the stabilityunder boundary deformations. We have tested that our packings are stableunder the most common strain deformations to isotropic packings by performinga uniaxial compression test, a simple shear and a pure shear test.A simple shear test implies a strain deformation ∆ ǫ = ∆ ǫ = 0, while therest of the strain components ǫ ij remain constant. A pure shear test is donewith ∆ ǫ = − ∆ ǫ = 0 and a uniaxial compression test along the 1-directionis performed by keeping the strain constant in ∆ ǫ = ∆ ǫ = 0, and ∆ ǫ = 0.Here, the strain ǫ ij , is determined from the imposed dimensions of the unitcell. For example, ǫ = ∆ L/L where ∆ L is the infinitesimal change in the 11direction and L is the size of the reference state.In all cases the packings are stable under strain perturbations indicative ofthe existence of a finite shear modulus (from the shear and compression tests)and a finite bulk modulus (from the compression test) of the jammed packings.A full investigation of the elasticity of the jamming phase diagram is left forfuture studies. It suffices to state that the numerically found states can beconsidered to be mechanically stable jammed states amenable to a statisticalmechanical description in the sense of the reversible branch of compaction inthe experiments of [18, 20, 42]. We repeat the above procedure 10 times with different random initial config-urations to get a better average of φ c . Figure 11 shows the results for a system of N = 1024. All of the packings shown in the present work as well as the code togenerate them are available at http://jamlab.org. Each data point correspondsto a single set of ( φ i , Γ , η, µ ) and is averaged over these 10 realizations. Weconsider 0 . ≤ φ i ≤ .
63, 10 − ≤ Γ ≤ − , 10 − ≤ η ≤ − , and 0 ≤ µ ≤ ∞ .The error of φ c obtained over the 10 realizations as shown in Fig. 11 is usually5 × − , larger than the precision of the split algorithm (2 × − ). In general,all numerically generated jammed states lie approximately within the predictedbounds of the phase diagram.The plot in Fig. 11a explores the dependence of the packings ( φ, Z ) on theinitial state φ i . (In the plot, φ refers to the φ c obtained in a split algorithm).In Fig. 11 we plot our results for a fixed quench rate Γ = 10 − and dampingcoefficient η = 10 − (except for the last orange curve on the right with η =290 − ) and for different initial states ranging from left to right (see Fig. 11for details) Each color here corresponds to a particular φ i (0.40 at the left to0.63 at the right), while each data point along a curve corresponds to a systemprepared at a different friction from µ = 0 at Z ≈ µ → ∞ at Z ≈ µ . While these lines are approximately horizontal, indicatingan approximate independence of Z on the preparation protocol φ i , we observedeviations specially around the RLP line for intermediate values of µ . Despitethis, the Z ( µ ) dependence is approximately valid as seen in Fig. 2, which isplotted using the results of Fig. 11a.We find that the packings prepared from the larger initial densities φ i closelyreproduce the RCP line of zero compactivity at φ RCP , confirming the predictionthat RCP extend along the vertical line from the J-point to the C-point.We find that the packings along the RCP line have equal geometrical coor-dination number z ≈ Z = 6 to Z ≈ µ as depicted in the volumelandscape picture of Fig. 5.In the limit of small initial densities (see curves for φ i = 0 . , .
53 and0 .
55 in Fig. 11a) we reproduce approximately the predictions of the RLP-line.These packings follow the theoretical prediction for infinite compactivity exceptfor deviations (less than 5%) in the coordination number for packings closeto the lower value of Z = 4, where the numerical curve seems to flatten outdeviating from the theory. Numerically, we find the lowest volume fraction at φ minRLP = 0 . ± . φ minRLP = 0 . h z = e − , thusthe minimum RLP deviates slightly from the theoretical prediction in the limit h z → Z correspond to infinitely rough spheres. Thedeviation of coordination number between theory and simulation (specially atlow volume fraction) could be from the system not achieving an isostatic stateat infinite friction. We were not able to generate packings with exactly Z = 4but our algorithm generates packings very close to this isostatic packings at Z ≈ . µ = 10 . In general, we find that while there are many states alongthe RLP line, the barriers of these states decrease as the coordination numberdecreases towards Z = 4, i.e. when the friction increases. Thus, the states atthe lower left part of the phase diagram are very difficult to equilibrate.Besides the states delimiting the phase space, we generate other packingswith intermediate values of φ i = 0 . , . , .
61 as shown in Fig. 11a. Interest-ingly, we find that these states (all obtained for fixed Γ = 10 − and η = 10 − )closely follow the predicted lines of isocompactivity as indicated in the figure.We find that the simulations from φ i = 0 . , . , .
61 correspond to compactiv-ities X = 1 . , . , and 1.16, respectively (measured in units of 10 − V g ).The constant h z weakly affects the finite compactivity lines. We find that h z = e − provides the best fit to the data for finite isocompactivity lines inFig. 11a. Thus, with reasonable approximation and for this particular protocol,30e identify the density of the initial state, φ i with the compactivity of thepacking, providing a way to prepare a packing with a desired compactivity. Wenote, however, that other values of h z produce approximately the same phasediagram boundaries (i.e., the RLP and RCP lines, as long as h z ≪
1) but donot fit the isocompactivity lines within as accurately. Thus, the identificationof the compactivity with φ i remains dependent on this particular value of h z used to fit the theory.We also test the dependence on the state of the packings with the compres-sion rate Γ and viscosity η . Both parameters should have similar effects sincethey slow down the dynamics of the grains. Figures 8b and 8c represent typ-ical preparations. In general, we reproduce the RLP line for slow quenches orfor large viscosities as seen in Fig. 11b (for φ i = 0 .
40 in black). In this case,the grains are allowed to slide and develop large transverse displacements andMindlin forces with the concomitant low φ and large compactivity. Thereforewe find a predominance of the Mindlin forces over the normal Herztian forcesas a characteristic of the lower volume fractions of RLP, having a pronouncedpath dependent structure with large compactivity in the packing. When thecompression rate is increased or the damping is reduced we find packings withhigher volume fractions as indicated in Fig. 11b. Indeed, analogous proto-col to generate these packings would be to change the damping coefficient byimmersing the particles in liquids of different viscosities.It is also interesting to investigate the packings generated by fixed frictionand varying quenched rate, since this situation corresponds to a given system.For a given µ we find approximately a common value of Z ( µ ) independent ofΓ and φ i , as discussed in Fig. 2, indicating that the curve Z ( µ ) is roughlyindependent of the protocol. Thus, a fixed µ corresponds to a horizontal line inthe phase diagram as indicated in Fig 11a by the dashed lines.Figure 8 shows a summary of the numerical protocols we use to generatethe packings in the phase diagram. The “easiest” packings to generate arealong the RLP line (Fig. 8b, c) while the most difficult are deep in the phasediagram with low coordination number (large friction) and high volume fraction,that is near the C-point (Fig. 8a). In general, we always start with a diluteunjammed (non-interacting gas) sample without friction and compressed slowlyuntil another dilute sample not jammed at a volume fraction φ i . If φ i is aroundor smaller than 0 .
55 then the RLP line is obtained. If we decrease by 3 ordersof magnitudes the compression rate using this initial density, we move insidethe phase diagram, though not substantially, as indicated schematically in thefigure.Therefore, in order to access the inner states around the C-point at highfriction, we took another route. We prepared a frictionless dilute unjammedsample at a higher φ i , which can be as close as possible to ∼ .
63 since friction-less spheres only jam at RCP. Then we switched on friction and got the packingsinside the phase diagram as indicated in Fig. 8a. We realize that switching onand off friction is impossible in experiments. However, a possible way to exper-imentally realize the path to the C-point is by resetting the path dependencyto zero along the preparation protocol: this path could be reproduced by stabi-31izing a frictional packing with fast compressions and resetting the interparticleforces by gently tapping. This will simulated the lost of path dependence thatwe simulate by the initial compression without friction and may allow to reachthe packings around the C-point.
At this point we do not have a theoretical explanation for why packings withthe same initial unjammed state φ i end up having the same compactivity whenjammed after compression. It is important to note that this finding is based onthe particular value we use for h z , such that a change in h z gives rise to differentlines of isocompactivity which may not fit the simulation results of Fig. 11a soclosely. Nevertheless, the close fit between theory and simulation deserves anexplanation.We may conjecture that φ i determines a type of disorder quenched in theinitial configuration that leads to systems with the same compactivity but dif-ferent volume fractions and coordination numbers, evidenced by our results.We can use this empirical result to control the compactivity of the packing, atleast for this particular protocol, defining a “granular compactimeter”. Labo-ratory measurements of compactivity usually involve indirect measures throughFluctuation-Dissipation relations between fluctuations and response functions[18, 19, 24]. However, there is no simple thermometer (or ”compactimeter”) tomeasure this observable directly in a packing. Furthermore, there is no sim-ple way to prepare a packing with a desire compactivity. Our results can beused to, at a minimum, define packings with equal compactivity. The empiricalidentification of X with φ i promotes the possibility of controlling X within thisparticular protocol, an important step in any thermodynamical analysis.The advance of a granular thermodynamics crucially depends on the inven-tion of a thermometer that can easily measure the compactivity, as in Fig. 12.A zero-th law of thermodynamics presents serious challenges in granular ma-terials since there is no straight forward way to define a compactivity bath orreservoir.The present theory predicts the dependence of the volume fraction, φ ( X ),on the compactivity, allowing for the development of a compactimeter: Such adevice consists of a known granular medium, having a given µ and, therefore, agiven coordination number. The compactimeter is inserted in a granular systemand both, system and compactimeter, are gently shaken to allow for equilibra-tion at a common compactivity (the compactimeter has a flexible membranethat allows volume transfer between both systems in contact). The volume ofthe grain in the compactimeter will adjust accordingly and the compactivitycan be read directly from the scale attached to it through the equation of state, X ( φ ), provided by the theory.The foundation of the above compactimeter relies on the validity of a zero-thlaw for granular thermodynamics that determines the equilibration at a givencompactivity of a system: Two systems in contact should mechanically equili-brate at the same compactivity. A recent experiment [70] found a materials-independent relationship between the average volume fraction and its fluctua-32ions in two equilibrated granular subsystems, which gives support to the zero-thlaw. Further test of the zero-th law is difficult to facilitate without theoreticalguidance. Thus, it will be useful to perform the following test, either numer-ically or experimentally to test not only the idea of equilibration but also thepossibility to describe granular matter under the V-ensemble. Such a test hasbeen labeled the ABC experiment as proposed by Edwards [45], which is nowpossible to perform using the phase diagram.Two packings are numerically prepared at points A and B in the Fig. 13 fordifferent µ , µ A < µ B . The systems are equilibrated at different compactivities X A < X B by using different φ i . They are then put into contact through aflexible membrane (or by simply putting the particles at the surface in contact)and allowed to mechanically equilibrate by gently shaking the system. Numer-ically, this can be simulated with the same compression methods as explainedin Section 7. The total volume is kept fixed and the volumes of the subsystemsshould change accordingly, V = V A + V B . If the AB system equilibrates at thesame compactivity, then it will follow the trajectory depicted in Fig. 13 towardsthe state C of equilibration along the isocompactivity line. There is no needto measure the compactivity of the final packing. By just measuring the finalvolumes V A and V B (or their volume fractions) will suffice. This will provide animportant test for the V-ensemble and compactivity to describe jammed matter,as well as the validity of the phase diagram.Another intriguing possibility would be to mix a frictionless A packing at theJ-point with an infinite friction B packing at the L-point, which can be achievednumerically. This AB system would provide the maximum difference betweenthe isolated packings, allowing study of the zero-th law with more accuracy.According to theory, the frictionless packing is independent of X . Therefore,it should remain at the J-point. The µ → ∞ packing should, in principle,equilibrate at any X or volume fraction along the G-line. Most probably, itwould stay put at the L-point since the A and B packings are joined by theisocompactivity X → ∞ line already. In any case, the resulting AB packingshould fall inside the phase diagram if the theory is correct. If it does notequilibrate (falling outside the phase diagram, we note that the two packingsare at the borderline), it could indicate the breakdown of the approach. Twoconclusions could be reached from the failure of such a test. Either the V-ensemble is wrong or the approximations of the V-ensemble theory are notcorrect. In the latter, more sophisticated theories should be developed. If theformer case occurs, then a different ensemble will have to be considered beyondisostaticity. The notion that RCP arises for X = 0, implying that there are no fluctua-tions at RCP, deserves some discussion. Furthermore, it is important to discussthe implications of the related coincidence of RCP and RLP at the frictionlesspoint. That is, the fact that for frictionless particles there is a unique state ofjamming and therefore the compactivity is irrelevant when µ = 0.33irst of all, we need to emphasize that these results should be interpretedat the mesoscopic level. Thus, there is one mesoscopic state at RCP and onlythe mesoscopic entropy vanishes at this point. These results provides an un-ambiguous interpretation of RCP only at the meso level. As discussed above,for a single mesoscopic state, there are many microscopic states average out inthe calculation of the average free volume function w ( z ) which is used in thepartition function Eq. (19) (see also Jamming I [49]). Therefore, we expectthat these microscopic states will contribute to the entropy of RCP (see Jam-ming III [51] for a discussion) as well as giving rise to other states perhaps withdifferent volume fraction. Below we elaborate on the existence of microscopicand mesoscopic packing states. This point is better discussed in terms of thevolume landscape of Fig. 5: all the jammed states are degenerated around themesoscopic ground state with the coordination number z = 6. These states haveslightly different volume fractions, leading to microscopic fluctuations which arecoarse-grained in the mesoscopic theory.The existence of microscopic states neglected at the mesoscopic level raisesan interesting analogy with recent work of Zamponi and Parisi [38] where thejammed states are considered as infinite pressure glassy states obtained after afast compression from a liquid state at finite temperature. They propose thata range of different volume fractions of jamming can be achieved according tothe initial state from where the quench is started, with the maximum possiblejammed density obtained when the initial density is the Kauzmann point ofthe ideal glass. In their representation (see Fig. 4 in [38]) a range of initialdensities in the liquid phase ϕ ∈ [ ϕ d , ϕ K ], where ϕ d is the density where manymetastable states first appear in the liquid phase as suggested by mean fieldmodel picture of the glass transition, and ϕ K is the Kauzmann density of theideal glass transition, gives rise to a final density of jamming, ϕ j ( ϕ ), obtainedafter a fast quench from ϕ . Two limiting densities of jamming are predicted:the minimum at ϕ j ( ϕ d ) = ϕ th and the maximum volume fraction obtained aftera quench from the Kauzmann point, ϕ j ( ϕ K ) = ϕ .This situation could be interpreted in term of neglecting the microscopicstates in our theory. It is also a question whether the states at infinite pressureobtained using fast quenches from a liquid state “a la Zamponi” are equivalentto the jammed state obtained in our approach for soft spheres. Our volumeensemble approach might differ from the energy/force ensemble used in the workof Zamponi-Parisi. Further investigations are required to reveal the possibledifferences. In Section 7.3 we have tested the predictions of the theory. In this Sectionwe attempt to test some of the assumptions of the theory, specially those re-lated to the quasiparticles of w ( z ) and the bounds in z and its distribution. Aquasiparticle is a mathematical entity that is, in principle, impossible to isolatefrom the structure of a real packing prepared either with MD computer simu-lations or experiments, rendering the measurement of their properties from real34ackings difficult. This is because real packings already contain the ensembleaverage, and we can not dissociate the ensemble, with the Boltzmann distribu-tion of states, compactivity and density of states, from the isolated statistics ofa quasiparticle.Quasiparticles can be tested by preparing random packings with prescribedgeometrical coordination number. However, below, we argue that we may obtainapproximate properties from real packings in the limit X → ∞ and h z → X → ∞ the Boltzmann factoris one, and the average in the partition function is flat over the configurationsup to the density of states. Furthermore, the density of states is found tobe a very fast decaying function of z . These two results implies that along theRLP-line of real packings (either numerical or experimental) we can obtain fullyrandom packings which could reveal the properties of the quasiparticles. Underthis approximation we can then test some of the approximations of the theory.However, the conclusions from this section remain approximate.We have proposed that the bounds of the geometrical coordination numberare Z ≤ z ≤ z is averagedover this range, we find that z ranges approximately between ( Z, z without the microscopic fluctuations but includingthe mesoscopic fluctuations. Once we average over its neighbors, we not onlyremove the microscopic fluctuations but also, partially, the mesoscopic ones. Inthe limiting case of no mesoscopic fluctuations, when X = 0 along the RCPline, we find a very narrow distribution at z ≈ Z and µ , afteraveraging z over a mesoscopic region of two particle diameters. The distributionis even narrower when z is coarse-grained over a region of four particle diameters.The X-ray tomography experiments of Aste [46] (Fig. 6a) reveal the trendpredicted by Eq. (7) between the inverse of the Voronoi volume and the numberof neighbors of a set of Voronoi cells with similar volumes. However, it is evidentfrom the figure that such number of neighbors is spanning a range of valuesbetween ∼ ∼
10. The reason why the Aste’s group [46] found a widerrange of z is possibly due to the fact they did not consider a range of coarse-graining.A question arise as how to identify the geometrical from the mechanical co-ordination number, even numerically. A detailed discussion of this topic appearsin Ref. [71]. Below we summarize the results. Strictly speaking, this distinc-tion is not possible to materialize in real packings. Again, we need to generaterandom packings from a geometrical point of view, without resorting to forces,and then study their distribution in real packings where forces are taken intoaccount. Nevertheless, we follow the following approximate scheme to try toobtain information of the quasiparticles from real packings.Even in numerical simulations, there is always a round off error associatedto measurements. Thus, two particles that may not be in contact (giving riseto a zero force) may be close enough to be considered as contributing to thegeometrical coordination. Indeed, it is known that g ( r ) has a singularity [72],35 ( r ) ∼ ( r − . − . , implying that there are many particles almost touching. Fol-lowing this consideration, we introduce a modified radial distribution function(RDF) g z ( r ) in order to approximately identify z and Z from real packings: g z ( r ) = 1 N R r N X i N X j = i Θ (cid:16) r ij r − R − (cid:17) Θ (cid:16) r + Rr ij − (cid:17) , r > R (30)where R is the radius of particle, N is the number of particles, r ij is the distanceof two particle’s centers, r ij = | ~r i − ~r j | , and Θ is the Heaviside step function.The RDF describes the average value of the number of grains in contact with avirtual particle of radius r , and the factor of R /r is the ratio of a real sphere’sarea and the virtual one’s. g z ( r ) measures the number of balls with their volume intersecting the surfaceof a sphere of radius r measured from the center of a given ball. When r = R in (30) we obtain the mechanical coordination number while the geometricalone is obtained for a small value ∆ r = r − R R = 0 for which we distinctly find asignature from computer simulations, unambiguously defining it at ∆ r = 0 . g z (∆ r ) of packings with various friction coefficient µ on the isostatic plane along the RCP and RLP lines respectively. Followingthe definition of Eq. (30), g z , with ∆ r = 0, should be directly equal to themechanical coordination number, Z , and should range from 4 to 6 along bothRCP and RLP (if h z ≪
1) lines which is confirmed by our numerical simulationsin Figs. 14 and 15, respectively.Furthermore, as shown in the figures, we find that g z (∆ r ) along the RCP line,increases slightly as ∆ r increases, and finally reach the same value of g z at ∆ r =0 .
04 as shown in Fig. 14. We identify the geometrical coordination number as z = g z (0 .
04) under the accuracy of the simulations. Therefore, we conclude thatall RCP states have approximately the same geometrical coordination number, z ≈
6, in agreement with the theory. In terms of the volume landscape of Fig. 5,the states along the RCP line are ground states with different friction at X = 0.Thus they should all have z = 6. We notice that g z (∆ r = 0 .
04) = 6 . > z from real packings since we are not measuring the quasiparticles directly, and,secondly, from the increasing of the coordination number slightly away fromthe critical point as the local volume fraction φ increases slightly as Z − Z c ∼ ( φ − φ c ) β . We note that Z is also slightly increased from the isostatic point.Along the RLP line, Fig. 15, we find that the geometrical coordinationnumber as extracted from g z (0 .
04) is very close to the mechanical one. Sincethe RLP line is at X → ∞ and h z ≪
1, all the states along RLP have z ≈ Z as we move along the line varying the friction coefficient. Thus, the numericalresults confirm the theory, and further confirm that the value of h z is very small.Next, we study the coarse-grained coordination number h z i l as defined as, h z i l = 1 N l R r N l X i N l X j = i Θ (cid:16) r ij r − R − (cid:17) Θ (cid:16) r + Rr ij − (cid:17) , r > R (31)36here N l is the number of the particles inside a coarse-grained spherical rangewith a radius of l . Figures 16 and 17 plot the PDF of h z i l for all the packingsalong the RCP and RLP lines, respectively. The distributions show a narrowshape, and the average value gives the value of g z ( r ), i.e., h z i l = g z ( r ) . (32)We find that all the P ( h z i l ) along the RCP line coincide at ∆ r = 0 .
04 as shownin Fig. 16b and 16c, demonstrating that all the states have approximately thesame average value of z (as discussed above), as well as the same distribution.The distribution gets narrower as the coarse-graining parameter l increase from l = 2 to l = 4 as shown in Fig. 16c. On the other hand Z changes as the frictionchanges following the isostatic condition, as seen in Fig. 16a. This in in goodagreement with the theory.Along the RLP line, the situation is analogous. As discussed above, z ≈ Z .This is clearly demonstrated in Fig. 17b which should be compared with theanalogous Fig. 16b for the RCP line. Figure 17a also shows how Z changeswith friction in the same way as in Fig. 16a. More importantly, we find thatthe bounds of h z i l are well approximated between the bounds proposed by thetheory ( Z, z ≈
6, with Z changing from 6 to 4 as friction is increased. On the otherhand, along the RLP line, we find that z ≈ Z . These two results confirm thetheory very well. Figure 18 summarizes all these finding indicating the rangesof z and Z along the iso- X lines of RCP and RLP, the iso- z vertical lines, andthe iso- Z horizontal lines in the phase diagram.
8. Outlook
We have presented several results in the disordered sphere packing problemand here will discuss their significance. We start with the results that we believeare more robust and may survive higher scrutiny, and follow with the resultsthat require further study plus a discussion on how to improve the theory.The final result of the phase diagram, the characterization of RCP and RLPusing the compactivity, the role of friction in the determination of the limitingpackings are interesting results, almost independent of the theory. The dis-tinction between the geometrical and mechanical coordination number is alsoimportant, allowing for the implementation of averages in the partition function.On the mathematical side, the microscopic volume function and the formula ob-tained for the Voronoi cell in Jamming I [49] are important and exact resultswhich are the foundation of all preceding work presented herein.While the predictions Eqs. (24) and (25) of RCP and RLP approximate verywell the known values, they are obviously not exact expressions, since they arebased on several approximations of the statistical theory of volumes. Thus, the37mplementation of the mesoscopic volume function and the probability distribu-tion of volumes on which these results are based are, perhaps, the points whereimprovements may be necessary. While the √ z formula reproduces the averagevolume fraction surprisingly well, it is based on approximations to calculate theprobability distribution of volumes, P > ( c ), as explained in Jamming I [49]. Morework is needed to obtain better approximations to such a distribution to cap-ture not only the mean value of the Voronoi volume but higher moments as well.A full discussion is performed in Jamming IV [52]. The present theory showsthe way to improve upon Jamming I and obtain an exact formulation of P > ( c ),at least to a prescribed value of the coordination shell. Such a formulation ispresently being developed and is based on a systematic study of quasiparticleswith fixed geometrical coordination number; it produces good fits of P ( W si ) forthe case of 2d, while 3d calculations are under way.Other approximations that require additional attention are the considera-tions leading to the density of states g ( z ). While the main results of RCP andRLP are nearly independent of the particular form of g ( z ), other quantitiessuch as the entropy and the equations of state for finite nonzero compactivityare more sensitive to it. By a direct measurement of the entropy, g ( z ) can becalculated and compared to the simple exponential form used in the presentstudy.Furthermore, we have explicitly removed crystal states from our calculationsby considering the upper bound of z = 6. A more general theory is needed thento investigate the transition from the RCP to the FCC as crystals starts to formin the system. Such a theory will include not only disordered states but alsopartially crystallized states and would give X = 0 at FCC following the fullEdwards partition function, Eq. (5). We investigate this idea in [69], where wefind that RCP can be interpreted as the “freezing point” in a first-order phasetransition between ordered and disordered packing phases.The distribution of coordination number, assumed to be a delta function, isanother approximation to be improved upon. Once again, this may not affectthe predicted RCP and RLP but deeper studies require also a distribution of z .A more general ensemble is being considered that predicts the distribution of z in good agreement with the existing data (see Jamming IV [52]).Finally, the extension to study the microstates requires more investigation.The identification of the elementary units as quasiparticles of fixed coordinationnumber is interesting and could eventually lead to more precise solutions of RCPand RLP.
9. Conclusions
In conclusion, using Edwards statistical mechanics we have elucidated someaspects of RLP and RCP in the disordered spherical packing problem. Thenumerical results suggest a way to experimentally test the existence of the pre-dicted packings. By allowing the grains to settle in liquids of varying density,the speed of the particles can be varied and a systematic exploration of thejamming phase diagram can be achieved.38eyond the elucidation of some questions in the sphere packing problem,other problems can now be addressed systematically from the point of view ofwhat we have learned from the phase diagram. This includes the investigation ofthe criticality of the jamming transition from frictionless to frictional systems byextending the phase space to ( φ, Z, σ ), and generating packings away from theisostatic plane by isotropically compressing the packings generated in the planeof Fig. 4. Previous studies have focused on the frictionless point and on fixedfriction. According to our results it is important to sample all the phase spacefor different friction and compactivities according to the preparation protocolsexplained above.Previous results find power law scaling [29, 30, 27, 32] where the stressvanishes as σ ∼ ( φ − φ c ) α . A possible scenario is depicted in the extendedphase diagram of Fig. 19 in the 3-dimensional space ( φ, Z, σ ), which could betested numerically and experimentally. We note in passing that, all the resultsin the present study refer to the isostatic plane at the jamming transition σ = 0or hard sphere limit. We expect that, as we compress packings of soft spheres,different planes of σ = constant = 0 will be obtained as shown in Fig. 19. Thesecompressed states are described not only by the compactivity X but also theangoricity A through the full partition function in the VF ensemble, Eq. (2).Our preliminary results indicate that there seems to be evidence of criticalityas the isostatic plane is approached when σ →
0, although with exponentsdependent on friction as well as the value of the compactivity.Additional topics to be addressed by means of the presented approach in-clude characterization of jamming in the phase space of configurations, theproblem of elasticity and Green’s function and the study of the pair correla-tion function g ( r ), the volume distribution P ( W si ), Voronoi volume, P ( W vor i ),and coordination number P ( z ) , which are within the reach of the present ap-proach. Other systems under consideration are two-dimensional systems [53],polydisperse systems [54], as well as the mean field case of infinite dimension.These are interesting cases, since there are many questions debated as to whatis the more dense state in higher dimension. Extensions to other systems suchas anisotropic particles are also within the reach of the present approach.Finally, a complete characterization of the VF-ensemble can now be per-formed through the measurement of other quantities as a function of the volumefraction, such as the force distribution P ( F n , F t ), and equations of state throughthe angoricity in the ( φ, Z, σ ) space. Some of these quantities can be obtainedanalytically from the ensembles allowing for interesting predictions.The phase diagram introduced here serves as a beginning to understand howrandom packings fill space in 3d. The comparative advantage of the presentapproach over extensive work done in the past, is in the classification of allpackings through X , Z and φ in the theoretical phase diagram from wherethese studies could be systematically performed. This classification guides thesearch for indications of jamming from a systematic point of view, through theexploration of all jammed states from µ = 0 to µ → ∞ . Acknowledgements . This work is supported by the National Science Foun-dation, CMMT Division and the Department of Energy, Office of Basic Energy39ciences, Geosciences Division. We are grateful to B. Bruji´c and A. Yupan-qui for inspirations, C. Briscoe for a critical reading of the manuscript and thehospitality of UFC, Fortaleza where part of this work was done.
References [1] T. C. Hales, http://arxiv.org/-abs/math.MG/9811078.[2] A. Thue, Forandlingerneved de Skandinaviske Naturforskeres, , 352-353(1892).[3] A. Thue, Christinia Vid. Selsk. Skr. , 1-9 (1910).[4] J. D. Bernal, and J. Mason, Nature , 910 (1960); J. D. Bernal, Nature,185 (1960); J. D. Bernal, Proc. Roy. Soc. London, Ser. A , 299 (1964).[5] Anonymous, Nature , 488 (1972).[6] G. D. Scott, and D. M. Kilgour, Brit. J. Appl Phys (J. Phys. D) , 863-866(1969); G. D. Scott, Nature, , 908 (1960).[7] J. L. Finney, Proc. Roy. Soc. London, Ser. A , 479 (1970).[8] J. D. Berryman, Phys. Rev. A , 1053-1061 (1983).[9] G. Y. Onoda, and E. G. Liniger, Phys. Rev. Lett. , 2727 (1990).[10] S. Torquato, T. M. Truskett, P. G. Debenedetti, Phys. Rev. Lett. , 2064(2000).[11] R. P. Behringer and J. T. Jenkins, eds., Powders & Grains 97 (Balkema,Rotterdam, 1997).[12] A. Coniglio A, A. Fierro, H. Hermann and M. Nicodemi, eds, Unifyingconcepts in granular media and glasses, (Elsevier, 2004).[13] S. F. Edwards and R.B.S. Oakeshott, Physica A , 1080-1090 (1989).[14] S. F. Edwards, The role of entropy in the specification of a powder, inGranular matter: an interdisciplinary approach (ed A. Mehta) 121-140,Springer-Verlag, New York, 1994.[15] R. Blumenfeld, and S. F. Edwards, Phys. Rev. Lett. , 114303 (2002).[16] S. F. Edwards, Physica A , 114 (2005).[17] M. P. Ciamarra, A. Coniglio, M. Nicodemi, Phys. Rev. Lett. , 158001(2006).[18] E. R. Nowak, J. B. Knight, E. BenNaim, H. M. Jaeger and S. R. Nagel,Phys. Rev. E , 1971 (1998). 4019] M. Schr¨oter, D. I. Goldman, H. L. Swinney, Phys. Rev. E , 030301(R)(2005).[20] J. Bruji´c, P. Wang, D. L. Johnson, O. Sindt, and H. A. Makse, Phys. Rev.Lett. , 128001 (2005).[21] E. Bertin, O. Dauchot and M. Droz, Phys. Rev. Lett. , 120601 (2006).[22] L. F. Cugliandolo, J. Kurchan, L. Peliti, Phys. Rev. E , 3898 (1997).[23] A. Barrat, J. Kurchan, V. Loreto, M. Sellitto, Phys. Rev. Lett. , 5034(2000).[24] H. A. Makse and J. Kurchan, Nature , 614 (2002).[25] A. Fierro, M. Nicodemi, M. Tarzia, A. de Candia, A. Coniglio, Phys. Rev.E , 061305 (2005).[26] A. Liu and S. R. Nagel, Nature 396, 21-22 (1998).[27] H. A. Makse, D. L. Johnson, L. M. Schwartz, Phys. Rev. Lett. , 4160(2000).[28] V. Trappe, V. Prasad, L. Cipelletti, P. N. Segre, D. A. Weitz, Nature ,(2001).[29] C. S. O’Hern, S. A. Langer, A. J. Liu, and S. R. Nagel, Phys Rev. Lett. , 075507 (2002).[30] L. E. Silbert, D. Ertas, G. S. Grest, T. C. Halsey, D. Levine, Phys. Rev. E , 031304 (2002).[31] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E ,011306 (2003).[32] H. Zhang and H. A. Makse, Phys. Rev. E , 178001 (2005).[34] K. Shundyak, M. van Hecke and W. van Saarloos, Phys. Rev. E , 010301(2007).[35] L. D. Landau and E. M. Lifshitz, Theory of Elasticity, (Pergamon, NY,1970).[36] R. D. Mindlin, J. Appl. Mech. (ASME) , (1949).[37] J. Bruji´c, S. F. Edwards, D. V. Grinev, I. Hopkinson, D. Bruji´c, and H. A.Makse, Faraday Discuss. , 207 (2003).[38] G. Parisi, and F. Zamponi, Rev. Mod. Phys. , 789-845 (2010).[39] B. D. Lubachevsky and F. H. Stillinger, J. Stat. Phys. , 561 (1990).4140] M. Skoge, A. Donev, F. H. Stillinger, and S. Torquato, Phys. Rev. E ,041127 (2006).[41] F. Krzakala and J. Kurchan, Phys. Rev. E , 021122 (2007).[42] P. Philippe, and D. Bideau, Europhys. Lett. , 677 (2002).[43] S. Henkes and B. Chakraborty, Phys. Rev. Lett. , 198002 (2005).[44] R. Blumenfeld, On entropic characterization of granular materials, in Lec-ture Notes in Complex Systems Vol. 8: Granular and Complex Materials,43-53, eds (2007).[45] H. A. Makse, J. Bruji´c, and S. F. Edwards, Statistical Mechanics of JammedMatter, in The Physics of Granular Media, edited by H. Hinrichsen and D.E. Wolf (Wiley-VCH, 2004).[46] T. Aste, M. Saadatfar, T. J. Senden, J. Stat. Mech., P07010 (2006).[47] F. Lechenault, F. da Cruz, O. Dauchot, and E. Bertin, J. Stat. Mech.,P07009 (2006).[48] L. D. Landau and E. M. Lifshitz, Statistical Physics (Pergamon, NY, 1970).[49] C. Song, P. Wang, Y. Jin, and H. A. Makse, Physica A (2010),doi:10.1016/j.physa.2010.06.043.[50] C. Song, P. Wang, and H. A. Makse, Nature , 629 (2008).[51] C. Briscoe, C. Song, P. Wang, and H. A. Makse, Physica A , 3978(2010).[52] P. Wang, C. Song, Y. Jin, and H. A. Makse, e-print arXiv:0808.2198.[53] S. Meyer, C. Song, Y. Jin, K. Wang, H.A. Makse, Physica A (2010),doi:10.1016/j.physa.2010.07.030.[54] M. Danisch, Y. Jin, and H. A. Makse, Phys. Rev. E. , 051303 (2010).[55] J. A. van Meel, B. Charbonneau, A. Fortini, and P. Charbonneau, Phy.Rev. E. , 061110 (2009).[56] R. C. Ball and R. Blumenfeld, Phys. Rev. Lett. , 115505-1 (2002).[57] R. Blumenfeld, Eur. Phys. J. B , 261 (2005).[58] J. Zhou, S. Long, Q. Wang, Science , 1631 (2006).[59] S. Alexander, Phys. Rep. , 65 (1998).[60] C. F. Moukarzel, Phys. Rev. Lett. , 1634 (1998).[61] S. F. Edwards and D. V. Grinev, Phys. Rev. Lett. , 5397 (1999).4262] S. Torquato, F. H. Stillinger, J. Phys. Chem B , 11849 (2001).[63] J. Bruji´c, C. Song, P. Wang, C. Briscoe, G. Marty, and H. A. Makse, Phys.Rev. Lett. , 248001 (2007).[64] F. H Stillinger, Science , 1935-1939 (1995).[65] H. A. Makse, N. Gland, D. L. Johnson, and L. M. Schwartz, Phys. Rev.Lett. , 5070 (1999);[66] H. A. Makse, N. Gland, D. L. Johnson, and L. M. Schwartz, Phys. Rev. E , 061302 (2004).[67] F. Zamponi, Nature , 606 (2008).[68] S. Mukhopadhyay, M. L. Gilchrist, and H. A. Makse (in preparation).[69] Y. Jin and H. A. Makse, Physica A (2010), doi:10.1016/j.physa.2010.08.010.[70] F. Lechenault and K. E. Daniels, e-print arXiv:1001.5411.[71] C. Song, P. Wang, and H. A. Makse, AIP Conf. Proc. 1227, 271 (2010).[72] A. Donev et al. , Phys. Rev. E , 011105 (2005).43 ij (ˆ s ) l i (ˆ s ) i j r ij ˆ s (a) θ ij Figure 1: The Voronoi volume is the light grey area (shown in 2d for simplicity). The limitof the Voronoi cell of particle i in the direction ˆ s is l i (ˆ s ) = r ij / θ ij . Then the Voronoivolume is proportional to the integration of l i (ˆ s ) over ˆ s as in Eq. (6). -7 -5 -3 -1 Z i = 0.40 i = 0.53 i = 0.55 i = 0.57 i = 0.59 i = 0.61 i = 0.63 Figure 2: Mechanical coordination number versus friction µ obtained in our numerical simu-lations explained in Section 7 for different preparation protocols characterized by the initialvolume fractions φ i indicated in the figure. The symbols and parameters used in these simu-lations are the same as in the plot of Fig. 11. a) Isostatic packing z = 6 µ = Z = 6 (b) Z = z = 6 µ = ∞ Figure 3: (a) Consider a frictionless packing at the isostatic limit with z = 6. In this case theisostatic condition implies also Z = 6 mechanical forces from the surrounding particles. (b)If we now switch on the tangential forces using the same packing as in (a) by setting µ → ∞ ,the particle requires only Z = 4 contacts to be rigid. Such a solution is guaranteed by theisostatic condition for µ → ∞ . Thus, the particle still have z = 6 geometrical neighbors butonly Z = 4 mechanical ones. Z RLP lineX RCP lineX=0 Figure 4:
Phase diagram of jammed matter: Theory . Theoretical prediction of thestatistical theory. All disordered packings lie within the yellow triangle demarcated by theRCP line, RLP line and G line. Lines of finite isocompactivity are in color. The grey area isthe forbidden zone where no jammed packings can exist. O disordered RCP
RLP z = 5 z = 6 z = 4 FCC W O RCP
RLP z = 5 z = 6 FCC W O RCP=RLP z = 6 FCC
Energy O microscopic states µ = ∞ finite µµ = 0 ! r i ! r i ! r i ! r i Isostatic plane packings Z = Z ( µ ) = 5 Z = (path dependent) (a) (b)(c)(b) (d) W mesoscopic state Figure 5:
Schematic representation of the volume landscape of jammed matter ( ~r i , W ) . The multidimensional coordinate ~r i represents the degrees of freedom: the particlepositions. Each dot represents a discrete mesoscopic jammed states at different z . It isimportant to note that for each meso state there are microscopic states with the same z .All disordered packings are in the yellow region (of an schematic shape) of the phase spacewhich corresponds to the isostatic plane of hard spheres at the jamming transition where ourcalculations are performed. Other ordered packings have lower volume, such as the FCC. (a)We represent the case of µ = ∞ . The states represent those along the G-line in Fig. 4 as thecompactivity varies from X = 0 (ground state) to X → ∞ at the RLP. The horizontal linesindicate packings at constant volume. The ground state of jammed matter for this frictioncoefficient has z = 6 and the highest volume states are found for z = 4. The arrow indicatesthe limits of integration in the partition function for this particular friction. (b) For anotherfinite µ , the space is delimited by above by a line of constant z = Z ( µ ). (c) For µ = 0 only theground state is available, giving rise to a single state. (d) The volume states states in (a,b,c)are separated by energy barriers represented by the third coordinate in the phase space. Theenergy barriers are path dependent due to friction between the particles. Nevertheless, thejammed configurations are well-defined in the isostatic plane, and the energy barriers representthe work done to go from one configuration to another. Only at the frictionless point, theenergy barriers are path independent. -4 -3 -2 -1 minRLP RCP Z =4.0 Z =4.5 Z =5.0 Z =5.5 S X / V g Figure 6:
Predictions of the equation of state of jammed matter in the ( X, φ, s ) space. Each line corresponds to a different system with Z ( µ ) as indicated. The projectionsin the ( φ, s ) and ( X, s ) planes show that the RCP ( X = 0) is less disordered than the RLP( X → ∞ ). The projection in the ( X, φ ) plane resembles qualitatively the compaction curvesof the experiments [18, 19, 20]. G Γ orderdisorder Figure 7: Dimensional analysis indicates that for finite shear modulus G of the particles, thesystem can crystallize or randomize according to the quenched rate Γ. This argument allowsfor a comparison with hard-sphere simulations that are done in the limit G → ∞ . In this caseonly Γ → ∞ avoids partial crystallization. .50 0.52 0.54 0.56 0.58 0.60 0.62 0.643.54.04.55.05.56.06.5 G-line (c)(b) J-Point unjammed ii > , < << , < , , Z i unjammed C-PointL-Point (a) Figure 8: Schematic representation of the protocol used to dynamically achieved the packingsin the phase diagram. (a) The region near the C-point (high friction, high volume fraction)is the most difficult to access. Here, we prepare an unjammed sample with no friction at φ i (=0.62 in the diagram) very close to RCP, and then we switch friction to compress to therequired frictional jammed state. The packing follow the paths indicated in the figure fortwo friction coefficients. On the other hand, the RLP-line is easily obtained by compressingunjammed packings with low φ i (=0.52 in the diagram) at (b) a given compression rate Γ fora given friction coefficient. (c) There is a small dependence on Γ, the faster the compressionthe deeper we enter in the phase diagram. -1 ( K P a ) time ( MD steps ) Stable packing Unstable packing Compression Z time ( MD steps ) Figure 9: Time evolution of stress (the pressure in the system) for two packings simulated asexplained in the text. The solid red line represents a packing with φ > φ c and dotted blackline represents a packing with φ < φ c , where φ = φ + 2 × − . The inset shows the timeevolution of the coordination number. φ = ∆ φ / ∆ φ = ∆ φ / φ = φ φ = φ − ∆ φ / φ − ∆ φ / φ − ∆ φφ StableUnstableStable? yesno
Figure 10: Flow chart for searching procedure of critical volume fraction φ c . .04.44.85.25.66.0 (b) Z (a) Z Figure 11:
Phase diagram of jammed matter: Simulations . Numerical simulationsdemonstrate how to dynamically access the theoretically found states. The numerical protocolis parameterized by ( φ i , Γ , η, µ ). (a) The plot shows the dependence of the final jammed states( φ, Z ) on φ i for a fix Γ = 10 − and η = 10 − (except for the orange ◮ curve which is for η = 10 − ) from left to right φ i = 0 .
40 (black (cid:4) ), 0.53 (red • ), 0.55 (violet N ), 0.57 (blue H ), 0.59 (green (cid:7) ), 0.61 (pink ◭ ) and 0.63 (orange ◮ ). Each equal-color line set represents adifferent φ i and the dashed lines join systems with the same friction µ . Solid lines represent thetheoretical results with h z = e − or z c = 0 .
01 for different compactivities measured in unitsof 10 − V g . From left to right X = ∞ (black solid line), 1.62 (blue solid line), 1.38 (green solidline), 1.16 (pink solid line) and 0.88 (orange solid line). The error bars correspond to the s.d.over 10 realizations of the packings. (b) The plot focuses on the dependence of ( φ, Z ) on (Γ , η )for two different φ i . Solid black symbols are for φ i = 0 .
40 and (Γ , η ) = (10 − , − ) (black (cid:4) ) and (10 − , − ) (black • ). Open red symbols are for φ i = 0 .
63 and (Γ , η ) = (10 − , − )(red (cid:3) ) and (10 − , − ) (red ◦ ). igure 12: Measuring the temperature of sand. Pictorial representation.Figure 13: Sketch depicting a possible ABC experiment to test the zero-th law of granularthermodynamics and the validity of the phase diagram for finite compactivities. For interme-diate frictions, the A and B systems should follow the arrows to C when put into contact andallow to interchange volume by gently shaking. -4 -3 -2 -1 g z ( r ) r = r - 0.5 RCP linepressure = 100K = 0.000001 = 0.0001 = 0.01 = 0.025 = 0.04 = 0.07 = 0.1 = 0.15 = 0.2 = 0.3 = 0.5 = 1.0 = 10 = 1000 = 100000 g z ( ) =Z g z ( ) >z Figure 14: g z (∆ r ) of packings with various friction coefficient µ along RCP line. We set2 R = 1. -4 -3 -2 -1 g z ( r ) r = r - 0.5 RLP linepressure = 100K = 0.000001 = 0.0001 = 0.01 = 0.025 = 0.04 = 0.07 = 0.1 = 0.15 = 0.2 = 0.3 = 0.5 = 1.0 = 10 = 1000 = 100000 g z ( ) =Z g z ( ) >z Figure 15: g z (∆ r ) of packings with various friction coefficient µ along RLP line. .5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.50123 P ( h z i l ) h z i l RCP linepressure = 100K = 0.000001 = 0.0001 = 0.01 = 0.025 = 0.04 = 0.07 = 0.1 = 0.15 = 0.2 = 0.3 = 0.5 = 1.0 = 10 = 1000 = 100000 ( a ) r = 0 l = 2 l = 2 ( b ) P ( h z i l ) h z i l l = 4 ( c ) P ( h z i l ) h z i l Figure 16: PDF of the coarse-grained coordination number h z i l for packings with variousfriction coefficient µ along the RCP line. (a) ∆ r = 0 and l = 2; (b) ∆ r = 0 .
04 and l = 2; (c)∆ r = 0 .
04 and l = 4. .5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.50123 ( a ) P ( h z i l ) h z i l RLP linepressure = 100K = 0.000001 = 0.0001 = 0.01 = 0.025 = 0.04 = 0.07 = 0.1 = 0.15 = 0.2 = 0.3 = 0.5 = 1.0 = 10 = 1000 = 100000 r = 0 l = 2 ( b ) P ( h z i l ) h z i l r = 0.04 l = 2 Figure 17: PDF of the coarse-grained coordination number h z i l for packings with variousfriction coefficient µ along the RLP line. (a) ∆ r = 0 and l = 2; (b) ∆ r = 0 .
04 and l = 2. z=Z=4 J-Pointz=Z=6 G-line4 z
6, Z=4 Z along RLP line z=Z along RCPlinez=64 Z Figure 18: Summary of the theoretical findings regarding the range of z and Z along thedifferent iso- z , iso- Z , and iso- X lines and the J, C, and L-points in the phase diagram. igure 19: Sketch depicting a possible extension of the phase diagram of Fig. 4 in the( φ, Z ) isostatic plane at the jamming transition to the full phase diagram in the ( φ, Z, σ )space to investigate the criticality of the jamming transition for all packings with any frictioncoefficient. c N i = 0.400, = 0, = 1e-6, = 1e-3 i = 0.630, = 0, = 1e-3, = 1e-4 i = 0.525, = 0, = 2e-9, = 1e-3 ii