Theory of random packings

0.28 MB·

by  Chaoming Song

Theory of random packings
Chaoming Song, Ping Wang, Herna´n A. Makse
LevichInstituteandPhysicsDepartment,CityCollegeofNewYork,NewYork,NY10031,US
0
1
0 Abstract. Wereviewarecentlyproposedtheoryofrandompackings.Wedescribethevolumefluctuationsinjammedmatter
2 throughavolumefunction,amenabletoanalyticalandnumericalcalculations.Wecombineanextendedstatisticalmechanics
approach’alaEdwards’(wheretheroletraditionallyplayedbytheenergyandtemperatureinthermalsystemsissubstitutedby
n
thevolumeandcompactivity)withaconstraintonmechanicalstabilityimposedbytheisostaticcondition.Weshowhowsuch
a
approachescanbringresultsthatcanbecomparedtoexperimentsandallowforanexploitationofthestatisticalmechanics
J
framework.ThekeyresultistheuseofarelationbetweenthelocalVoronoivolumeoftheconstituentgrainsandthenumber
9
ofneighborsincontactthatpermitsasimplecombinationofthetwoapproachestodevelopatheoryofrandompackings.We
2
predict thedensityof randomloosepacking (RLP)andrandomclosepacking(RCP)incloseagreement withexperiments

anddevelopaphasediagramofjammedmatterthatprovidesaunifyingviewofthedisorderedhardspherepackingproblem
]
andfurthersheddinglightonadiversespectrumofdata,includingtheRLPstate.Theoreticalresultsarewellreproducedby
t
f numericalsimulationsthatconfirmtheessentialroleplayedbyfrictionindeterminingboththeRLPandRCPlimits.Finally
o
wepresentanextendeddiscussionontheexistenceofgeometricalandmechanicalcoordinationnumbersandhowtomeasure
s
bothquantitiesinexperimentsandcomputersimulations.
.
t
a Keywords: granularmatter,randomclosepacking,statisticalmechanics
m PACS: 81.05.Rm,83.80.Fg
-
d
n I.STATISTICAL MECHANICSOF JAMMED MATTER
o
c
ConventionalStatisticalMechanicsusestheergodichypothesistoderivethemicrocanonicalandcanonicalensembles,
[
based on the quantities conserved, typically the energy E [1]. Thus the entropy in the microcanonicalensemble is

1 S(E)=k log d (E H(p,q))dpdq,whereH(p,q)istheHamiltonian.Thisbecomesthecanonicalensemblewith
B
v exp[ H(¶ S/R¶ E)].−Experiments [2, 3, 4, 5] indicate that systematically shaken granular materials show reversible
8 −
behavior, and the analogue of the conserved quantity is the volume V, thus the micro-canonical ensemble or V-
6
4 ensembleis[6,7,8,9,10,11,12]:
5
1. W (V)=exp[S(V)/l ]=Z d (V−W(~ri))Q jam(~ri)d~ri, (1)
0
0 whereQ jam(~ri)isafunctionthatdefinesthejammedconfiguration.Asaminimumrequirementthejammingfunction
1 Q jam(~ri) should ensure touching grains, and obedience to Newton’s force laws.~ri denotes the particle positions in
: the system and W(~r) is the volume function defining the volume associated with each grain (see below). This
v i
givesa canonicalensembleofexp[ W(¶ S/¶ V)]. Just as¶ E/¶ S=T isthe temperaturein equilibriumsystem, the
i
X temperature-like variable in granula−r systems is the compactivity X =¶ V/¶ S. In Eq. (1), l is the analogue of the
r Boltzmannconstant.
a
Thermodynamicanalogiesmayilluminatemethodsforattemptingtosolvecertainproblems,butinevitablyfailat
somepointintheirapplication.Themodeofthisfailureisaninterestingphenomenon,illustratedbythecompaction
experimentsof the groupsof Chicago, Texas, Paris and Schlumberger[2, 3, 4, 5]. They have shown that reversible
states exist along a branch of compaction curve where statistical mechanics is more likely to work. Conversely,
experimentsalso showeda branchof irreversibilitywhere the statistical frameworkis notexpectedto work. Poorly
consolidated formations, such as a sandpile, are irreversible and a new “out-of-equilibrium” theory is required to
describethem.Belowwefocusonatheoreticaldescriptionofthereversiblebranchofthecompactioncurvefocusing
onatheoryoftherandomclosepackedstate.
ThecanonicalpartitionfunctionintheV-ensembleisthestartingpointofthestatisticalanalysisofjamming:
Z(X)= g(W)e−W/X Q jamdW, (2)
Z
whereg(W)isthedensityofstatesforagivenvolumeW.
From Eq. (2) we identify three minimal steps in developing analytical solutions which are discussed in the next
sections.SectionIIdiscussestheneedforavolumefunctionintermsofthecontactnetwork.SectionIIIdiscussesthe
needforaproperdefinitionofjammedstatethatallowsonetodefineQ .SectionIVdiscussesthedensityofstates.
jam
FinallyinSectionV we explainthegeometricalandmechanicalcoordinationnumbersandhowtomeasurethemin
SectionV.A,andweconcludeinSectionVI.
II. VOLUMEFUNCTION
Whileitisalwayspossibletoquantifythetotalvolumeofthesystem,itisunclearhowtotreatthevolumefluctuations
atthegrainlevel.ThefirststeptostudytheV-ensembleistofindthevolumeW(~r)associatedtoeachparticle~r that
i i
successfullytilesthesystem.Thisisanalogoustotheadditivepropertyofenergyinequilibriumstatisticalmechanics.
Initial attempts included a model volume function under mean-field approximation [6], the work of Ball and
Blumenfeld[13] andsimplerversionsintermsofthefirstcoordinationshell[14].Thesedefinitionsareproblematic
since some are not additive, others present problems in polydisperse systems or are proportional to coordination
contrarytoexpectation.InRef.[15,16]wehavefoundananalyticalformofthevolumefunctioninthree-dimensions
anddemonstratedthatitistheVoronoivolumeofaparticlei:
1 r 3
Wvor= min ( ij ) ds, (3)
i 3I (cid:18)sˆrˆij>0 2sˆ rˆij (cid:19)
· ·
where~r isthevectorfromthepositionofparticleitothatofparticle j,theintegrationisdoneoverallthedirections
ij
sˆforminganangleq with~r asinFig.1a,andRistheradiusofthegrain.Rwillbesettounityforsimplicity.While
ij ij
thisformulamayseemcomplicated,ithasasimpleinterpretationdepictedinFig.1a.
The Voronoi construction is additive and successfully tiles the total volume. Prior to this result, there was no
analyticalformulatocalculatetheVoronoivolumeintermsofthecontactnetworkr .Afurthersimplificationarises
ij
whenweconsiderisotropicsystems.Thenthevolumefunctionreducestotheorientationvolume,withouttheaverage
oversˆ.Wedefinethereducedfreeorientationalvolumefunctionas
Ws V
ws i − g, (4)
≡ Vg
withWs V 1 min rij 3,seeFig.1a(V istheparticlevolume).Thisequationallowstheoreticalanalysisin
i ≡ g(cid:16)2R sˆ·rˆij>0sˆ·rˆij(cid:17) g
theV-ensemblesinceitreducesthecomplicateddefinition(3)toamoreamenable“one-dimensional”volumewhich
canbetreatedanalytically.
ThenextstepistodevelopatheoryofvolumefluctuationstocoarsegrainWs overamesoscopiclengthscale.We
i
callthisthequasi-particleapproximation.Itcouldbeconsideredaswellasamean-fieldapproximation,althoughmean
fieldissupposedtobeexactininfinitedimensions.Theapproximationsusedinthepresenttheoryaresupposedtoget
betterasthedimensionincreases,butwecannotclaimthatthetheoryisexactinlargedimensions.Thus,wepreferto
callourapproximation“quasi-particle”inthespiritofLandauandthequasiparticlesas“coordinons”.
Thecoarseningreducesthedegreesoffreedomtoonevariable,thecoordinationnumberofeachgrain,anddefines
anaveragevolumefunctionwhichismoreamenabletostatisticalcalculationsthanEq.(3)asshownin[16]:
2√3
w(z)= ws = , (5)
i
h i z
validformonodispersehardsphereswherezisthegeometricalcoordinationnumber.FornowonweassumeV =1
g
for simplicity.The available volumeper grain is inverselyproportionalwith the coordinationnumber,in agreement
withtheX-raytomographyexperiments(seeFig.6in[17]wherethevolumefractionisf 1=w+1).
−
III.DEFINITIONOF JAMMINGVIAQ :ISOSTATIC ENSEMBLE
jam
ThedefinitionoftheconstraintfunctionQ isintimatelyrelatedtotheproperdefinitionofajammedstate,witha
jam
minimumrequirementofmechanicalequilibrium.Inanattempttodefinethejammedstatesinarigorousmathematical
way,Torquatoandcoworkershaveproposedthreecategoriesofjamming[18]:locallyjammed,collectivelyjammed
and strictly jammed based on geometrical constraints. Unfortunately this definition cannot be easily extended to
frictional systems since it is based on geometry and does not include the contact forces. Other approaches based
onminimaofthepotentialenergylandscapealsofailsincesuchapotentialdoesnotexistforfrictionalgrainsdueto
theirpath-dependency.
Thenwedefineanalternativeapproachtocharacterizejammingforthegeneralcaseoffrictionalgranularmatter.In
[16]weproposetheisostaticcondition[19,20]asapossibleformulationofjamming.Theisostaticconditionimplies
a mechanical coordination number to be Z =2d, where d is the dimension, for frictionless spherical particles and
Z=d+1forinfinitelyroughparticles(withinterparticlefrictioncoefficientm ¥ ).Numericalsimulations[21,22],
→
experiments[23]andtheoreticalwork[20]suggestthatatthejammingtransitionthesystembecomesexactlyisostatic.
However, no rigorous proof of this statement exist. It should be noted that Z =6 is a necessary but not sufficient
condition for isostaticity [20]. A kinematic condition has to be satisfied as well, which refers to the location of the
centeroftheparticlesthataredeterminedonlybythelengthofthevectorsthatjointhecenterofcontactingparticles.
In other words the equations of equilibrium needs to be independent of each other. While this condition has been
provedto exist in frictionless packings[20], the problem of frictional packings, even with infinite friction, remains
open. Thus, to be precise, the only point where we can claim isostaticity is the frictionless point. In any case, we
considerthatweareabletoextendtheisostaticconditiontoinfinitefrictionalpackingsaswellatZ=4.Interpolating
betweenthetwolimits,thereexistpackingsoffinitem ;thecoordinationnumbersmoothlyvariesbetweenZ(m =0)=6
andZ(m ¥ ) 4[24,16].
→ →
WenoticethatwhiletherelationshipbetweenfrictionandZmaynotbeunique,thetheoryisonlybasedonZ.Thus,
given the mechanical coordination number we predict the state of the packing with the compactivity. Additionally,
we conjecture that m determines Z(m ) if we follow certain protocols as discussed in [16]. These protocols imply
compressionofapackingfromanunjammedstatetojammingbyfollowingonesinglecontinuouspath.Thatis,the
pathdependentshearforcesarenotresetbyasuddenchangeinthepreparationpath.Itshouldbestatedthatthereare
otherprotocols,like shearcycling,thatstartwith a givenZ andcan producepackingsthatcontinuouslycompactify
untilRCP(andevenbeyond).Thisisdonebyeffectivelychangingthepathfollowedbytheshearforcesateverycycle.
Thus,forshearcyclingprotocols,a uniquerelationmaynotbeexpectedbetweenZ and m .However,thetheoretical
results as derived below are still valid in shear cycling experiment since they concern the relation between Z and
volume.Infact,ashearexperimentmaybetheeasiestwaytoobtainthepackingsattheRLPlineasexplainedbelow.
Assumingthata systemofhardspheresisisostatic atthejammingtransition,Eq.(2)canbewrittenintermsofz
andthemean-fieldEq.(5)canbeusedinthesingle-particle(ormorepreciselyquasi-particle)partitionfunction:
6
Ziso(X)= e−w(z)/Xg(z)dz. (6)
ZZ
Wenoticethatinprinciplethecoordinationnumberofeachparticletakesonlyintegervalues.However,thecoordina-
tionnumberin(6)impliesacoarsegrainingoverseveralparticles,andthereforecantakenon-integervalues.Thusthe
useofanintegralin(6)insteadofasumisjustified
Thisensembleis referredto asthe Isostatic-ensemble.Note thatthe upperlimit of integrationis z=6 [16]. This
impliesthatonly disorderedpackingsare included.The solutionof such a partitionfunctionfor monodispersehard
sphereshasbeendonein[16]revealingthephasediagramdepictedinFig.1b.
Thisphasediagrampredictsaseriesofimportantresults,suchasthevalueofRCPatX =0,
f =6/(6+2√3), (7)
RCP
andthelowestdensityoftheRLPatX =¥ ,
f =4/(4+2√3), (8)
RLP
in close agreement with experiments. The diagram restricts the possible packings to the yellow triangle in Fig. 1b,
rangingfromfrictionlesssystemswithZ=6,toinfinitelyroughgrainsintheZ=4granularlineorG-line.
frictionless point
6.0 (J-point)
5.5 RXL P 1line R liCneP
Z 5.0 C A
X=0
4.5
B C
4.0
G-line
0.54 0.56 0.58 0.60 0.62 0.64
(a) (b) (c)
FIGURE1. (a)SchematicsoftheVoronoivolumeandtheorientationalvolumeassociatedwithparticlei.Theboundaryofthe
Voronoicell(shownintwo-dimensionsforsimplicity)correspondstotheirregularpentagoninblackwhichdefinesWvor.Thelimit
i
oftheVoronoicellofparticleiinthedirectionsˆistheminimumofrij/2cosq ij overalltheparticlesinthepacking,asindicated.
ThisdefinestheorientationalvolumeWis whichisthevolumeofthesphereofradiusrij/2cosq ij definedbythedashredcircle
inthefigure.(b)Phasediagramofjamminginthehardsphereplaneundertheisostaticassumption.Allthedisorderedpackings
fallwithintheyellow triangledemarcated bytheRCPandRLPlinesandtheG-line.Theisocompactivitylinesareincolor.(c)
Generalizationofthephasediagramtothespace(Z,f ,p).
IV.DENSITY OF STATES
A difficult problem is the determination of the density of states g(W) in Eq. (2). For the simplest case of the Iso-
ensemblefromEq(6),thedensityofstatesreducesto
g(z)=(h )z, (9)
z
whereh isasmallmicroscopicconstantarisingduetothediscretevolumespaceofconfigurations[16].Thesituation
z
isanalogoustothediscretenessoftheconfigurationspaceimposedbytheHeisenberguncertaintyprincipleinquantum
mechanics.The formulais analogousto thefactorh d forthe densityofstates in equilibriumstatistical mechanics.
−
Whilethedegreesoffreedom p,q arecontinuous,theuncertaintyprincipleimposesthediscreteness(D p,D q)in
i i
theconfigurationalspacegiven{byD p}D q h.Thisconsiderationallowsfortheapproximatesolutionexplainedinthe
∼
abovesectionanddepictedinthephasediagramofFig.1b.
V. GEOMETRICALVERSUS MECHANICALCOORDINATIONNUMBER
ItisimportanttonotethatthederivationofthevolumefunctioninSectionIIimpliesnothingaboutthevalueofthe
contactforces;thevolumefunctionrepresentsthecontributionarisingpurelyfromthegeometryofthepacking.Thus,
the coordinationnumberz appearingin Eq. (5) is the geometrical coordinationnumberrelated to volume,which is
differentfromthe mechanicalcoordinationnumberZ that countsthe numberof contactsper particle with non-zero
forcerelatedtotheisostaticconditionandforcenetwork.
Havingacknowledgeda differencebetweenthe geometricalcoordinationnumberz in Eq.(5)andthe mechanical
coordinationnumberZ which countsonlythe contactswith non-zeroforces,below we discuss the boundsof z and
howtomeasureit.
Sincesomegeometricalcontactsmaycarrynoforce,thenwehave:
Z z. (10)
≤
Toshowthis,imagineapackingofinfinitelyrough(m ¥ )sphereswithvolumefractioncloseto0.64.Theremust
→
be z=6 nearestneighborsaroundeachparticleon theaverage.However,the mechanicalbalancelaw requiresonly
Z=4contactsperparticleonaverage,implyingthat2contactshavezeroforceanddonotcontributetothecontact
forcenetwork.
SuchasituationispossibleasshowninFig.2:startingwiththecontactnetworkofanisostaticpackingoffrictionless
spheres having Z =6 and all contacts carrying forces (then z=6 also as shown in Fig. 2a), we simply allow the
existence of tangential forces between the particles and switch the friction coefficient to infinity. Subsequently, we
solve the force and torque balance equations again for this modified packing of infinitely rough spheres but same
Isostatic packing
µ=0 z=6 Z=6 µ=∞ z=6 Z=4
(a) (b)
FIGURE2. (a)Considerafrictionlesspackingattheisostaticlimitwithz=6.Inthiscasetheisostaticconditionimpliesalso
Z=6mechanicalforcesfromthesurroundingparticles.(b)Ifwenowswitchonthetangentialforcesusingthesamepackingasin
(a)bysettingm ¥ ,theparticlerequiresonlyZ=4contactstoberigid.Suchasolutionisguaranteedbytheisostaticcondition
form ¥ .Thus→,theparticlestillhavez=6geometricalneighborsbutonlyZ=4mechanicalones.
→
geometricalnetwork,asshownin Fig.2b[Noticethattheshearforceiscomposedofanelastic Mindlincomponent
plustheCoulombconditiondeterminedbym .Thuswhenm ¥ ,theelasticMindlincomponentstillremains].
→
Theresultingpackingismechanicallystableandisobtainedbysettingtozerotheforcesoftwocontactsperball,on
average,tosatisfythenewforceandtorquebalanceconditionfortheadditionaltangentialforceatthecontact.Such
asolutionisguaranteedtoexistduetotheisostaticcondition:atZ=4thenumberofequationsequalsthenumberof
forcevariables.Despitemechanicalequilibrium,givingZ=4,therearestillz=6geometricalcontactscontributing
tothevolumefunction.
Therefore,weidentifytwotypesofcoordinationnumber:thegeometricalcoordinationnumber,z, contributingto
thevolumefunctionandthemechanicalcoordinationnumber,Z, measuringthecontactsthatcarryforcesonly.This
distinctioniscrucialtounderstandthesumoverthestatesandtheboundsinthepartitionfunction.
We have established a lower bound of the geometrical coordination in Eq. (10). The upper bound arises from
consideringtheconstraintsinthepositionsoftherigidhardspheres.Forhardspheres,theNdpositionsoftheparticles
areconstrainedbytheNz/2geometricalconstraints, ~r =2R,ofrigidity.Heredisthedimension.Thus,thenumber
ij
| |
ofcontactssatisfiesNz/2 Nd,andzisboundedby:
≤
z 2d. (11)
≤
Noticethatthisupperboundappliestothegeometricalcoordination,zandnottothemechanicalone,Z,anditisvalid
foranysystemirrespectiveofthefrictioncoefficient,fromm =0 ¥ .
Inconclusion,themechanicalcoordinationnumber,Z,rangesf→rom4to6asafunctionofm ,andprovidesalower
bound to the geometrical coordinationnumber,while the upper bound is 2d. A granular system is specified by the
interparticlefrictionwhichdeterminestheaveragemechanicalcoordinationatwhichthesystemisequilibrated,Z(m ).
ThepossiblemicrostatesintheensembleavailableforthissystemfollowaBoltzmanndistributionEq.(2)forstates
satisfyingthefollowingbounds:
Z(m ) z 2d=6. (12)
≤ ≤
V.A.How tomeasure the geometrical coordinationnumber
Measuringthegeometricalcoordinationnumbercanbeatrickytask,inprinciple.Attheonset,itisthecoordination
numberofasinglequasiparticle.Whatwemeasureinarealpacking(numericallyorexperimentallygenerated)isan
ensemble average of many quasiparticle accordingto the partition function Eq. (6). Thus, rigorouslyspeaking, it is
notpossibletoisolateaquasi-particleandmeasureitspropertiesinarealpacking.Beyondthiscausticandsomehow
pessimisticremark,yetrigorous,belowweofferlightattheendofthetunnelbyusingthetheoreticalpredictionsto
defineanapproximative,yetaccurate,waytomeasurethegeometricalcoordination.
Figure3summarizesthepredictionsofthetheoryregardingthebehaviourofzandZ forthepackingsinthephase
diagram.First,wecanthinkthatthequasiparticlebehaviourisrevealedwhenthesystemhasinfinitecompactivityand
6.5
J-Point
z=Z=6
6.0
5.5 along RLP line along
z=Z RCP
Z line
5.0 iso-Z
z=6
4.5 iso-z 46Z66
4.0
L-Point G-line C-Point
z=Z=4 46z66, Z=4 z=6, Z=4
3.5
0.52 0.54 0.56 0.58 0.60 0.62 0.64 0.66
FIGURE3. SummaryofthetheoreticalfindingsregardingtherangeofzandZalongthedifferentiso-z,iso-Z,andiso-X lines
andtheJ,C,andL-pointsinthephasediagram.
behaves like a non-interacting“gas” of quasiparticles. This is the behaviour found along the RLP line in the phase
diagramofFig.1b.Indeed,alongthisline,thepackingshavethehighestentropy[25]andthereforetheyarethemost
likelytobefoundnumericallyandexperimentally.Indeed,thisiswhatisfoundinournumericalsimulations[16,25].
Alongthislinewe expectthatthemechanicalandgeometricalcoordinationsarethesame,z Z.Thisresultcomes
≈
aboutsinceatinfinitecompactivityweareexploringthestateswiththehighestvolumeorlowervolumefraction.In
the partitionfunctionEq.(6), these are states with z Z. Itshouldbe noticedthat the densityof states playsa role
≈
in this argument. We are assuming that the density decays very rapidly owning to a small constant h 0 in Eq.
z
→
(9). This is of course a very reasonable assumption, since h is related to the discretization of the volume space of
z
configurations,anditislikeaPlanckconstantofgranularmatter.However,wemayevenrelaxthisconsiderationand
allow this constantto be of the order 1. The aboveargumentis still correct, but the only differenceis that the RLP
limit is notat infinite compactivitybut in the limit of X 0 . Thatis, the highest volumecorrespondsto negative
−
temperaturestates,whichareevenhotterthanX +¥ .A→nextendedexplanationofthispointisgivenin[25]andits
→
followuppaperincond-mat.
The second interesting aspect of the theory is the prediction that along the RCP line in Fig. 3, the geometrical
coordination number is constant and equal to the maximum coordination, which produces the minimum Voronoi
volume:z=6.Atthesametime,themechanicalcoordinationnumbervariesfrom4to6aswereducefrictiontozero.
Thisresultis explainedwith theanalysisofferedinFig. 2andexplainedabove.The conclusionis thatthe packings
alongtheRCPlinearegeometricallythesame(withz=6allofthem)buttheydifferinthevalueoftheforcesbetween
theparticlesasZ variesfrom6to4.Thispredictioncanbealsoappliedtootherpackingswithlowerf .Foragiven
volumefraction,thepackingsalongtheiso-zverticallineinFig.3shouldhaveapproximatelythesamezbutdifferent
Zfrom4tothemaximumZgivenbytheRLPline.
Using this theoretical result, it is easy then to define the geometrical coordination number for real packings and
propose a clear way to measure it. The idea is to inflate the particles infinitesimally by a D r value and measure
the contacting particles. By setting D r = 0, we clearly measure the mechanical coordination. By considering an
infinitesimally small D r we should measure the geometrical one. The question is to know what value of D r to use.
Hereiswherethetheorycomeshandy.WeknowthatalongtheRCPlinethepackingsarethesamegeometrically.So,
whateverthedefinitionofcoordinationnumberweuse,itshouldsatisfythatafteragivenD r,weshouldfindthesame
packingstructureforthepackingsalongtheRCPlineandbelowthegivenD r,thestructureshouldchangereflecting
thedifferentvaluesofZ fordifferentpackingsintheRCPline.
Using these considerations, we identify the geometrical coordination as follows. Two particles that may not be
in contact (giving rise to a zero force) may be close enough to be considered as contributing to the geometrical
coordination.Indeed,itisknownthattheradialdistributionfunctiong(r)hasasingularity,g(r) (r 0.5) 0.5[24],
−
∼ −
implyingthat there are manyparticlesalmost touching.We introducea modified radialdistributionfunction(RDF)
8.0 RCP line
pressure = 100K
7.5 gz(0.04) >z === 000...000001000101
7.0 = 0.025
gz(0)=Z = 0.04
6.5 == 00..017
= 0.15
)r 6.0 == 00..23
= 0.5
(gz5.5 == 11.00
= 1000
5.0 = 100000
4.5
4.0
3.5
3.0
-4 -3 -2 -1 0
10 10 10 10 10
r = r - 0.5
FIGURE4. gz(D r)ofpackingswithvariousfrictioncoefficientm alongRCPline.Weset2R=1.
8.0 RLP line
pressure = 100K
7.5 = 0.000001
gz(0.04) >z == 00..000101
7.0 gz(0)=Z == 00..00245
6.5 == 00..017
= 0.15
)r 6.0 == 00..23
= 0.5
(gz5.5 == 11.00
= 1000
5.0 = 100000
4.5
4.0
3.5
3.0
-4 -3 -2 -1 0
10 10 10 10 10
r = r - 0.5
FIGURE5. gz(D r)ofpackingswithvariousfrictioncoefficientm alongRLPline.
g (r)inordertoapproximatelyidentifyzandZfromrealpackings:
z
g (r)= 1 R2(cid:229)N (cid:229)N Q rij 1 Q r+R 1 , r>R (13)
z N r2 i j=i (cid:16)r−R− (cid:17) (cid:16) rij − (cid:17)
6
whereRistheradiusofparticle,Nisthenumberofparticles,r isthedistanceoftwoparticle’scenters,r = ~r ~r ,
ij ij i j
andQ istheHeavisidestepfunction.TheRDFdescribestheaveragevalueofthenumberofgrainsincontac|tw−ith|a
virtualparticlewhichhasbeeninflateduptoaradiusr R,andthefactorofR2/r2 istheratioofarealsphere’sarea
andthevirtualone’s.Withoutthenormalizationfactor≥R2/r2,Eq.(13)isthesamedefinitionofcoordinationasused
byTorquatoandZamponiintheiranalysisofinfinitepressurejammedhardsphereglasses[26,27].Thisfactorisnot
crucialforouranalysis(asthemainconstraintisthatwhenr=RweshouldgetZ)asexplainedbelow,butweargue
thatitisusefulsinceweneedtopropernormalizebythefactthatweareinflatingtheballs.PleasenoticethatEq.(29)
inSupplementaryInformationSectionof[16]containsatypo.ThecorrectdefinitionisEq.(13).
g (r)measuresthenumberofballswiththeirvolumeintersectingthesurfaceofasphereofradiusrmeasuredfrom
z
thecenterofagivenball.Whenr=Rin(13)weobtainthemechanicalcoordinationnumberwhilethegeometrical
oneisobtainedforasmallvalueD r= r R =0forwhichwedistinctlyfinda signaturefromcomputersimulations,
unambiguouslydefiningitatD r=0.04f2−oRrt6hesystemsizeusedbyfollowingthepackingsalongtheRCPline.
Figures 4 and 5 plot the g (D r) of packings with various friction coefficient m along the RCP and RLP lines
z
respectively. Following the definition of Eq. (13), g , with D r = 0, should be directly equal to the mechanical
z
coordinationnumber,Z,andshouldrangefrom4to6alongbothRCPandRLP(ifh 1)lineswhichisconfirmed
z
≪
byournumericalsimulationsinFigs.4and5,respectively.
More importantly,as shown in the figures, we find thatg (D r) along the RCP line is exactly the same for all the
z
packingswhenD r>0.04asshowninFig.4.ForD r<0.04,Fig.4showsthatthepackingshavedifferentmechanical
coordinationnumbersfrom4to6.Thisisexactlywhatthetheorypredicted.Allthesepackingshaveactuallythesame
geometricalstructureevidencedwhenD r>0.04butwithdifferentmechanicalcoordinationswhichappearsonlyina
differenceobservedforD r<0.04.Basedonthisanalysis,wethendefinethegeometricalcoordinationnumberasthe
oneappearingatD r=0.04.We identifythegeometricalcoordinationnumberasz=g (0.04)undertheaccuracyof
z
thesimulationsandforthisparticularsystemsizeN=10,000(wenoticethoughthatthisvaluemaydependonsystem
size).
Itisimportanttonotethatintermsoftheradialdistributionfunction,g(r),nothingreallyhappensatD r=0.04:that
is,thereisnopeaking(r)at0.04andthesecondpeakafterthefirstcoordinationshellappearsforlargerD ring(r).
On the otherhand,we clearly see a peakin Fig. 4 at thisvalue. We pointoutthatthe peakat 0.04is the byproduct
of the normalizationfactor R2/r2 in Eq. (13). This factor is suggested since we need to renormalize by the area of
the virtualsphere. However,other factors, for instance R3/r3 in Eq. (13)would producea peak in Fig. 4 locatedin
anotherposition.Theimportantfactisnotthelocationofthepeak,butthefactthatabove0.04allthefunctionsinFig.
4coincide.Thisisthebasisofthedefinitionof0.04asthelocationtodefinethegeometricalcoordinationnumber.At
thisposition,thereisnopeakintheg(r).Indeed,thesecondpeaking(r)beyondthefirstcoordinationshellappears
much further around D r 1 and are indeed also identified by g (D r) as can be seen in Fig. 4 and Fig. 5 as well.
z
≈
Thereforethepeakassociatedwiththegeometricalcoordinationnumberisnotrevealedfromthestructureing(r).It
hasamoresubtlemeaningasexplainedabove.
It is also important to note that in experiments there is always an uncertainty in measuring the position of the
particles.Accordingtoouranalysisansmalluncertaintyofafewpercentwillrenderthemechanicalcoordinationinto
thegeometricalone.Thusexperimentswillbe verydifficultto differentiatebetweenzandZ. A possiblesolutionto
thisproblemistousecomplementaryfluorescenttechniques[23,28, 29]toobtaineda signalwhentheparticlesare
carryingaforceandnottorelyongeometricalreconstructionofindex-matchedimagesorX-raytomography.Another
route would be to obtain approximate coordinates from experiments and then use them as input into a Molecular
Dynamicssimulationstoobtaintheexactforcebalanceforeachparticle.Thislastapproachmayprovidethefinalway
toaccuratelymeasurethecoordinationsoftheparticleswithaccuracy.Inourwebsitehttp://www.jamlab.orgweoffer
thecomputercodestoperformMDsimulationswithHertz-Mindlinforcesaswellthecodetocalculatetheentropyof
thepackings.
The theoretical analysis is also confirmed in the RLP packings. Along the RLP line, Fig. 5, we find that the
geometricalcoordinationnumberasextractedfromg (0.04)isveryclosetothemechanicalone.SincetheRLPline
z
isatX ¥ andh 1,allthestatesalongRLPhavez Zaswemovealongthelinevaryingthefrictioncoefficient.
z
→ ≪ ≈
Thus,thenumericalresultsconfirmthetheory.
Inconclusion,aprescriptiontomeasurethegeometricalcoordinationisthefollowing:Firstweidentifyatheoretical
waytodefineit.Forinstance,hereweusethetheoreticalpredictionthatallthepackingsalongRCPhavethesamez.
Withtheproperdefinitionofcoordinationnumberforaninflatedparticle,Eq.(13),wecalculatethecoordinationasa
functionofD r.ThisidentifiesD r=0.04asthepositiontoobtainthegeometricalz.Wethenexploreanypacking(not
only at RCP line) and apply Eq. (13) at D r=0.04 and obtain z. Figure 3 summarizesthe theoreticalpredictionsof
valuesofzandZforallthepackingsinthephasediagram.Ifthereaderstillhasdoubtsaboutthedifferencebetweenz
andZweofferafinalmorevividwaytounderstanditintermsofthefamouskissingnumberconjecturefromNewton
andGregoryduetoaremarkofA.Coniglio(privatecommunication):amechanicalcontactislikeaFrenchkisswhile
ageometricalcontactisanyotherinconsequentialkiss.
VI. CONCLUSIONS
Inconclusion,usingEdwardsstatisticalmechanicswehaveelucidatedsomeaspectsofRLPandRCPinthedisordered
spherical packing problem. The phase diagram introduced here serves as a beginning to understand how random
packingsfillspaceinthreedimensions.Thecomparativeadvantageofthepresentapproachoverextensiveworkdone
in the past, is in the classification of all packings throughX, Z and f in the theoretical phase diagram from where
thesestudiescouldbesystematicallyperformed.Thisclassificationguidesthesearchforindicationsofjammingfrom
asystematic pointofview,throughtheexplorationofalljammedstatesfrom m =0to m ¥ . Ourresultsnotonly
→
applytopackingsatthejammingtransitioninthelimitofhardspheres,butmayalsobeextendedtoageneralphase
diagramassketchedinFig.1ctoincludestateswithfinitenonzeropressure.Suchstatesaredescribedbyanangoricity
in addition to the compactivity as developed here. Extensions to other dimensions, polydisperse systems and other
shapesofparticleslike ellipsoidsandsphero-cylindersarebeingworkedoutwith the goalof developinga unifying
thermodynamicviewofthephysicsofpackings.
ACKNOWLEDGMENTS
ThisworkissupportedbytheNationalScienceFoundation,CMMTDivisionandtheDepartmentofEnergy,Officeof
BasicEnergySciences,GeosciencesDivision.
REFERENCES
1. L.D.LandauandE.M.Lifshitz,StatisticalPhysics(Pergamon,NY,1970).
2. E.R.Nowak,J.B.Knight,M.L.Povinelli,H.M.JaegerandS.R.Nagel,PowderTechnol.94,79(1997).
3. P.Philippe,andD.Bideau,Europhys.Lett.60,677(2002).
4. J.Brujic´,P.Wang,C.Song,D.L.Johnson,O.Sindt,andH.A.Makse,Phys.Rev.Lett.95,128001(2005).
5. M.Schro¨ter,D.I.Goldman,andH.L.Swinney,Phys.Rev.E71,030301(R)(2005).
6. S.F.Edwards,inGranularmatter:aninterdisciplinaryapproach(edA.Mehta)121-140(Springer-Verlag,NewYork,1994).
7. R.BlumenfeldandS.F.Edwards,Phys.Rev.Lett.90,114303-1(2003)
8. F.Lechenault,F.daCruz,O.Dauchot,andE.Bertin,J.Stat.Mech.,P07009(2006).
9. M.P.Ciamarra,A.Coniglio,andM.Nicodemi,Phys.Rev.Lett.97,158001(2006).
10. H.A.MakseandJ.Kurchan,Nature415,614(2002).
11. P.Wang,C.Song,andH.A.Makse,NaturePhysics2,526-531(2006).
12. C.Song,P.Wang,andH.A.Makse,Proc.Nat.Acad.Sci.102,2299-2304(2005).
13. R.C.BallandR.Blumenfeld,Phys.Rev.Lett.88,115505-1(2002).
14. H.A.Makse, J.Brujic´,andS.F.Edwards,inThePhysicsofGranular Media,editedbyH.HinrichsenandD.E.Wolf
(Wiley-VCH,2004).
15. C.Song,P.Wang,andH.A.Makse,http://arxiv.org/abs/0808.2186
16. P.Wang,C.Song,andH.A.Makse,Nature453,629-632(2008).
17. T.Aste,M.Saadatfar,andT.J.Senden,J.Stat.Mech.,P07010(2006).
18. S.TorquatoandF.H.Stillinger,J.Phys.ChemB105,11849(2001).
19. S.Alexander,Phys.Rep.296,65(1998).
20. C.F.Moukarzel,Phys.Rev.Lett81,1634(1998);ibid,Gran.Matter3,41(2001).
21. H.A.Makse,D.L.Johnson,andL.M.Schwartz,Phys.Rev.Lett.84,4160(2000).
22. C.S.O’Hern,S.A.Langer,A.J.Liu,andS.R.Nagel,PhysRev.Lett.88,075507(2002).
23. J.Brujic´,C.Song,P.Wang,C.Briscoe,G.Marty,andH.A.Makse,Phys.Rev.Lett.98,248001(2007).
24. L.E.Silbert,D.Ertas,G.S.Grest,T.C.Halsey,andD.Levine,Phys.Rev.E65,031304(2002).
25. C.Briscoe,C.Song,P.Wang,andH.A.Makse.Phys.Rev.Lett.101,188001(2008).
26. M.Skoge,A.Donev,F.H.Stillinger,andS.Torquato,Phys.Rev.E74,041127(2006).
27. G.Parisi,andF.Zamponi,tobepublishedinRev.Mod.Phys.Arxivpreprintatwww.arxiv.org/abs/0802.2180(2008).
28. J.Brujic´,S.F.Edwards,D.V.Grinev,I.Hopkinson,D.Brujic´,andH.A.Makse,FaradayDiscuss.,123,207(2003).
29. J.Brujic´,S.F.Edwards,I.Hopkinson,H.A.Makse,PhysicaA327,201(2003).