Minimal Energy Clusters of Hard Spheres with Short Range Attractions

We calculate the ground states of hard-sphere clusters, in which n identical hard spherical particles bind by isotropic short-ranged attraction. Combining graph theoretic enumeration with basic geometry, we analytically solve for clusters of n (cid:1) 10 particles satisfying minimal rigidity constraints. For n (cid:1) 9 the ground state degeneracy increases exponentially with n , but for n > 9 the degeneracy decreases due to the formation of structures with > 3 n (cid:2) 6 contacts. Interestingly, for n ¼ 10 and possibly at n ¼ 11 and n ¼ 12 , the ground states of this system are subsets of hexagonal close-packed crystals. The ground states are not icosahedra at n ¼ 12 and n ¼ 13 . We relate our results to the structure and thermodynamics of suspensions of colloidal particles with short-ranged attractions.

The structures of small clusters of atoms or particles bear directly on problems central to materials science and condensed matter physics, including nucleation, glass formation, and self-assembly [1][2][3]. The first step towards understanding the thermodynamics of a particular cluster system is to calculate the ground states as a function of particle number n [4]. Such minimal-energy clusters have been catalogued for many different potentials, the most studied being the Lennard-Jones potential [5]. But hardsphere clusters-collections of n nonoverlapping spheres with an isotropic, short-ranged, attractive pair potentialare conspicuously absent from the catalogue [6], because global optimization methods are poorly suited for strictly nonoverlapping particles.
Such clusters may yield insights into the self-assembly and nonequilibrium behavior of colloidal particles. Although the bulk phase behavior of colloidal spheres interacting through a short-ranged depletion or DNAmediated interaction has been well studied [7,8], clusters have not. The structures of energetically stable hard-sphere clusters (''inherent structures'' [9]) that are not commensurate with an equilibrium crystal may provide some clues about barriers to nucleation of colloidal crystals [10,11] or structural motifs in hard-sphere glasses [12]. Furthermore, colloids are natural building blocks for nanomaterials [13,14], and the clusters we consider are analogous to ''colloidal molecules'' [15]. The enumeration of hardsphere clusters determines the set of minimally rigid colloidal molecules that can be formed by self-assembly.
In this Letter, we consider the ground states of clusters that consist of identical spherical particles with repulsive, nonoverlapping cores and a pairwise-additive attractive potential with a range much smaller than the particle diameter [16]. The range of the attraction is short enough that the total potential energy of a cluster is linearly proportional to the number of contacts. Multiple ground states are therefore possible. Our analysis avoids optimization entirely and instead combines graph theory with geometry to analytically solve for clusters satisfying minimal rigidity constraints (!3 contacts per particle, ! 3n À 6 total contacts).
Mathematically, the clusters we enumerate correspond to minimally rigid finite sphere packings. Physically, they represent all possible colloidal molecules that can be formed from spherical particles with no barriers to bond angle rotation. Our packings include as subsets many different structures previously observed and described in the literature, such as minimal-second-moment clusters [17], colloidal clusters observed through capillary-driven assembly [13,18], and Janus particle clusters [19].
In what follows, we outline our method for enumerating packings, then discuss the results, their geometrical interpretation, and their application to colloidal systems. We focus on the ground state degeneracy, which has some striking features; the free energy of these packings at finite temperature will be analyzed in future work.
Our procedure for enumerating packings has two steps. First we use graph theory to construct all possible n-particle configurations, then we use geometry to determine which configurations correspond to minimally rigid packings. We distinguish between two types of packings: iterative packings, in which all possible m particle subsets with !3m À 6 contacts also correspond to minimally rigid packings, and noniterative packings or seeds. The majority of packings at small n are iterative. We have derived an analytical formula that quickly solves for all iterative packings; however, new geometrical rules must be derived for seeds.
The only previous effort for enumerating hard-sphere clusters that we are aware of [20] started from two seeds, a tetrahedron and an octahedron, and iteratively constructed higher order packings by combining packings of lower n with single particles. This procedure excludes any structure that combines smaller packings larger than a single particle, or that contains a different seed. As we show, both of these possibilities arise frequently as n increases.
Graph theory produces the set of possible packings.-A configuration of n particles can be described by an n Â n adjacency matrix, A. A ij ¼ 1 if the ith and jth particles touch, and A ij ¼ 0 if they do not. Given the nðn À 1Þ=2 possible contacts between n particles, there are 2 nðnÀ1Þ=2 possible A's. However, most of these are isomorphic due to particle labeling degeneracy, and thus represent the same configuration. Enumerating nonisomorphic A's constructs a complete, nonredundant set of configurations; while much smaller, the number of nonisomorphic A's still grows exponentially with n [21].
Rigidity constraints further restrict the set of possible packings. Rigidity requires (i) there be at least 3 contacts per particle, and (ii) there be at least as many contacts as internal degrees of freedom-thus there must be at least 3n À 6 contacts. Imposing these constraints eliminates all but one A for n 5 particles, though above n ¼ 5, rigidity alone does not distinguish A's corresponding to sphere packings.
Algebraic formulation.-Solving for packings requires determining the distances between spheres, if a solution exists in R 3 , and checking there are no overlaps. Each element, A ij , is associated with an interparticle distance, r ij ¼ jjz i À z j jj, which is the distance between particles whose centers are located at For A's with 3n À 6 contacts, there are precisely as many equations as unknowns. The particle configuration encoded by each A is specified by the distance matrix D, whose elements D ij ¼ r ij . If any D ij < 2r, the particles overlap and the structure is unphysical. If a continuous set of D corresponds to a given A, the structure is not rigid. The fundamental question is to find an efficient method for mapping A ! D.
Geometric solutions.-Numerical approaches for solving these equations cannot be guaranteed to converge to all D, and existing algebraic geometric approaches do not scale practically with n [22]. Instead, we use basic geometry to construct rules associating patterns of 1's and 0's in A's with either a given relative distance, D ij , or an unphysical configuration (in which case no D ! 2r exists).
First we determine if a given A is iterative by searching for subgraphs of A corresponding to lower-n seeds. The elements of D corresponding to these lower order seeds are known, and the remaining distances in D can be solved using the following property: 2 points are fixed in threedimensional space if they can be related to a common triangular base. Let there exist two particles i, j whose interparticle distance, r ij , is unknown. If there also exist 3 particles, k, p, q, with known interparticle distances (r kp , r kq , r pq ), and if the distances between i, j and the 3 particle base (r ip , r ik , r iq , r jp , r jk , r jq ) are also known; then there exists an analytical relationship for r ij [23,24] By applying this rule sequentially to all unknown distances in an A, we can solve for all relative distances of iterative packings. For a packing to be valid, the distances must satisfy the triangle inequality for each i, j, k (r ij r ik þ r kj ), all bases must lead to the same fr ij g, and all r ij ! 2r.
Any A that cannot be eliminated or solved in this fashion is a potential new seed. In this case we use geometry to analytically solve for the unknown distances. This is a reasonable proposition if there are only a few noniterative A's. At n ¼ 9 there are only 5 noniterative A's, but at n ¼ 10, there are 126. This number is still a small fraction of the total number of matrices (about 750 000 at n ¼ 10), so while overall our procedure offers a considerable computational simplification, the limiting step to the analytical enumeration of packings remains the derivation of new geometrical rules to solve for or eliminate seed candidates.
The packings.-For each packing of n 10 we have analytically solved for D [23,25]. Figure 1 shows the packing structures up to n ¼ 7. Structures for n ¼ 8, 9, 10 can be found in the supplementary information [26] and in [23]. Table I summarizes the results. The number of packings grows rapidly for n > 6. For 4 < n < 9 there is a single new seed at each n corresponding to a convex deltahedron, but for n ! 9 the number of new seeds grows rapidly. Many of these are nondeltahedral.
FIG. 1. Minimally rigid packings for n ¼ 2 À 7, with point groups indicated at upper right in Schoenflies notation. Structures at n ¼ 8, 9 can be found in the supplementary information [26], and at n ¼ 10 can be found in [23].
Each rigid packing corresponds to one D, but if a packing is chiral, D is associated with two enantiomeric structures, and thus with two states. We therefore identify the point group of each packing in order to determine chirality [23,27]. The number of states listed in Table I includes all enantiomers. The number of chiral structures also increases dramatically at n ¼ 9, where 27 out of 50 packings are chiral.
Up to n ¼ 9 every packing has exactly 3n À 6 contacts, so that all packings have the same potential energy. Figure 2 shows that for n 9 the ground state degeneracy increases exponentially. But at n ¼ 10 this trend changes due to a small number of packings that can have 25 ¼ 3n À 5 contacts (all other 10 particle packings have 3n À 6 contacts). There exist 3 such packings, each containing octahedra [ Figs. 3(b)-3(d)]. These three structures are the ground states at n ¼ 10.
Also at n ¼ 9 we find the first example of a nonrigid packing with 3n À 6 ¼ 21 contacts. Although minimally rigid, this packing has an internal rotational mode that allows it to flex without forming or breaking bonds [ Fig. 3(a)]. Adding another particle (red) onto this seed yields one of the rigid n ¼ 10 particle packings with 3n À 5 ¼ 25 contacts [ Fig. 3(b)].
Interestingly, these special packings are all subunits of the hexagonally close-packed lattice, being combinations of face-sharing tetrahedra and octahedra (Fig. 3). The structure shown in Fig. 3(d) is a subunit of both the HCP and face-centered cubic lattice. One can continue to build on to the HCP structures to create a 3n À 4 contact structure at n ¼ 11 and a 3n À 3 structure at n ¼ 12 [ Figs. 3(e) and 3(f)]. The 29-contact packing (3n À 4) at n ¼ 11 can be formed from either Fig. 3(b) or Fig. 3(c), illustrating that the family of sphere packings does not obey a strict tree structure. At n ¼ 12 the 33 contact packing (3n À 3) results from adding one more sphere to form a truncated triangular dipyramid. In all these cases the structures are commensurate with HCP.
Because we have not enumerated all the states for n ! 11, we can only conjecture that the structures we show are the ground states. What is clear is that the ground state degeneracy drops dramatically at n ¼ 10 and likely continues to oscillate with n, increasing rapidly in stretches where the maximal number of contacts, as a function of n, remains unchanged (as in n 9), and decreasing or re-  3 (color online). Structures of HCP packings. (a) Nine particle nonrigid new seed, with nonrigid motion (corresponding to a twisting of the square faces) shown as black arrows. (b)-(d) Ten particle ground state packings with 3n À 5 ¼ 25 contacts. (e) Conjectured eleven particle ground state (3n À 4 ¼ 29 contacts). (f) Conjectured twelve particle ground state (3n À 3 ¼ 33 contacts). Packings with either octahedra joined by 3 edges or half-octahedra create faces (shown in blue) that permit adding a sphere with 4 contacts (red). The joining of m octahedra by one edge (d) also yields >3n À 6 contact packings.  [20]. Chiral structures are counted as one packing. The number of ground states includes only the minimal-energy (maximal contact) structures, but does include enantiomers as separate states.
n Packings from [20] Packings (current study) New seeds Non-rigid packings  Ground states   5  1  1  0  0  1  6  2  2  1  0  2  7  4  5  1  0  6  8  1 0  1 3  1  0  1 6  9  3 2  5 0  4  1  7 7  10  113  223  8  4  3 maining small when this functional form increases. The latter can occur when an iterative n particle packing is formed by adding m spheres with >3m contacts to a minimally rigid (n À m)-particle packing. We have not yet established whether the ground states continue to be commensurate with a lattice packing for n > 12. But clearly the icosahedron is not the ground state at n ¼ 12, nor is an icosahedron with a central sphere the ground state at n ¼ 13. A twelve-sphere icosahedron has only 3n À 6 ¼ 30 contacts, and in a thirteen-sphere icosahedron the outer spheres would not be close enough to interact with each other. In equilibrium we expect finite collections of attractive colloidal spheres to form the structures we enumerate here. For n 9 the number of observable ''colloidal molecules'' with equivalent internal energy should exponentially increase, but for n ¼ 10 and, we conjecture for larger n also, there are a very small number of observable packings at low temperatures. The free energies of the packings will vary; for example, the nonrigid structure at n ¼ 9 should have a higher vibrational entropy than the other structures. We therefore expect that there will be a rich set of thermodynamic structural transitions above n ¼ 9.
Our results may connect to the nonequilibrium behavior of bulk systems. The icosahedral energy minimum for van der Waals clusters has long been argued [28] to imply icosahedral order for bulk systems; but as we have shown, icosahedra are not energy minima for hard-sphere clusters. Thus, we do not expect icosahedral order to be a structural motif in hard-sphere gels or glasses. Our results suggest that the cluster behavior of attractive hard spheres could be qualitatively different from that of many other potentials, including Lennard-Jones; it may prove fruitful to reexamine experiments [11] and simulations of bulk hardsphere systems in light of these results . Furthermore, the stability of HCP-like clusters at small n may influence the nucleation of colloidal crystals.
Finally, the methods we employ to enumerate clusters have application beyond the thermodynamics of colloidal sphere clusters: for example, they might be generalized to enumerate packings in a closed container, which would allow for explicit tests of the Edwards conjecture [29] for the entropy of a granular material.
We thank Guangnan Meng for a stimulating collaboration, and John Lee for consultations and contributions to the computer code. We also acknowledge support from the MRSEC program of the National Science Foundation under Grant No. DMR-0820484, the NSF Division of Mathematical Sciences, and DARPA under contract BAA 07-21.