Full Text

Turn on search term navigation

ARTICLE

Received 19 Feb 2013 | Accepted 25 Jun 2013 | Published 23 Jul 2013

Finding the optimal random packing of non-spherical particles is an open problem with great signicance in a broad range of scientic and engineering elds. So far, this search has been performed only empirically on a case-by-case basis, in particular, for shapes like dimers, spherocylinders and ellipsoids of revolution. Here we present a mean-eld formalism to estimate the packing density of axisymmetric non-spherical particles. We derive an analytic continuation from the sphere that provides a phase diagram predicting that, for the same coordination number, the density of monodisperse random packings follows the sequence of increasing packing fractions: spheres ooblate ellipsoids oprolate ellipsoids odimers ospherocylinders. We nd the maximal packing densities of 73.1% for spherocylinders and70.7% for dimers, in good agreement with the largest densities found in simulations. Moreover, we nd a packing density of 73.6% for lens-shaped particles, representing the densest random packing of the axisymmetric objects studied so far.

DOI: 10.1038/ncomms3194

Mean-eld theory of random close packings of axisymmetric particles

Adrian Baule1,2, Romain Mari1, Lin Bo1, Louis Portal1 & Hernn A. Makse1

1 Levich Institute and Physics Department, City College of New York, New York, New York 10031, USA. 2 School of Mathematical Sciences, Queen Mary University of London, London E1 4NS, UK. Correspondence and requests for materials should be addressed to H.A.M. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 4:2194 | DOI: 10.1038/ncomms3194 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms3194

Understanding the properties of assemblies of particles from the anisotropy of their building blocks is a central challenge in materials science13. In particular, the shape

that leads to the densest random packing has been systematically sought empirically417, as it is expected to constitute a superior glass forming material1. Despite the signicance of random packings of anisotropic particles in a range of elds like self-assembly of nanoparticles, liquid crystals, glasses and granular processing18, there is yet no theoretical framework to estimate their packing density. Thus, random packings of anisotropic particles are typically investigated on a case-by-case basis using computer simulations, which have shown, for example, that elongated shapes like prolate ellipsoids and spherocylinders can pack considerably denser than the random close packing (RCP) fraction of spheres at fRCPE0.64. These shapes exhibit a maximum in the packing fraction for aspect ratios (length/ width) close to the sphere46.

Table 1 summarizes the empirical ndings for maximal densities and highlights a further caveat of simulation and experimental studies: The protocol dependence of the nal close-packed (or jammed) state leading to a large variance of the maximal packing fractions found for the same shape. This observation can be explained using the picture of a rugged energy landscape from theories of the glass phase19. Different algorithms get stuck in different metastable basins of the energy landscape, reaching different nal packing states.

Here we present a mean-eld approach to systematically study the packing fraction of a class of anisotropic shapes with rotational symmetry, which can therefore guide further empirical studies. Explicit results are obtained for axisymmetric particles like dimers, spherocylinders and lens-shaped particles, and we discuss generalizations to other shapes like tetrahedra, cubes and irregular polyhedra. Furthermore, we derive an analytic continuation of the spherical RCP which provides a phase diagram for these and other anisotropic particles like oblate and prolate ellipsoids. We rst dene the Voronoi volume of a non-spherical particle on which our calculation is based, and show that it can be calculated analytically for many different shapes by a decomposition of the shape into overlapping and intersecting spheres, which we organize into interactions between points, lines and anti-points. We then develop a statistical mean-eld theory of the Voronoi volume to treat the particle correlations in the packing. This geometric mean-eld approach is complemented by a quantitative estimation of the variation of the average contact number with the particle aspect ratio. The predicted packing density is interpreted as an upper bound of the empirically obtained packings.

ResultsVoronoi boundary between non-spherical objects. We consider rotationally symmetric objects for which the aspect ratio a is dened as length/width, where the length is measured along the symmetry axis. In the following, we focus on the region 0oao2, where the largest densities are found17. Our description of packings relies on a suitable tessellation of space into non-overlapping volumes20. We use the standard Voronoi convention21,22, where one associates with each particle the fraction of space that is closer to this particle than to any other one. This denes the Voronoi volume Wi of a particle i, which depends on the congurations x r;^t of all particles (including

position r and orientation ^t). The total volume V occupied by N particles is V PN

Table 1 | Overview of packing fractions from simulations and experiments.

Shape /max Aspect ratio at/max

Reported z

Spherocylinder5 0.653 1.5M&M candy6 0.665 0.5 9.8 Spherocylinder14 0.689Spherocylinder8 0.694 1.4 Spherocylinder4 0.695 1.4 8.6 Dimer 0.697 1.4 8.0 Dimer12 0.703 1.4 Spherocylinder15 0.703 1.5 Spherocylinder9 0.704 1.4Oblate ellipsoid6 0.707 0.6 9.6 Dimer (theory) 0.707 1.3 8.74 Spherocylinder10 0.708 1.5 9.1 Prolate ellipsoid6 0.716 1.5 9.6 Spherocylinder17 0.722 1.5 8.7 Spherocylinder (theory) 0.731 1.3 9.5 Lens-shaped particle(theory)

0.736 0.8 9.2

General ellipsoid6 0.735General ellipsoid7 0.74 10.7 Tetrahedron13 0.76 12 Tetrahedron16 0.763Tetrahedron11 0.7858

The maximal packing fraction f and reported coordination number z at f of random packings of spherocylinders, dimers, ellipsoids and tetrahedra, determined from simulations and experiments. The aspect ratio is dened for rotationally symmetric objects. Some simulations do not report z. Results are separated by the symmetry of the object (rotationally symmetric and asymmetric) and ordered by packing fraction. From the available empirical data we cannot conclude whether spherocylinders pack better than dimers or ellipsoids of revolution, for instance.

I d^cli^c3: 2

The VB between two equal spheres is identical to the VB between two points and is a at plane perpendicular to the separation vector (Fig. 2a)20. Finding the VB for more complicated shapes is a challenging problem in computational geometry, which is typically only solved numerically23. We approach this problem analytically by considering a decomposition of the non-spherical shape into overlapping spheres. The VB is then determined as follows: Every segment of the VB arises due to the Voronoi interaction between a particular sphere on each of the two particles reducing the problem to identifying the correct spheres that interact. This identication follows an exact algorithm for a large class of shapes obtained by the union and intersection of spheres, which can be translated into an analytical expression of the VB as outlined in Fig. 3 for dimers, spherocylinders and lens-shaped particles.

For instance, a dimer is the union of a pair of spheres (Fig. 2b). The dimers VB is thus a composition of maximal four different surfaces depending on the relative orientation of the dimers dened by four points at the centre of each sphere (Fig. 3a).

i 1 Wifx1; :::; xNg, and the packing fraction

of monodisperse particles of volume Va and aspect ratio a follows as f NVa/V. In order to determine Wi one has to know the

Voronoi boundary (VB) between two particles i and j, which is the hypersurface that contains all points equidistant to both

particles (Fig. 1 for spherocylinders). The VB of the volume Wi along ^c, denoted by li^c, is the minimal one in this direction

among all possible VBs of each particle j in the packing. It is formally obtained by the global minimization20:

li^c min

j: s 4 0 srj;^tj;^c; 1 where srj;^tj;^c denotes the VB along ^c between particles i and j

with relative position rj and orientation ^tj (Fig. 1). The Voronoi volume follows then exactly as the orientational integral,

1 3

2 NATURE COMMUNICATIONS | 4:2194 | DOI: 10.1038/ncomms3194 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms3194 ARTICLE

I d^cli^c3

1 3

I d^c li^c3

1 3

I d^c Z

Figure 1 | Parametrization of the VB. The VB in blue, denoted by srj;

^tj; ^c,

along a direction ^c between two spherocylinders of relative position rj and orientation ^tj.

1 dc c3pc: 3

In the last step, we have introduced the probability density p(c), which contains the probability to nd the VB at c in the direction ^

c. The lower integration limit c ^c is the minimal value of the

boundary along ^c, which corresponds to the hard-core boundary of the particle in that direction. We introduce the cumulative distribution function (CDF) P(c) via the usual denition pc ddc Pc. Substituting the CDF in equation (3) and

performing an integration by parts leads to the volume integral

Wz Z dc Pc; z; 4 where we indicate the dependence on z. In a geometric picture20,

P(c,z) is interpreted as the probability that N 1 particles are

outside a volume O centred at c (see Fig. 4), as otherwise they would contribute a shorter VB. This leads to the denition

O^c;^t Z drYc sr;

^t;^cYsr;^t;^c; 5

where Y(x) denotes the usual Heavyside step. We refer to O as the Voronoi excluded volume, which extends the standard concept of the hard-core excluded volume Vex considered by

Onsager in his theory of elongated equilibrium rods28 (Fig. 4).

The dependence of P(c,z) on O has been treated at a mean-eld level in Song et al.20 and has been derived from a theory of correlations using liquid state theory in29 for high-dimensional sphere packings. In both cases it provides a Boltzmann-like exponential form Pc; z / exp R

s(rj,tj,c)

^ ^ c^

The extension to trimers is straightforward (Fig. 2c). Likewise, n overlapping spheres lead to compositions of n surfaces. A spherocylinder is a dense overlap of spheres of equal radii and the VB interaction is identical to that between four points and two lines (Fig. 2d). The interactions then simplify into lineline, linepoint, and pointpoint interactions, which generally lead to a curved VB for non-parallel orientations (Fig. 3b).

The Voronoi decomposition used for dimers and sphero-cylinders can be generalized to arbitrary shapes by using a dense lling of spheres with unequal radii24. However, even if it is still algorithmically well dened, this procedure may become practically tedious for dense unions of polydisperse spheres. Alternatively one can apply specialized algorithms to compute numerical VBs between curved line segments25. Here we propose an analytically tractable approach: Convex shapes can be approximated by intersections of a nite number of spheres. An oblate ellipsoid, for example, is well approximated by a lens-shaped particle, which consists of the intersection of two spheres; an intersection of four spheres is close to a tetrahedra, and six spheres can approximate a cube. This is illustrated in Fig. 2eh, and the corresponding algorithms outlined in Fig. 3c. The main insight is that the effective Voronoi interaction of these shapes is governed by a symmetry: Points map to anti-points (as the interactions between spheres is inverted; Fig. 3c). The VB of ellipsoid-like objects arises from the interaction between four anti-points and four points in two dimensions (Fig. 3c) or lines in three dimensions, and thus falls into the same class as spherocylinders. For cubes, the effective interaction is that of twelve lines, eight points and six anti-points (Fig. 2g). Analytic expressions of the VB for dimers and spherocylinders are calculated in the Supplementary Methods.

A statistical theory for Voronoi volume uctuations. We turn the above formalism into a mean-eld theory to calculate the volume fraction of a packing of monodisperse non-spherical objects. In order to take into account multi-particle correlations in the packing, we use a statistical mechanics treatment where the overall volume is expressed in terms of the average Voronoi volume Wz: V NWz (ref. 20) characterized by the average

coordination number z, which denotes the mean number of contacting neighbours in the packing. This approach is motivated by the observation that, as N-N, packings exhibit reproducible

phase behaviour, which is characterized by only few observables such as f and z (ref. 26). Our statistical mechanics framework is based on the Edwards ensemble approach, which considers the volume as a Hamiltonian of the system and attempts to nd the minimum volume27. Here W is given as the ensemble average of Wi over all particles in the packing: W Wi

h ii. We obtain

therefore from equation (2):

1 3

Ocdrrr; z

in the limit

N-N, where r(r, z) is the density of spheres at r.

The crucial step is to generalize this result to anisotropic particles. Following Onsager28, we treat particles of different orientations as belonging to different species. This is the key assumption to treat orientational correlations within a mean-eld approach. Thus, the problem for non-spherical particles can be mapped to that of polydisperse spheres for which P factorizes into the contributions of the different radii30. We thus obtain the factorized form:

Pc; z exp Z d

8 >

Oc;^t

drrr;^t; z

9 >

;

; 6

where rr;^t; z is the density of particles with orientation ^t at r.

Next, we assume an approximation of this density in terms of contact and bulk contributions, which is motivated by the connection with the radial distribution function in spherical theories in both high and low dimensions20,29. The contact contribution relies on the condition of contact between the two particles of a given relative position r and orientation ^t, which denes the contact radius r ^r;^t: r is the value of r for which the

two particles are in contact without overlap. In the case of equal spheres the contact radius is simply r ^r;^t 2a. For non-

spherical objects, r ^r;^t depends on the object shape and the

relative orientation (Supplementary Methods). Using r ^r;^t we

NATURE COMMUNICATIONS | 4:2194 | DOI: 10.1038/ncomms3194 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms3194

Object shape Decomposition Effective Voronoi interaction

Sphere

Dimer

Trimer

Spherocylinder

Ellipsoid

Tetrahedron

Cube Six spheres Twelve lines, eight points,

six anti-points

Irregular polyhedron Unequal spheres Points, lines, anti-points

Figure 2 | Decomposition of various shapes and effective Voronoi interactions. Arbitrary object shapes can be decomposed into unions and intersections of spheres. (ad) Union of spheres. The VB between two such objects is equivalent to the VB between the point multiplets at the centre of the spheres, as shown for four basic shapes. (eg) Intersection of spheres. The VB between such intersections is equivalent to that between multiplets of anti-points at the centre of the spheres, indicated by crosses, and, in addition, lines at the edges of the intersections, shown as points in ed. The additional lines arise due to the positive curvature at the singular intersections, resulting in edges that point outwards from the particle rather than inwards. In the case of dimers and trimers shown in b,c, the curvature is negative and the edges do not inuence the VB. The generalization to (f) tetrahedra-like, (g) cubes and (h) irregular polyhedra-like shapes is straightforward. Note that the VBs drawn in eh are only qualitative.

4 NATURE COMMUNICATIONS | 4:2194 | DOI: 10.1038/ncomms3194 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

One sphere Two points

Two spheres Four points

Three spheres Six points

Two lines and four points

N spheres

Two spheres

Two lines and four anti-points

Four spheres Six lines, four points, four anti-points

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms3194 ARTICLE

Construct separation lines from the spherical decomposition

Separate different interactions between two objects

Dimer

Ellipsoid (lens-shaped particle)

x^ x^

Spherocylinder

Figure 4 | The Voronoi excluded volume and surface. (a) The hard-core repulsion between the two objects denes the hard-core excluded volume Vex (enclosed by a dashed blue line): This volume is excluded for the centre of mass of any other object. Packings of rods in the limit a-N can be described by a simple random contact equation based on Vex (ref. 51). We introduce the Voronoi excluded volume O (enclosed by a dashed red line), which is the basis of our statistical theory of the Voronoi volume. The volume O, equation (5), is excluded by the condition that no other particle should contribute a VB smaller than c in the direction ^c, which denes the

CDF P(c, z). (b) Taking into account the hard-core exclusion leads to the effective Voronoi excluded volume V* (indicated as red volume), which is excluded for bulk particles. Likewise, the overlap of Vex and O excludes the surface S* (thick green line) for all contacting particles. The volumes are shown here for a single orientation ^t. (c) The 3D plot corresponding to b:

The central particle is in brown, Vex is indicated in blue, V* in red, and S* in green.

Figure 3 | Analytical solution to determine the VB for non-spherical objects. (a) The VB between the two objects of a given relative position and orientation consists of the VBs between particular spheres on each of the two objects. The spheres that interact are determined by separation lines given as the VBs between the spheres in the lling. For dimers, there is one separation line for each object, tesselating space into four areas, in which only one interaction is correct. The pink part in a, for example, is the VB between the two upper spheres. (b) The dense overlap of spheres in spherocylinders leads to a line as effective Voronoi interaction at the centre of the cylindrical part. This line interaction has to be separated from the point interactions due to the centres of the spherical caps as indicated. Overall, the two separation lines for each object lead to a tessellation of space into nine different areas, where only one of the possible lineline, linepoint, pointline, and pointpoint interactions is possible. The yellow part in b, for example, is due to the upper point on spherocylinder 1 and the line of 2. Regions of line interactions are indicated by blue shades. (c) The spherical decomposition of ellipsoid-like shapes is analogous to dimers, only that now the opposite sphere centres interact. We indicate this inverted interaction by a cross at the centres of the spheres and refer to these points as anti-points. In addition, the positive curvature at the intersection point leads to an additional line interaction, which is a circle in 3D (a point in 2D) and indicated here by two points. The separation lines are then given by radial vectors through the intersection point/line. The Voronoi interaction between the two ellipsoids is thus given by two pairs of two anti-points and a line, which is the same class of interactions as spherocylinders. The different point and line interactions are separated analogous to spherocylinders, as shown.

Here we have explicitly written the dependence of r on W, which is important to interpret equation (4) as a self-consistent equation to obtain the volume fraction of the packing. The free-volume per particle in the bulk depends specically on Wz as

r 1=Wz Va:

The CDF thus factorizes into two contributions: A contact term:

PCc; z exp szS c

f g; 9

and a bulk term:

PBc exp rWV c

;

can separate bulk and contact terms in rr;^t; z as in20,29:

rr;^t; z

14p rYr r r;

^t szdr r ^r;^t :

The prefactor 1/4p is the density of orientations, which we assume isotropic. The symbols r and s(z) stand for the average

free-volume of particles in the bulk and the average free-surface of particles at contact, respectively, which are discussed further below. The approximation equation (7) corresponds to considering a pair distribution function as a delta function modelling the contact particles plus a constant term modelling the particles in the bulk29, which are thus considered as a uniform structure. These assumptions are further tested in the Methods section.

Substituting equation (7) into equation (6) leads to our nal result for the CDF:

Pc; z exp rWV c szS c

NATURE COMMUNICATIONS | 4:2194 | DOI: 10.1038/ncomms3194 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms3194

such that

Pc; z PCc; z PBc: 11 The volume V* is the volume excluded by O for bulk particles and takes into account the overlap between O and the hard-core excluded volume Vex: V O O\Vex

h i^t, where :::

h i^t denotes an

orientational average. Likewise, S* is the surface excluded by O for contacting particles: S @Vex\O

h i^t, where @Vex denotes the

boundary of Vex. The volumes Vex and O as well as the resulting V* and S* are calculated in the Supplementary Methods and shown in Fig. 4 for spherocylinders.

The surface density s(z) is a measure for the available surface for contacts when the packing is characterized by an average coordination number z. We evaluate this density by simulating random local congurations of one particle with z non-overlapping contacting particles and determining the average available free-surface. This surface is given by S*(cm), where cm is the minimal contributed VB among the z contacts in the direction ^c. Averaging over many realizations with a uniform distribution of orientations and averaging also over all directions ^c provides the surface density in the form,

1 S cm

h i

h i^c

: 12

In this way, we can only calculate s(z) for integer values of z. For fractional z that are predicted from our evaluation of degenerate congurations in the next section, we use a linear interpolation to obtain Wz.

Equations (4) and (8) lead to a self-consistent equation for the average Voronoi volume Wz in the form: Wz FWz .

Analytic expressions for V* and S* can be derived in the spherical limit in closed form, where also the self-consistency equation can be solved exactly20. For non-spherical shapes we resort to a numerical integration to obtain V* and S*. Equation (4) can then be solved numerically, which yields Wz, and subsequently the

equation of state for the volume fraction versus coordination number, fz; a Va=Wz, in numerical form (denoting

explicitly the dependence on a).

Variation of the coordination number with aspect ratio. In this purely geometric theory of the average Voronoi volume, the packing fraction is given as f(z, a), with z and a free parameters, in principle. In practice, z is xed by the symmetry properties of

the object shape, z(a), and the physical condition of mechanical stability, requiring force and torque balance on every particle. Under the assumption of minimal correlations, these conditions typically motivate the isostatic conjecture based on Maxwells counting argument31: z 2df, with df the number of degrees of

freedom, giving z 6 for fully symmetric objects (spheres), z 10

for rotationally symmetric shapes like spherocylinders, dimers and ellipsoids of revolution6, and z 12 for shapes with three different

axis like aspherical ellipsoids and tetrahedra13. While the isostatic conjecture is well-satised for spheres, packings of non-spherical objects are in general hypoconstrained with zo2df, where z(a)

increases smoothly from the spherical value for a41 (ref. 6). The fact that these packings are still in a mechanically stable state can be understood in terms of the occurrence of stable degenerate congurations (Fig. 5), which reduce the effective number of degrees of freedom32. However, the observed variation z(a) could not be explained quantitatively so far. Here we deduce the relation z(a) by evaluating the probability of nding these degenerate congurations to provide a prediction of f(a) in close form.

In a degenerate conguration, force balance already implies torque balance, as the net forces are aligned with the inner axis of the particle (Fig. 5). This implies that there is redundancy in the set of force and torque balance equations for mechanical equilibrium as force and torque balance equations are not linearly independent. Our evaluation of these degenerate congurations is based on the assumption that a particle is always found in an orientation such that the redundancy in the mechanical equilibrium conditions is maximal. This condition allows us to associate the number of linearly independent equations involved in mechanical equilibrium with the set of contact directions. Averaging over the possible sets of contact directions then yields the average effective number of degrees of freedom ~dfa, from which the coordination number follows as

za 2

~dfa (Methods).

The results for z(a) are shown in Fig. 6a for prolate ellipsoids of revolution, spherocylinders, dimers and lens-shaped particles. We are able to recover the observed continuous transition as a function of a from the isostatic coordination number for spheres, z 6 at a 1, to the isostatic value z 10, for aspect ratios above

E1.5. The trend compares well to known data for ellipsoids6 and spherocylinders10,17. In particular, our approach explains the decrease of z for higher aspect ratios observed in simulations of spherocylinders10,17: For large a, the most probable case is to have

f1 + f2

f3 +f4

Figure 5 | Quantitative method to calculate z(a). (a) A 2D sketch of a spherocylinder with a random conguration of contact directions rj. The associated forces are along directions ^

nj normal to the surface (indicated in red) and torques are along rj ^

nj. From these directions, one can determine if mechanical equilibrium has some redundancy, that is, if force and torque balance equations are not linearly independent. The conguration shown has no redundancy: The equivalent situation in three dimensions would show force and torque balance equations as ve different constraints (the most general case for a 3D particle would be six constraints, but the torque along the axis of a spherocylinder is always vanishing, due to its rotational symmetry). (b) Here the spherocylinder is rotated. With the same contact directions rj as in a, the contact force directions ^

nj are now modied. We explore the space of possible orientations for the spherocylinder, and try to nd congurations which maximize redundancy in the mechanical equilibrium conditions. (c) As an example, this orientation exhibits some redundancy: All the contacts are on the spherical caps of the spherocylinder. Therefore, f1 f2 and f3 f4 are aligned with the

spherocylinder axis and the condition of force balance automatically implies torque balance. If this is the orientation of the spherocylinder for which redundancy is maximal, we associate the number of linearly independent equations (that is, the effective number of degrees of freedom) from the mechanical equilibrium condition with the set of contact directions {rj} and perform an average over the possible sets of {rj}. This yields the averaged effective number of degrees of freedom ~dfa for a spherocylinder having an aspect ratio a and the coordination number follows as za 2

~dfa. Note that

for non-convex shapes like dimers, the resulting z is the number of contacting neighbours, not the number of contacts, which can exceed the former.

6 NATURE COMMUNICATIONS | 4:2194 | DOI: 10.1038/ncomms3194 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms3194 ARTICLE

contacts only on the cylindrical part of the particle, so that all normal forces are coplanar reducing the effective number of degrees of freedom by one. Consequently, z-8 as a-N, as we obtain in Fig. 6a. This decrease is specic to spherocylinders, and not observed for dimers or ellipsoids, since the normal forces are not coplanar.

Phase diagram of non-spherical particles. Our calculation leads to a close theoretical prediction for the packing density f(a)

f(z(a), a) which does not contain any adjustable parameters. Figure 6b shows the prediction for dimers, spherocylinders and lens-shaped particles. For spherocylinders, results in the literature on f(a) vary greatly (Table 1), but all show a peak at around aE1.31.5, which is captured by our formalism. We predict the maximum density of spherocylinders at a 1.3 with a density

fmax

0.731 and that of dimers at a 1.3 with fmax 0.707. We

have also calculated the packing fraction of the lens-shaped particles of Fig. 3c, which yields fmax 0.736 for a 0.8. This

shape represents the densest random packing of an axisymmetric shape known so far.

We further investigate packings of non-spherical objects in the zf representation. This change in perspective allows us to characterize packings of differently shaped objects in a phase diagram. By plotting z(a) against f(a) parametrically as a function of a, we obtain a phase diagram for jammed anisotropic

particles in the zf plane (Fig. 6c). In the same diagram, we also plot the equation of state obtained with the present theory in the case of spheres in20: fsphz z=z 2 3

p , which is valid

between the two isostatic limits of frictionless spheres z 6 and

innite frictional spheres at z 4. Surprisingly, we nd that both

dimer and spherocylinder packings follow an analytical continuation of these spherical packings. This result highlights that the spherical random branch can be continued smoothly beyond the RCP in the zf plane.

The analytical continuation of RCP is derived by solving the self-consistent equation (4) close to the spherical limit (Supplementary Methods):

fz 1 o1

1 g1o1 z z 1

h i

Mb Mz

g2o1 z z 1

Mb Mz

h i

Mv Mz

here o1 1= 3

p denotes the spherical free-volume at RCP dened as o1 1=fsph 1 evaluated at z 6 as calculated in20,

z 6 is the spherical isostatic value, and the functions g1,2 can be

expressed in terms of exponential integrals. The dependence of equation (13) on the object shape is entirely contained in the geometrical parameters Mb, Mv, and Mz: Mb and Mv quantify the rst-order deviation from the sphere at a 1 of the objects

0.74

FCC

0.72

0.70

Z([afii9825] )

0.68

[p10]([afii9825])

0.66

Lens-shaped particles theory

Ellipsoids theory Dimers theory Spherocylinders theory

0.64

RCP

Lens-shaped particles theory

0.62

Dimers theory Spherocylinders theory

0.5 1.0 1.5 2.0 2.5

0.60 0.5 1.0 1.5 2.0

[afii9825]

[p10]

[afii9825]

Spherocylinder M&M candy Spherocylinder Spherocylinder Spherocylinder

Dimer

Dimer Spherocylinder Spherocylinder Oblate ellipsoid Spherocylinder Prolate ellipsoid

Spherocylinder

Abreu et al.5 Donev et al.6

Lu et al.14

Jia and Williams8 Williams and Philipse4

Schreck and OHern, personal communication

Faure et al.12

Kyrylyuk et al.15

Bargiel et al.9

Donev et al.6 Wouterse et al.10

Donev et al.6

Zhao et al.17

RCP

Z([p10])

0.62 0.66 0.70 0.74

RCP

Spherical random branch

Dimers theory Spherocylinders theory Analytic continuation Eq.(13)

4 0.55 0.60 0.65 0.70 0.75

Figure 6 | Theoretical predictions for packings of particles with various shapes. (a) The function z(a) determined by evaluating the probability of degenerate congurations. Both spherocylinders and dimers increase up to just below the isostatic value z 10. For dimers, z(a) is the number of

contacting neighbours, not the number of contacts, as a single contacting particle can have more than one contacting point. For spherocylinders, z reduces to 8 for large a, as the forces acting on the cylindrical part are coplanar and reduce the effective degree of freedom. We also include the results from our method for prolate ellipsoids of revolution and lens-shaped particles. (b) The predicted packing fraction f(a) of spherocylinders, dimers and lens-shaped particles compared with simulation results of maximal densities from the literature. We predict the maximal packing fraction of spherocylinders fmax

0.731 at a 1.3 and of dimers fmax 0.707 at a 1.3, demonstrating that spherocylinders pack better than dimers. For the lens-shaped particles we

obtain fmax 0.736 at a 0.8. (c) By plotting z vs f we obtain a phase diagram for smooth shapes. We observe that the spherical random branch fsph,

which ends at the RCP point at (0.634,6) (ref. 20), in fact continues smoothly upon deformation into dimers and spherocylinders as predicted by our theory. The spherocylinder continuation provides a boundary for all known packing states of rotationally symmetric shapes. Inset: The continuations from RCP. For a given value of z, the densest packing is achieved by spherocylinders, followed by dimers, prolate ellipsoids and oblate ellipsoids. Note that the continuations for spherocylinders and dimers are almost identical.

NATURE COMMUNICATIONS | 4:2194 | DOI: 10.1038/ncomms3194 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms3194

hard-core boundary and its volume, respectively, while Mz measures the rst-order change in the coordination number upon deformation of the sphere. The resulting continuations z(f)

obtained by inverting equation (13) for different object shapes are

plotted in the inset of Fig. 6c.

For the smooth shapes considered, we nd generally that denser packing states are reached for higher coordination numbers. For a given value of z, spherocylinders achieve the densest packing, followed by dimers, prolate ellipsoids and oblate ellipsoids, as seen in the inset of Fig. 6c. We observe that the densest packing states for dimers and spherocylinders found in simulations lie almost exactly on the continuation, while the one of the ellipsoids deviate considerably.

Comparison with empirical data. Table 1 indicates that there is a nite range of densities for random jammed packings according to the particular experimental or numerical protocol used (denoted as a J-line in the case of jammed spheres19,33). On the other hand, our mean-eld theory predicts a single density value and Fig. 6 indicates that our predictions are an upper bound of the empirical results. We interpret these results in terms of current views of the jamming problem developed in the limiting case of spheres, where the question of protocol-dependency of packings has been systematically investigated.

Random close packings can be considered as innite-pressure limits of metastable glass states, which was shown theoretically in refs 19,3436 and conrmed in computer simulations in ref. 37. Indeed, there exist a range of packing fractions named as [fth, fGCP] following the notation of mean-eld Replica Theory (RT)19. Here fGCP stands for the density of the ideal glass close packing and is the maximum density of disordered packings, while fth is the innite-pressure limit of the least dense metastable states. In RT, the states [fth, fGCP] are all isostatic.

From the point of view of simulations, the well-known LubachevskyStillinger (LS) protocol33 provides this range of packings for different compression rates. The densities [fth,

fGCP] are achieved by the corresponding compression rates (from large to small) [gth, gGCP-0]. Compression rates larger than gth all end to fth. The threshold value gth corresponds to the relaxation time 1/gth of the least dense metastable glass states. The denser states at GCP are unreachable by experimental or numerically generated packings, as it requires to equilibrate the system in the ideal glass phase, a region where the relaxation time is innite. In general, large compression rates lead to lower packing fractions. This picture was investigated for sphere packings in33,38 and it is particularly valid for high-dimensional systems where crystallization is avoided19.

Random close packings are also known to display sharp structural changes3943 signalling the onset of crystallization at a freezing point fc (ref. 18). All the (maximally random) jammed states along the segment [fth, fGCP] can be made denser at the

cost of introducing some partial crystalline order. Support for a order/disorder transition at fc is also obtained from the increase of polytetrahedral substructures up to RCP and its consequent decrease upon crystallization44. In terms of protocol preparation like the LS algorithm, there exists a typical time scale tc corresponding to crystallization. Crystallization appears in LS18,19,41 if the compression rate is smaller than gc 1/tc,

around the freezing packing fraction45. A possible path to avoid crystallization and obtain RCP in the segment [fth, fGCP] is to

equilibrate with g4gc to pass the freezing point, and eventually setting the compression rate in the range [gth, gGCP-0] to

achieve higher volume fraction.

As the present statistical mechanics framework is based on the Edwards ensemble approach27, our prediction of the packing

density fEdw corresponds to the ensemble average over the conguration space of random states at a xed coordination number. As the volume has the role of the Hamiltonian, the energy minimization in equilibrium statistical mechanics is replaced in our formalism by a volume minimization: The highest volume fraction for a given disordered system is achieved in the limit of zero compactivity. Therefore, the present framework provides a mean-eld estimation of such a maximal volume fraction (minimum volume) of random packings with no crystallization. As we perform an ensemble average over all packings at a x coordination number, the obtained volume fraction fEdw corresponds to the one with the largest entropy (called largest complexity in RT) along [fth, fGCP]. This point

needs not to be fth, and in general it is a larger volume fraction. Thus, fth o fEdw o fGCP.

The above discussion can be translated to the present case of non-spherical particles. In this case, unfortunately, there is no detailed study of the protocol-dependent packing density as done by19,33,38 for spheres. However, the survey of the available simulated data obtained by different groups (Table 1 and Fig. 6b,c) can be interpreted analogously as for spheres. In the case of spherocylinders, packings have been obtained in the range[0.653, 0.722] (these minimum and maximum values have been obtained in ref. 5 and in ref. 17, respectively, see Table 1). Our predicted density is 0.731, representing an upper bound to the simulated results. In the case of dimers, there are two simulations giving a density of 0.697 (C. F. Schreck and C. S. OHern (2011) personal communication) and 0.703 (ref. 12), which are both smaller than and very close to our prediction 0.707. Thus, our prediction is interpreted as the upper limit in the range of packings observed with numerical algorithms. Under this scenario, which is consistent with analogous three-dimensional (3D) spherical results, packings may exist in the region [fth, fEdw], and our theory is a

mean-eld estimation of fEdw. This region is very small for spheres but the above evidence indicates that non-spherical particles may pack randomly in a broader range of volumes. The present framework estimates the upper bound for such a range.

DiscussionWe would like to stress that our analytic continuation is non-rigorous and appears as the solution of our mean-eld theory for rst-order deviations in a from the sphere using suitable approximations. The shapes of dimers, spherocylinders, ellipsoids are then all shown to increase the density of the random packing to rst-order. In the case of regular (crystal) packings, recent mathematically rigorous work has shown in fact that for axisymmetric particles any small deformation from the sphere will lead to an increase in the optimal packing fraction of the crystal46. This appears only in 3D and is related to Ulams conjecture stating that the sphere is the worst case scenario for ordered packings in 3D (ref. 47). A full mathematical proof of this conjecture is still outstanding, but so far all computer simulations verify the conjecture. In particular, recent advances in simulation techniques allow to generate crystal packings of a large variety of convex and non-convex objects in an efcient manner48,49. The extensive study of de Graaf et al.48 has extended the verication of Ulams conjecture to the rst eight regular prisms and antiprisms, the 92 Johnson solids, and the 13 Catalan solids. The verication for regular n-prisms and n-antiprisms can be extended to arbitrary n using this method, providing an exhaustive empirical verication of the conjecture for these regular shapes. We remark that a random analogue of Ulams packing conjecture has been proposed and veried for the Platonic solids (apart from the cube) in simulations16. The results presented here support the random version of Ulams conjecture

8 NATURE COMMUNICATIONS | 4:2194 | DOI: 10.1038/ncomms3194 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms3194 ARTICLE

and might help in investigating this conjecture further from a theoretical point of view.

We believe that our decomposition of various shapes into intersections and overlaps of spheres will be a useful starting point for a systematic investigation of this issue. Our approach can be systematically continued beyond the axisymmetric shapes considered here. For instance, in Fig. 2eh, we have 2,3,6,n anti-points to describe ellipsoids and polyhedra of increasingly varying complexity. The challenge would be to implement our algorithm to calculate the resulting Voronoi excluded volumes that appear in our mean-eld theory. For this, one might also consider a fully numerical evaluation using, for example, graphics hardware25.

Methods

Quantitative method to calculate z(a). Mathematically, we can write the local mechanical equilibrium on a generic non-spherical frictionless particle having k contacts dened by their location rj, normal ^

nj and force fj^

nj, as:

^ is a minimum. Note that N^t depends on ^ti, as only the absolute

direction of contact points are chosen, and thus rotating particle i affects the direction and normal of these contacts with respect to particle i. This situation is described in Fig. 5c, which includes a two-dimensional (2D) sketch of a 3D degenerate conguration that we observe often in our procedure. In this case, the rank is reduced by one unit, and the probability of occurrence of such a situation is large at small aspect ratio, as it just requires that there is no contact on the cylindrical part of the inner particle.

Within our assumptions, we explore the space of possible contact directions for one particle, given a local contact number k, and aspect ratio a. We then extract the average effective number of degrees of freedom ~df a; k, which is

the average over the contact directions of the minimal value of rankN^t:

~df a; k h min^trankN

^ i fr ; ... ;r g, where :::

h i r ; ;r

f g N 1 R

J :::dr1:::drk

denotes the average over contact directions. This average is limited to a subset J of all possible {r1yrk} such that mechanical equilibrium (equation 14) is possible with positive forces, as expected for a packing of hard particles. This corresponds geometrically to sets {r1,y,rk} which do not leave a hemisphere free on the unit sphere. Finally, the normalization N is the volume of J. For a packing with a

coordination number distribution Qz(k), with average z, the effective df is:

~df a P

k Qzk

~df a; k, and the average z follows as za 2

~df a. In our

evaluation, we use a Gaussian distribution for Qz(k), with variance 1.2 and average z, consistent with simulations50. Overall, z(a) is thus the solution of the following self-consistent relation:

za 2 X

f1 ... fk

n1 ::: ^

r1 ^

^ n1 ::: rk ^

N f 0; 14

where N is a df k matrix. A local degenerate conguration has a matrix N

such that rankN o mindf ; k. We base our evaluation on two assumptions: (i)

Contact directions around a particle in the packing are uncorrelated, and (ii) Given one set of contact directions, a particle i is found in an orientation ^ti such that the

redundancy in the mechanical equilibrium conditions is maximal, that is, rankN

: 15

The way we look for the orientation ^t on the unit sphere showing the lowest rank is simply by sampling it randomly with a uniform distribution (106 samples).

Qzk min

^ fr ; ... ;r g

^ t rankN

[afii9825]=1.1, [afii9835]c=0.22 [afii9825]=1.1, [afii9835]c=0.8 [afii9825]=1.1, [afii9835]c=1.51

[afii9825]=1.5, [afii9835]c=0.22 [afii9825]=1.5, [afii9835]c=0.8 [afii9825]=1.5, [afii9835]c=1.51

[afii9825]=2.0, [afii9835]c=0.22 [afii9825]=2.0, [afii9835]c=0.8 [afii9825]=2.0, [afii9835]c=1.51

100

101

102

103

100

101

102

103

100

101

102

103

100

101

102

103

100

101

102

103

100

101

102

103

100

101

102

103

100

101

102

103

100

101

102

103

1.2 1.4 1.6 1.8 2.0 1.2 1.4 1.6 1.8 2.0 1.0 1.2 1.4 1.6 1.8 2.0

1.6 1.8 2.0 2.2 2.4 1.4 1.6 1.8 2.0 2.2 2.4 1.0 1.2 1.4 1.6 1.8 2.0

2.0 2.2 2.4 2.6 2.8 3.0 3.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 1.0 1.2 1.4 1.6 1.8 2.0

c/a c/a c/a

Figure 7 | Comparison of the CDF with simulation data. We plot the theoretical predictions (solid lines) for P(c, z) (black), PB(c) (red) and PC(c, z) (green) with the corresponding CDFs sampled from simulated congurations (symbols) of spherocylinders. For each aspect ratio a 1.1, 1.5, 2.0, we plot results for

three values of the polar angle ycA[0,p/2]. We generally observe that the three CDFs agree quite well in the regime of small c values, which provides the dominant contribution to the average Voronoi volume Wz. The same plots are shown on a linear scale in the Supplementary Fig. S1. The error bars denote

the root mean square error of the nite-size sampling.

NATURE COMMUNICATIONS | 4:2194 | DOI: 10.1038/ncomms3194 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms3194

The computation of the rank is done via a standard Singular Value Decomposition of N^t, which is here numerically accurate for aZ1.05.

Test of the approximations of the theory. We perform a comprehensive test of the different approximations of the theory using computer simulations of spherocylinder packings (Supplementary Note 1). From the generated congurations at the jamming point we obtain the CDF P(c, z), where z is also an observable of the simulation determined by the jamming condition. P(c, z) contains the probability that the boundary of the Voronoi volume in the direction ^c is found at a value larger than c and is determined as follows. We select an orientation ^c relative to the orientation ^ti of a chosen reference particle i. A large number of particles in the packing contribute a VB along ^c with particle i. We determine all these different VBs denoted by srj;^tj;^c. The boundary of the Voronoi volume in the direction ^c

is the minimum cm of all positive VBs:

cm min

j: s 4 0 srj;^tj;^c; 16 where rj and ^tj are the relative position and orientation of particle j with respect to the reference particle i. Determining this minimal VB for all particles i in the packing yields a list of cm values for a given ^c (which is always relative to the orientation ^ti). The CDF P(c,z) simply follows by counting the number of values larger than a specied c.

Owing to the rotational symmetry of the spherocylinders, the orientational dependence of P(c, z) is reduced to P(c, yc; z), where yc is the polar angle of the orientation ^c in spherical coordinates. Moreover, due to inversion symmetry it is sufcient to select only yc 2 0; p=2 . Therefore, we choose three yc values to cover

this range: yc 0.22, 0.8, 1.51. We also use the rotational symmetry to improve the

sampling of P(c, z): We x yc to one of the three values, but select a number of azimuthal angles at random. As the packing is statistically isotropic for all azimuthal angles, the resulting cm value for these directions can all be included in the same ensemble. We consider three different aspect ratios a 1.1, 1.5, 2.0

of the spherocylinders to capture a range of different shapes. The results are plotted in Fig. 7.

We test the two main approximations considered in the theory: (a) The derivation of P(c, z) using a liquid like theory of correlations as done in Songet al.20 Jin et al.29 leading to the exponential form of equation (8). (b) The factorization of this CDF into contact and bulk contributions as in equation (11). This approximation neglects the correlations between the contacting particles and the bulk. In Fig. 7, we test these approximations by comparing theory and simulations for three different CDFs: P(c, z), PB(c) and PC(c, z), equations (8)(10).

In order to determine the PB(c) from the simulation data we need to take the contact radius r ^rj;^tj between particle i and any particle j into account. The

minimal VB, cm, is determined from the contributed VBs of particles in the bulk only, that is, particles with rj 4 r ^rj;^tj. Likewise, PC(c) is determined from the

simulation data by only considering VBs of contacting particles with rj r ^rj;^tj.

Following this procedure, we have tested these approximations with the computer generated packings. We nd (Fig. 7): (i) The contact term PC is well approximated by the theory for the full range of c; (ii) For small values of c the bulk distribution PB is well approximated by the theory, and deviations are observed for larger c; (iii) The full CDF P(c) agrees well between the computer simulations and the theory, especially for small c. The small values of c provide the dominant contribution in the self-consistent equation to calculate the average Voronoi volume equation (4), and therefore to the main quantity of interest, the volume fraction of the packing. This can be seen by rewriting equation (4) as

Wz Va I d^cZ

1c ^c

References

1. Glotzer, S. C. & Solomon, M. Anisotropy of building blocks and their assembly into complex structures. Nat. Mater. 6, 557562 (2007).

2. Damasceno, P. F., Engel, M. & Glotzer, S. C. Predictive self-assembly of polyhedra into complex structures. Science 337, 453457 (2012).

3. Ni, R., Gantapara, A. P., de Graaf, J., van Roij, R. & Dijkstra, M. Phase diagram of colloidal hard superballs: from cubes via spheres to octahedra. Soft Matter 8, 88268834 (2012).

4. Williams, S. & Philipse, A. Random packings of spheres and spherocylinders simulated by mechanical contraction. Phys. Rev. E 67, 051301 (2003).

5. Abreu, C., Tavares, F. & Castier, M. Inuence of particle shape on the packing and on the segregation of spherocylinders via Monte Carlo simulations. Powder Technol. 134, 167180 (2003).

6. Donev, A. et al. Improving the density of jammed disordered packings using ellipsoids. Science 303, 990993 (2004).

7. Man, W. et al. Experiments on random packings of ellipsoids. Phys. Rev. Lett. 94, 198001 (2005).

8. Jia, X. M. G. & Williams, R. A. Validation of a digital packing algorithm in predicting powder packing densities. Powder Technol. 174, 1013 (2007).9. Bargiel, M. Geometrical properties of simulated packings of spherocylinders. Computational ScienceICCS2008 5102, 126135 (2008).

10. Wouterse, A., Luding, S. & Philipse, A. P. On contact numbers in random rod packings. Granular Matter 11, 169177 (2009).

11. Haji-Akbari, A. et al. Disordered, quasicrystalline and crystalline phases of densely packed tetrahedra. Nature 462, 773777 (2009).

12. Faure, S., Lefebvre-Lepot, A. & Semin, B. in Ismail, M., Maury, B. & Gerbeau, J.-F. (eds) ESAIM: Proceedings Vol. 28, 1332, (2009).

13. Jaoshvili, A., Esakia, A., Porrati, M. & Chaikin, P. M. Experiments on the random packing of tetrahedral dice. Phys. Rev. Lett. 104, 185501 (2010).14. Lu, P., Li, S., Zhao, J. & Meng, L. A computational investigation on random packings of sphere-spherocylinder mixtures. Science China 53, 22842292 (2010).

15. Kyrylyuk, A. V., van de Haar, M. A., Rossi, L., Wouterse, A. & Philipse, A. P. Isochoric ideality in jammed random packings of non-spherical granular matter. Soft Matter. 7, 16711674 (2011).

16. Jiao, Y. & Torquato, S. Maximally random jammed packings of platonic solids: hyperuniform long-range correlations and isostaticity. Phys. Rev. E 84, 041309 (2011).

17. Zhao, J., Li, S., Zou, R. & Yu, A. Dense random packings of spherocylinders. Soft Matter 8, 10031009 (2012).

18. Torquato, S. & Stillinger, F. H. Jammed hard-particle packings: From Kepler to Bernal and beyond. Rev. Mod. Phys. 82, 26332672 (2010).

19. Parisi, G. & Zamponi, F. Mean-eld theory of hard sphere glasses and jamming. Rev. Mod. Phys. 82, 789845 (2010).

20. Song, C., Wang, P. & Makse, H. A. A phase diagram for jammed matter. Nature 453, 629632 (2008).

21. Aurenhammer, F. Voronoi diagrams - a survey of a fundamental geometric data structure. ACM Computing Surveys 23, 345405 (1991).

22. Okabe, A., Boots, B., Sugihara, K. & Nok Chiu, S. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams (Wiley-Blackwell, 2000).

23. Boissonat, J. D., Wormser, C. & Yvinec, M. in Effective Computational Geometry for Curves and Surfaces (eds Boissonnat, J. D. & Teillaud, M.) 67 (Springer, 2006).

24. Phillips, C. L., Anderson, J. A., Huber, G. & Glotzer, S. C. Optimal lling of shapes. Phys. Rev. Lett. 108, 198304 (2012).

25. Hoff, K., Culver, T., Keyser, J., Lin, M. & Manocha, D. Fast computation of generalized voronoi diagrams using graphics hardware. in SIGGRAPH 99 Conference Proceedings. (Computer Graphics) 277286 (ACM SIGGRAPH Assoc Computing Machinery, 1999).

26. Makse, H. A., Bruji, J. & Edwards, S. F. in The Physics of Granular Media (eds Hinrichsen, H. & Wolf, D. E.) (Wiley-VCH, 2004).

27. Edwards, S. F. & Oakeshott, R. B. S. Theory of powders. Physica A 157, 10801090 (1989).

28. Onsager, L. The effects of shape on the interaction of colloidal particles. Ann. N.Y. Acad. Sci. 51, 627659 (1949).29. Jin, Y., Charbonneau, P., Meyer, S., Song, C. & Zamponi, F. Application of Edwards statistical mechanics to high-dimensional jammed sphere packings. Phys. Rev. E 82, 051126 (2010).

30. Danisch, M., Jin, Y. & Makse, H. A. Model of random packings of different size balls. Phys. Rev. E 81, 051303 (2010).

31. Alexander, S. Amorphous solids: their structure, lattice dynamics and elasticity. Phys. Rep. 296, 65236 (1998).

dc Pc;^c; z; 17

as the CDF is trivially unity for c values smaller than the hard-core boundary c ^c.

The main contribution to the integral then comes from c values close to c ^c due

to the decay of the CDF.

Systematic deviations in our approximations arise in the bulk distribution PB for larger values of c, but, interestingly, the slope of the decay still agrees with our theory. Overall, the comparison highlights the mean-eld character of our theory: Correlations are captured well up to about the rst coordination shell of particles, after which theory and simulations diverge, especially for the bulk term. The agreement is acceptable for the nearest neighbour-shell, but is incorrect for the second neighbours. Beyond this shell, bulk particles are affected in a nite range by correlations that we do not address, as we assume a uniform distribution of the density of these particles; this is a typical assumption in a mean-eld theory. The additional unaccounted correlations lead to a slightly higher probability to observe the VB at intermediate c values in the simulation, compared with our theory. However, these deviations from simulations are small. For instance, Fig. 7 indicates that for a typical value a 1.5 and polar angle yS 0.22, the numerically measured

CDF P(c,z) at a relative large value c/a 2 is of the order of 10 3, while the theory

predicts this probability at a slightly larger value of c/a 2.07. This small

discrepancy is not relevant, as such a value of the probability is negligibly small in the calculation of the volume fraction in equation (4). Thus, because of this small probability to nd the VB with values larger than c/a 2, the deviations expected

from our approximations are small. These results indicate that, overall, the theory captures the distribution of VBs in the region of small c, which is the relevant region in the calculation of the volume fraction.

The neglected higher-order correlations in the upper coordination shells can only decrease the volume fraction in the calculation leading to smaller packing densities. Following this analysis, we interpret our predicted packing fractions as upper bounds for the empirically found ones, which is indeed observed in Fig. 6b,c.

10 NATURE COMMUNICATIONS | 4:2194 | DOI: 10.1038/ncomms3194 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms3194 ARTICLE

32. Donev, A., Connelly, R., Stillinger, F. H. & Torquato, S. Underconstrained jammed packings of nonspherical hard particles: Ellipses and ellipsoids. Phys. Rev. E 75, 051304 (2007).

33. Skoge, M., Donev, A., Stillinger, F. H. & Torquato, S. Packing hyperspheres in high-dimensional euclidean spaces. Phys. Rev. E 74, 041127 (2006).34. Krzakala, F. & Kurchan, J. Landscape analysis of constraint satisfaction problems. Phys. Rev. E 76, 021122 (2007).

35. Mari, R., Krzakala, F. & Kurchan, J. Jamming versus glass transitions. Phys. Rev. Lett. 103, 025701 (2009).

36. Biazzo, I., Caltagirone, F., Parisi, G. & Zamponi, F. Theory of amorphous packings of binary mixtures of hard spheres. Phys. Rev. Lett. 102, 195701 (2009).

37. Hermes, M. & Dijkstra, M. Jamming of polydisperse hard spheres: The effect of kinetic arrest. Europhys. Lett. 89, 38005 (2010).

38. Chaudhuri, P., Berthier, L. & Sastry, S. Jamming transitions in amorphous packings of frictionless spheres occur over a continuous range of volume fractions. Phys. Rev. Lett. 104, 165701 (2010).

39. Anikeenko, A. V. & Medvedev, N. N. Polytetrahedral nature of the dense disordered packings of hard spheres. Phys. Rev. Lett. 98, 235504 (2007).

40. Radin, C. Random close packing of granular matter. J. Stat. Phys. 131, 567573 (2008).

41. Jin, Y. & Makse, H. A. A rst-order phase transition denes the random close packing of hard spheres. Physica A 389, 53625379 (2010).

42. Klumov, B. A., Khrapak, S. A. & Morll, G. E. Structural properties of dense hard sphere packings. Phys. Rev. B 83, 184105 (2011).

43. Kapfer, S. C., Mickel, W., Mecke, K. & Schrder-Turk, G. E. Jammed spheres: Minkowski tensors reveal onset of local crystallinity. Phys. Rev. E 85, 030301 (2012).

44. Anikeenko, A. V., Medvedev, N. N. & Aste, T. Structural and entropic insights into the nature of the random-close-packing limit. Phys. Rev. E 77, 031101 (2008).

45. Cavagna, A. Supercooled liquids for pedestrians. Phys. Rep. 476, 51124 (2009).46. Kallus, Y. & Nazarov, F. In which dimensions is the ball relatively worst packing? Preprint at http://arxiv.org/abs/1212.2551

Web End =http://arxiv.org/abs/1212.2551 (2012).

47. Gardner, M. The Colossal Book of Mathematics: Classic Puzzles, Paradoxes, and Problems (Norton, 2001).

48. de Graaf, J., van Roij, R. & Dijkstra, M. Dense regular packings of irregular nonconvex particles. Phys. Rev. Lett. 107, 155501 (2011).

49. de Graaf, J., Filion, L., Marechal, M., van Roij, R. & Dijkstra, M. Crystal-structure prediction via the Floppy-Box Monte Carlo algorithm: Method and application to hard (non)convex particles. J. Chem. Phys. 137, 214101 (2012).

50. Wang, P., Song, C., Jin, Y. & Makse, H. A. Jamming II: Edwards statistical mechanics of random packings of hard spheres. Physica A 390, 427455 (2011).

51. Philipse, A. The random contact equation and its implications for (colloidal) rods in packings, suspensions, and anisotropic powders. Langmuir 12, 11271133 (1996).

Acknowledgements

We gratefully acknowledge funding by NSF-CMMT and DOE Ofce of Basic Energy Sciences, Chemical Sciences, Geosciences and Biosciences Division. We are grateful to C.F. Schreck and C.S. OHern for discussions and for providing simulated data on 3D packings of dimers. We are also grateful to F. Potiguar for discussions, T. Zhu for simulations and M. Danisch for theory. We also thank F. Zamponi, P. Charbonneau andY. Jin for discussions on the interpretation of protocol-dependent packings.

Author contributions

A.B., R.M., L.B., L.P. and H.A.M. designed research, performed research and wrote the paper.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Baule, A. et al. Mean-eld theory of random close packings of axisymmetric particles. Nat. Commun. 4:2194 doi: 10.1038/ncomms3194 (2013).

NATURE COMMUNICATIONS | 4:2194 | DOI: 10.1038/ncomms3194 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 11

Word count: 9884

Show less

Abstract

Translate

Finding the optimal random packing of non-spherical particles is an open problem with great significance in a broad range of scientific and engineering fields. So far, this search has been performed only empirically on a case-by-case basis, in particular, for shapes like dimers, spherocylinders and ellipsoids of revolution. Here we present a mean-field formalism to estimate the packing density of axisymmetric non-spherical particles. We derive an analytic continuation from the sphere that provides a phase diagram predicting that, for the same coordination number, the density of monodisperse random packings follows the sequence of increasing packing fractions: spheres <oblate ellipsoids <prolate ellipsoids <dimers <spherocylinders. We find the maximal packing densities of 73.1% for spherocylinders and 70.7% for dimers, in good agreement with the largest densities found in simulations. Moreover, we find a packing density of 73.6% for lens-shaped particles, representing the densest random packing of the axisymmetric objects studied so far.

Details

Title

Mean-field theory of random close packings of axisymmetric particles

Author

Baule, Adrian; Mari, Romain; Bo, Lin; Portal, Louis; Makse, Hernán A

Pages

2194

Publication year

2013

Publication date

Jul 2013

Publisher

Nature Publishing Group

e-ISSN

20411723

Source type

Scholarly Journal

Language of publication

English

DOI

https://doi.org/10.1038/ncomms3194

ProQuest document ID

1411095210

Mean-field theory of random close packings of axisymmetric particles

Jump to:

Full Text

Abstract

Details

Suggested sources