Jamming in two-dimensional packings
Sam Meyer, Chaoming Song, Yuliang Jin, Kun Wang, Hernán A. Makse
aa r X i v : . [ c ond - m a t . s o f t ] N ov Jamming in two-dimensional packings
Sam Meyer a , Chaoming Song b , Yuliang Jin b , Kun Wang b , Hern´an A. Makse b, ∗ a ´Ecole Normale Sup´erieure de Lyon, Universit´e Claude Bernard, Lyon I, France 695014 b Levich Institute and Physics Department, City College of New York, New York, NY, USA10031
Abstract
We investigate the existence of random close and random loose packing limits intwo-dimensional packings of monodisperse hard disks. A statistical mechanicsapproach— based on several approximations to predict the probability distribu-tion of volumes— suggests the existence of the limiting densities of the jammedpackings according to their coordination number and compactivity. This resulthas implications for the understanding of disordered states in the disk packingproblem as well as the existence of a putative glass transition in two dimensionalsystems.
Keywords:
Random packings; Volume function; Disordered system
1. Introduction
The concept of jamming is a common feature of out of equilibrium systemsexperiencing a dynamical arrest ranging from emulsions, colloids, glasses andspin glasses, as well as granular materials [1]. For granular matter, it is arguedthat a statistical mechanical description can be used, with volume replacing en-ergy as the conservative quantity [2]. In this framework, a mesoscopic model hasbeen presented [3], allowing the development of a thermodynamics for jammingin any dimension. Here we develop this theoretical approach to investigate theexistence of disordered packings in two-dimensional systems composed of equal-sized hard disks [4]. The existence of amorphous packings in 2d is a problem ofdebate in the literature: two dimensional systems are found to crystallize veryeasily since disordered packings of disks are particularly unstable [5].In two dimensional Euclidean space, the hexagonal packing arrangement ofcircles (honeycomb circle packing) has the highest density of all possible planepackings (ordered or disordered) with a volume fraction φ hex = π √ ≃ . φ rcp ≈ .
82 [4] while the lower limit (random loose packing, RLP) has ∗ Corresponding author
Email address: [email protected] (Hern´an A. Makse)
Preprint submitted to Elsevier November 18, 2010 een little investigated and its existence has not been treated so far. The ten-dency of 2d packings to easily crystallize has led to consider bidisperse systemswhich pack at a higher RCP volume fractions of φ rcp ≈ .
84 [6, 7].In parallel to studies in the field of jamming— which consider the packingproblem as a jamming transition approached from the solid phase [8, 6, 7]—other studies attempt to characterize jamming approaching the transition fromthe liquid phase [9, 10]. Here, amorphous jammed packings are seen as infinitepressure glassy states [9]. Therefore, the existence of disordered jammed struc-tures (of frictionless particles) is related to the existence of a glass transition in2d [9], a problem that has been debated recently [11].Here we treat the disordered disk packing problem with the statistical me-chanics of granular jammed matter [2]. The formal analogy with classical sta-tistical mechanics is the following: the microcanonical ensemble, defined byall microstates with fixed energy, is replaced by the ensemble of all jammed mi-crostates with fixed volume. Hence, the appropriate function for the descriptionof the system is no longer the Hamiltonian, but the volume function , W i , givingthe volume available to each particle unit such that the total system volume is V = P i W i [2].The aim of the present work is to develop the model presented in [3] for thecalculation of the volume function in the case of 2d packings. The validity ofthe hypothesis employed in [3] are discussed, and they are modified accordingto the properties of 2d packings. We use our results to study the nature of theRLP and RCP limit in 2d through an elementary construction of a statisticalmechanics that allows the study of the existence of a maximum and minimumattainable density of disordered circle packings. We find that amorphous pack-ings can pack between the density limits of ∼ ∼ . Volume function The volume function is the key of the system’s statistics: its flat averageover all jammed configurations determines the total volume. The most naturalway of dividing the system is called the Voronoi diagram, which can be seen inFig. 1a. Each grain’s region is the part of the space closer to this grain thanto any other, so that the volume is clearly additively partitioned. The majordrawback of this construction is that, so far, there was no analytical formulafor the Voronoi volume of each cell, such that attempts have been made to useother constructions [15].The Voronoi volume of a particle i can be written as: W vor i = 12 I ( min ˆ s · ˆ r ij > r ij s · ˆ r ij ) ds, (1)where ~r ij is the vector from the position of particle i to j , the average is over allthe directions ˆ s forming an angle θ ij with ~r ij as in Fig. 1a. This formula hasa simple interpretation depicted in Fig. 1a. For consistency of notation withprevious work, we will use the words ”volume” and ”surface” in 2d, althoughthey correspond to ”surface” and ”length” respectively.The volume function defined in terms of the particle coordinates is of nouse, since it does not permit the calculation of the partition function. To solvethis problem, we calculate an average free volume function based on the envi-ronment of the particle, referring to a coarse-graining over a certain mesoscopiclength scale. We assume a probability distribution for the positions of the near-est neighbors as well as the other particles. After averaging over the probabilitydistribution we obtain an average mesoscopic free volume function represent-ing quasiparticles in the partition function. Considering isotropic amorphouspackings allows for removal of the orientational averaging.
3. Probability distribution of volumes
Using the notation of Fig. 1b, we see that the microscopic volume function isentirely defined by the parameter c = min[ r/ cos θ ]. The calculation of the aver-age free volume function, w , requires knowledge of the probability distributionof this parameter. That is: w ≡ hW vor i i /V g − − Z ∞ c =1 ( c − dP > dc dc, (2)where P > ( c ) represents the inverse cumulative distribution, i.e. the probabilitythat all balls verify r ij / cos θ ij > c , which is calculated under the followinghypothesis:(1) The cumulative distribution P > ( c ) is made of two contributions, one ofthe ”contact” balls ( i.e. touching the considered grain), called P C> , and one ofthe other (background) balls, P B> . These probabilities can be understood asthe probabilities of a particle in contact (resp. background) for being situated3utside the grey zone on Fig. 1b, and therefore not contributing to the Voronoivolume c .(2) The probability distributions are those of a large number of particles at agiven density, and of negligible size giving rise to Boltzmann-like distributions: P B> ( c ) = exp( − ρV ∗ ( c ))and P C> ( c ) = exp( − ρ S S ∗ ( c ) . Here V ∗ and S ∗ represent a free volume and surface respectively towards c , i.e . V ∗ ( c ) = R Θ( c − r ˆ r · ˆ s ) d~r and S ∗ ( c ) = H Θ( c − r · ˆ s ) ds where the integrals coverrespectively the space and the unity sphere. The densities, ρ ( w ) and ρ S ( z ), aremean free-volume and free-surface densities, respectively.In 2d, the volume of the grain (with 2 R = 1) is V g = π/
4. The free volumedensity (inverse of the free volume per particle) is ρ ( w ) = N/ ( N V g φ − − N V g ) =1 / ( V g w ). Then, V ∗ ( c ) = h c − i arccos(1 /c ) + c r − c . The main assumption here is that the packing structure is uniform, thus the pairdistribution function is assumed to be a delta function at contact plus a constantfor larger distances. This assumption is an oversimplification, and more realisticbackground could be considered, such as peaks in the pair distribution at thenext nearest neighbor sites.For the surface contribution, we have: S ∗ ( c ) = 2 Z arccos(1 /c )0 dθ = 2 arccos(1 /c ) . The surface density ρ S = 1 / h S i is the inverse of the average surface left freeby z contact balls (see Fig. 2a). As a rough approximation, one can assume itis proportional to z , but because of the size of one ball, there is an ”excluded-surface” effect, so that the exact value is determined by numerical simulations.It consists in setting sequentially and randomly z non-overlapping circles ofradius 1 at the surface of the unity circle (Fig. 2a). The closest ball to theconsidered direction ˆ s defines the free angle. The free surface is then twice thisangle. Its average value is the mean free-surface h S i .Results are shown in Fig. 2b. Important deviations from the linearity in z arenot surprising, since each contact ball occupies an important surface ( z max = 6in 2d), and strong finite-surface effects are expected. For z = 5, in around 41%of the trials, the fifth ball cannot be set because there is not enough space, andwe take into account only the 59% remaining trials.For 3 ≤ z ≤
4, we will use the linear dependency ρ S ( z ) = z − . π as fitted inFig. 2. Obviously a fitting for 1 ≤ z ≤ V ∗ ( c ) and the circle occupied by the closest ”contact-grain” ( i.e. the excludedzone for the center of any other grain because of its presence) for a value of c = 1 .
2. In fact, the values of c < . φ RCP = 0 .
89, to be compared with the reportedvalue of 0 .
82. The nearer grains are exceedingly taken into account.Therefore, the volume term is modified by substituting the free volume V ∗ by V ∗ − ∆ V ∗ which represents the free volume minus the part occupied bythe closest surface grain. The meaning of this change is that the contributionsare no longer independent, and depend on two parameters c B and c C . Thedistribution P > ( c ) is the probability that both c C and c B be higher than c .Figure 3b shows the overlap of the contact grain parameterized by c C and thebackground volume parameterized by c B , defining ∆ V ∗ ( c B , c C ). The analyticalformula of ∆ V ∗ is determined by geometrical calculations.The probability density is P ( c ) = − dP > /dc = − P C> ( c ) · dP B> /dc − P B> ( c ) · dP C> /dc. The meaning of the latter equality is that, to realize c , we must have either c C or c B equal to c, and the other higher. The background probability depends on c C : P B> ( c B | c C ) = exp[ − ρ ( V ∗ ( c B ) − ∆ V ∗ ( c B , c C ))] , and P B ( c B | c C ) = − ddc B P B> ( c B | c C ) . Then, P ( c ) = Z ∞ c C = c P C ( c C ) P B ( c | c C ) dc C + P C ( c ) P B> ( c | c ) == − ddc Z ∞ c C = c P C ( c C ) P B> ( c | c C ) dc C = − ddc P > ( c ) . From (2), we integrate by parts using the latter equality. The boundaryterm [( c − P > ( c )] vanishes, since the limits of integration correspond to c = 1and c → ∞ , with P ( c → ∞ ) = 0. We obtain for the average volume function5rom Eq. (2): w = 2 Z ∞ c =1 c Z ∞ c dP C> dc C exp[ − ρ ( w ) (cid:0) V ∗ ( c ) − ∆ V ∗ ( c, c C ) (cid:1) ] dc C dc (3)with P C> ( c C ) = exp[ − ρ S ( z ) S ∗ ( c C )].
4. Free volume function
Equation (3) is a self-consistent equation to obtain w ( z ), which cannot besolved exactly, therefore a numerical integration of (3) is necessary to obtain w vs z . For various values of z , we integrate Eq. (3) numerically, and we thencalculate a fitting of the results (Fig. 4a). We obtain the free volume functionand the local density φ − i = w + 1 (Fig. 4b): w ( z ) = 0 . − . z, φ i ( z ) = 11 . − . z . (4)Generally speaking, we would expect w to be roughly proportional to 1 /z ,with w → z → ∞ . However, the statement z → ∞ has little meaningwhen we plot a figure for 3 < z < z is 6.
5. Statistical mechanics
Equation (4) plays the role of a ”Hamiltonian” of the system. Each jammedconfiguration corresponds to some ”volume level” in analogy with energy levelsin Hamiltonian systems. From the formal analogy with classical formulas, thecanonical partition function is [2]: Z ( Z, X ) = Z g ( w ) e − w/X dw, (5)where X is the (reduced) compactivity (normalized by the volume of the spheres)and g ( w ) is density of states for a given volume w . We remind that the com-pactivity is the equivalent of temperature in the Edwards statistics, and it is ameasure of the system’s looseness.Since the volume w is now directly related to coordination number z throughEq. (4), we can compute g ( w ) by replacing variable, g ( w ) = R P ( w | z ) g ( z ) dz, where P ( w | z ) is the conditional probability and g ( z ) is the density of states forgiven z [3].At this point a distinction has to be emphasized: we refer to z as the geo-metrical coordination number since it is purely defined by the particle positions.On the other hand, there is the mechanical coordination number, Z , defined bythose geometrical contacts that carry a non-zero force. Z is then defined by themechanical constraint leading to the isostatic condition. A packing is isostaticwhen the number of contact forces equals the number of force and torque bal-ance equations [3]. For example, for a packing of N infinitely rough particles6n d dimensions, each mechanical contact carries one normal force and d − d force balanceequations and d ( d −
1) torque balance equations. The isostatic condition re-quires that dN Z = dN + d ( d − N , or Z = d + 1. On the other hand, forfrictionless packings, frictional forces or tangential forces do not exist and thetorque balance equations are not taken into account. The isostatic condition inthis case leads to the relation Z = 2 d . For a system with a finite interparticlefriction µ , Z ( µ ) interpolates between both limits [7, 3]. It should be noted thatthe above calculations of the isostatic condition are based on the mechanicalcoordination number Z , rather than the geometrical coordination number z ,because a geometrical contact does not necessarily provide a mechanical con-straint. For instance, two particles are free to rotate with respect to each otherif the contact between them is geometrical and does not carry any tangentialforces.Obviously, z must be larger than Z for the mechanical condition to be satis-fied. The different volume levels of a packing can be understood in the followingway: the friction coefficient sets a mechanical constraint on Z , but the systemcan explore all geometrical levels z > Z . Additionally, z is bounded by the max-imal coordination number for a random packing, which is 2 d = 4, since there are z/ d particle coordinates. In relation with the discussion ofisostaticity, it is believed that above this value, the system is partially crystal-lized. Therefore we obtain: g ( w ) = R Z P ( w | z ) g ( z ) dz. Since w ( z ) from Eq. (4)is a coarse-grained free volume and independent of the microscopic positions ofparticles, we have P ( w | z ) = δ ( w − w ( z )).The expression of the density of states is obtained by considering that thestates are collectively jammed. Therefore, the space of configuration is discretesince we cannot continuously obtain one configuration from another. Assuminga typical distance between configurations as h z , we obtain g ( z ) ∝ ( h z ) z , theexponent z arising since there are z position constraints per particle in thejammed state compared to the free (gas) state. Such a formula is analogous tothe factor h − d for the density of states in traditional statistical mechanics, where h is the Planck constant, which arises because of the uncertainty principle, i.e. because of the discreteness of the elementary volume of phase space.Substituting into Eq. (5), we get : Z ( Z, X ) = Z Z ( h z ) z e − w ( z ) /X dz (6)To establish the maximum and minimum densities, we consider the limitsof zero and infinite compactivity, respectively. The ground state of jammedmatter, is analogous with the limit T →
0. The only accessible state is z = 4,corresponding to the random close packing. For this state, from Eq. (4) weget a fixed value of the volume fraction for any coordination number Z ∈ [3 , µ ∈ [0 , ∞ ]: φ rcp ( Z ) = 11 . − . × ≈ . , (7)7he lower density appears for X → ∞ when the Boltzmann factor is unityin Eq. (6), and we obtain the densities of RLP (assuming h z ≪ φ rlp ( Z ) = 1 Z ( Z, ∞ ) Z Z . − . z ( h z ) z dz ≈≈ . − . Z , Z ∈ [3 , . (8)This leads to the diagram in the plane ( φ, Z ) plotted in Fig. 4c defining thepossible jammed configurations. On the right part of the vertical line definedby Eq. (7), no disordered packing can exist. To the left of Eq. (8), the packingsare not mechanically stable.Between these two lines, we plot the lines of constant finite compactivity. Fora finite value of the compactivity, the equation φ ( Z ) is calculated by numericalintegration. In the figure, three curves are plotted, respectively X = 5 . − , X = 10 − and X = 10 − for h z = 0 .
01. The compactivity increases from theright ( X = 0) to the left ( X → ∞ ). The limit µ → ∞ , Z → φ minrlp = φ rlp ( Z = 3) = 11 . − . × ≈ .
775 (9)The value of φ rlp depends on the mechanical coordination number, contrar-ily to the value of φ rcp . The shape of the diagram is similar in 3d [3], and thisis in agreement with the wide range of reported values for RLP, in contrastwith RCP [4]. On a horizontal line given by a system with fixed Z , packings ofdifferent volume fractions can be achieved by applying different quench rates orcompression speeds during the preparation protocol. Slow compressions achieveloose packings (and high compactivities). The obtained predictions for the den-sity of RCP are close to the experimental values while we predict the existenceof a RLP density.
6. Summary
In summary, we have used a model of volume fluctuations to develop a sta-tistical mechanics of granular matter in 2d. From a quantitative point of view,we have seen how it lies on several approximations, that can appear too rough.The main difference with the 3d case is the need of taking into account properlythe correlations in the probability distribution of volume through the consider-ation of point (3) above. Indeed, if we take P C > and P B > to be independent asconsidered in [3] we find φ rlp = 0 .
84 and φ rcp = 0 .
89, both values above the ex-perimental value of RCP. The results, although not allowing exact predictions,are situated in the right order of magnitude for the limiting volume fractions.Due to the several approximations of the theory, the resulting limiting densi-ties have to be considered as bounds to the real values. Improvements can beachieved by taking into account the size and shape of the disks, as well as exact8numeration to calculate P > ( c ), which can be done at least to a prescribed coor-dination shell of particles, in a brute force analysis analogous to the Hales proofof the Kepler conjecture, currently being undertaken. Altogether, the presentframework seems to be successful in describing at least qualitatively the generalfeatures of jammed granular matter in 2d providing evidence of the existence ofRCP and RLP and their density value. These results suggest that a putativeideal glass transition may also exist in frictionless hard disk as discussed in [9]. References [1] A. Coniglio, A. Fiero, H. J. Herrmann, M. Nicodemi eds.
Unifying Conceptsin Granular Media and Glasses (Elsevier, Amsterdam, 2004).[2] S. F. Edwards and R. B. S. Oakeshott, Physica A , 629 (2008).[4] J. G. Berryman, Phys. Rev. A , 1053 (1983).[5] C. Brito and M. Wyart, Europhys. Lett. , 149 (2006).[6] C. S. O’Hern, S. A. Langer, A. J. Liu and S. R. Nagel, Phys. Rev. Lett. , 111 (2002).[7] T. Unger, J. Kert´esz, and D. E. Wolf, Phys. Rev. Lett. , 178001 (2005).[8] H. A. Makse, D. L. Johnson and L. M. Schwartz, Phys. Rev. Lett. , 4160(2000).[9] G. Parisi, and F. Zamponi, J. Chem. Phys. , 144501 (2005);arXiv:cond-mat/08022180[10] B. D. Lubachevsky and F. H. Stillinger, J. Stat. Phys. , 561 (1990);M. Skoge, A. Donev, F. H. Stillinger, and S. Torquato, Phys. Rev. E ,041127 (2006).[11] Y. Brumer and D. Reichman, J. Phys. Chem. B , 6831 (2004); A.Donev, F. H. Stillinger, and S. Torquato Phys. Rev. Lett. , 225502 (2006);M. Tarzia, J. Stat. Mech.: Theory and Experiment, P01010 (2007).[12] Y. Jin and H. A. Makse (2010), arXiv:1001.5287v1 [cond-mat.soft].[13] M. D. Shattuck (2006), arXiv:cond-mat/0610839v1 [cond-mat.soft].[14] D. Aristoff and C. Radin (2009), arXiv:0909.2608v1 [cond-mat.soft].[15] R. C. Ball and R. Blumenfeld, Phys. Rev. Lett. , 115505 (2002).9a)(b) Figure 1: (a) The limit of the Voronoi cell of particle i in the direction ˆ s is r ij / θ ij . Thenthe Voronoi volume is proportional to the integration of ( r ij / θ ij ) over ˆ s as in Eq. (1).(b) The particle contributing to the Voronoi volume along ˆ s is located at ( r, θ ). The dark grayregion is the considered grain ( r < R ), and in white the excluded zone for the center of anyother grain ( r < R ). For a given c = r/ cos θ , the light grey area is the region of the plane( r ′ , θ ′ ) where r ′ / cos θ ′ < c . S Numerical value Second-order fitting Linear fit for 3 Figure 2: (a) Free surface, in an example of z = 4 contact balls, defined by the angle of theclosest grain to ˆ s , with h S i = 2 θ ∗ . (b) Simulation results for 1 ≤ z ≤ 5, with a second-orderpolynomial fitting, and linear fitting for 3 ≤ z ≤ ρ S = ( z − . /π Figure 3: (a) For c = 1 . 2, the star is the center of the closest contact-grain, occupyingthe circular region printed with a pattern. The free volume is printed in grey. It isalmost completely overlapped by the surface contribution. (b) The closest contact-balldepends on c C , the free volume on c B , and ∆ V on both. Here it is the intersectionbetween the region in light grey and the patterned region. w Z w : Numerical value Linear fitting : 0.437 - 0.049 z (b) Z Volume fraction Calculated from (a) (c) Figure 4: (a) w ( z ) curve, with a linear fitting: w ( z ) = 0 . − . z . (b) Volume fractionaccording to the second Eq. (4). (c) Prediction of the model. The thick curves are the limitof the diagram at X = 0 and X → ∞ . We show several curves of constant compactivity, X .The curves are plotted for (from right to left) : X = 5 . − , X = 10 − and X = 10 − . Thehorizontal lines are both constant Z given by an arbitrary µ (dotted) and µ → ∞ (thickinferior limit of the diagram).(thickinferior limit of the diagram).