Spontaneous Patterns in Disk Packings


Boris D. Lubachevsky
Ron L. Graham
Frank H. Stillinger
Lucent Technologies  University of California Lucent Technologies 
Bell Labs Innovations 
at San Diego 
Bell Labs Innovations
700 Mountain Avenue 
La Jolla, CA 92093 
700 Mountain Avenue 
Murray Hill, NJ 07974 
     Murray Hill, NJ 07974 
bdl@bell-labs.com graham@ucsd.edu fhs@bell-labs.com



Using a computational procedure that imitates tightening of an assembly of billiard balls, we have generated a number of packings of n equal and non-equal disks in regions of various shapes for n in the range 10 to 10,000. Our experiments revealed various geometric patterns and regular features, including polycrystalline textures with "rattlers'' typically trapped along the grain boundaries. We also observed series of packings with similar patterns for n taking on increasing sequences of values, n = n(1), n(2),... n(k),....

Billiards algorithm

Packings of non-overlapping disks or cylinders can be mechanically generated by tightening the container boundary [7]. A variant of this procedure has been implemented on a computer [4] as a discrete event "billiards'' simulation algorithm [3]. Rather than shrinking the container, in the billiards simulations we uniformly expand the disks. Fig. 1 illustrates the work of the billiards simulation algorithm while packing 2000 equal disks in a square with periodic boundary (torus). The initial stage at time t = 0 is represented on the top square in the figure where 2000 points are randomly scattered. To each point an initial velocity vector is randomly assigned (not shown). Some points lie outside the square; they are periodic images of the corresponding points inside. When t > 0, the points grow into disks, and all disks at each time t have common diameter d = Et. The growth continues until the configuration "jams,'' at which time we expect to have a packing.

At t = 0 disks do not overlap because their sizes are zero. For t > 0 disks move along straight lines with given velocities; their motions may conflict with each other. At an instance of such a two-disk conflict, we simulate an elastic collision of these disks assuming their masses are equal. At a collision both disks change their velocity vectors so as to exchange momenta and energies according to known mechanical laws.Note that the evolution of the disk configuration does not change if we proportionally change both the disk expansion speed E and linear disk velocities.We normalize the simulation input by assuming the unity mean initial disk velocity.

Figure 1: Successive stages in an instance of disk packing by billiards algorithm.

Figure 2: A packing of 2000 disks in a square with periodic boundary. The packing is obtained under a fast disk  expansion, E = 100.The packing consists of crystalline grains with many rattlers represented as unshaded disks concentrated along the grain boundaries. Monovacancies occur within the hexagonally packed grains.

A number of packings have been produced by randomly varying initial configuration of points, their initial velocity, and the expansion speed E. Packings generated under a high E (in Fig. 2), a moderate E (in Fig. 3), and a low E (in Fig. 4) were obtained. For high E the pattern is a combination of hexagonally packed fragments, "grains,''with each grain having generally different orientation. Irregularity concentrates along the grain boundaries, as do the "rattlers,'' the disks that are not rigidly fixed in their positions but trapped in the cages formed by their fixed rigid neighbors or boundary walls if any.

Figure 3: A packing of 2000 disks in a square with periodic boundary. The packing is obtained under a moderately fast disk expansion, E = 3.2, and consists of grains that are larger than those in Fig. 2.As in Fig. 2, the rattlers concentrate along the grain boundaries but their number and that of the monovacancies is smaller than in Fig.2.

The size of the grains increases (and their number for a fixed n decreases) as the E decreases. For a sufficiently small E, a single hexagonal "crystal'' emerges. The array of hexagonally packed disks contains interesting structural deviations and insertions.Sometimes the symmetry of the packing coupled with a high value of its density suggests that perhaps we have achieved the optimum, i.e., the packing of the largest possible disk diameter and density. This might be the case in Fig. 4. (The density is the fraction of the region area covered by disks.)

Figure 4: A packing of 2000 disks in a square with periodic boundary obtained under a slow disk expansion, E = 10-3. If the monovacancy near the center is filled with the 2001-st disk, the obtained packing seemingly  becomes perfectly symmetric. Might that be the optimum packing of 2001 equal disks in a square with periodic boundary? Its experimentally computed density (when the 2001-st disk is inserted) is 0.901635...

Frustrated packings
56 equal disks stacked hexagonally in 8 alternating columns, 7 disks in a column, constitute the optimum packing in a rectangle with periodic boundary conditions if the ratio of the height of the rectangle to its width is (1 + 1/48)1/2 = 1.01036... We frustrate the packing by increasing the width of the rectangle and making it equal to the height. In the new shape, previously jammed disks become loose and can grow further to exploit the extra space. As the disks grow, a new jamming is achieved with a larger disk diameter.

Figure 5: The best found packing of 56 disks in a square with periodic boundary conditions. (The square perimeter is dashed). Pairs of disks that are in contact are indicated by line segments  connecting their centers. Triangles formed by such lines are shaded. The shading helps to see the structure of the packing: the plane is tiled with 4×7 blocks, each block having pattern " 'Z," shaded on it.

In the example of a frustration in Fig. 6 we begin with a perfectly hexagonally packed configuration of 10,800 equal disks. Note that we can always perfectly pack 3k2 equal disks in a regular hexagon with periodic boundary conditions, as explained in the figure. k disks will fit along each hexagon side. Here k = 60.

We uniformly decrease the size of all the disks but one in the center. The diameter of the latter is increased. As a result all the disks become loose. We tighten the configuration by simulating expansion with very small rate E. During the expansion, while diameters of all disks are increasing, we keep constant the ratio of the diameter of the large disk to that of other disks.

Figure 6: Displacement vectors of centers of 10800 disks packed in a regular hexagon. Boundary conditions are periodic as shown in the right top corner box where the same letter, a or b, marks identical points. If all disks were of the same size ("pure'' crystal), the packing would be perfectly hexagonal. Displacements from this perfect order are caused by the central "impurity'' disk being 20% larger than the others. Displacement vector of the central disk is 0 by definition, each other displacement vector begins at the perfect order position and its length is magnified 20 fold to enhance the resolution.

As a result, a frustration with a "crack'' centered at the larger disk is developed. To enhance visually the structure of the frustration, we display disk displacements. Each disk has the "ideal'' position of its center p0 in the original non-frustrated packing and the final position of the center p1 in the frustrated packing.Vector p1 - p0 is the displacement. Accumulation of the round off error may cause the entire configuration to shift during the expansion. We "calibrate'' the displacement vectors by subtracting from each of them the displacement of the larger disk. The latter displacement becomes zero by definition.

Repeated patterns in packings in equilateral triangle, 
square, and circle with hard walls

The optimum packing of D(k) = k(k+1)/2 equal disks in an equilateral triangle is the hexagonal arrangement [6] and the optimality holds for all k = 1, 2, .... Are the triangle numbers D(k) the only such lucky sequence? We conjecture the existence of an infinite number of sequences so that each sequence has its own well-defined pattern of optimal packings.

Figure 7: The conjectured densest packing of np(k) = 256 disks inside an equilateral triangle, where p = 5 and k = 3 and np(k) is defined by formula (1). The densest packings of n =np(k) disks for all checked values  of p and k, have this pattern consisting of one triangle of side (k+1)p-1 and 2p+1 alternating triangles of side k  with p-1 rattlers that are "falling off'' the larger triangle.

For each p = 0, 1, 2,... consider sequence 

np(k) =  D((k+1)p-1) + (2p+1)D(k),   k = 1, 2, ...               (1)

For p = 0 sequence np(k) is identical with the sequence of triangle numbers D(k) with known optimal packings. For each p > 0, we conjecture the optimal packing of np(k) disks as the pattern that consists of n-p+1 solid disks and p-1 rattlers; it includes one triangle of side (k+1)p-1 and 2p+1 alternating triangles of side k each as shown in Fig. 7 for the case p = 5 and k = 3.

A k×k orthogonal arrangement is a "natural'' square packing in the same spirit as hexagonal disk arrangement is a natural packing in an equilateral triangle. This square arrangement becomes non-optimal for large k, though. What we observe in the best found packings in the square is an interplay between two patterns: square and hexagonal. For a sufficiently large number of disks n the hexagonal pattern becomes dominant.With the billiards simulation algorithm we were able to examine details of this interplay for small n.

The following sequences were recently identified [1, 5] as candidates for pattern repetitions: k2-3, k2 -2,k2-1. These are "frustrations'' of the "natural'' square pattern. Some members of these sequences are presented in Fig. 8. In this and the following figures black dots identify contacts (between two disks in a pair or between a disk and the boundary). Sometimes, disks are seemingly in contact, like disks 2 and 15 in the packing of 22 disks in Fig. 8, but no dot is present in the contact site. This indicates a positive gap between the disks. The gap may be not discernible on the picture, but can be verified computationally. Because a pair of rows and a pair of columns are "frustrated'' in best packings of k2-2 disks (middle row in Fig. 8), the two pairs of indices can identify the packing, like (3,5; 2,4).

Sequences k(k+1) and k2+ [k/2] were also identified (where [x] denotes largest integer that is not exceeding x); they can be considered as "frustrations'' of a hexagonal pattern adjusted to the square boundaries. Some of their members are presented in Fig. 9. Note that the best packing of 10 disks does not follow the common pattern of sequence k2+ [k/2].

As in the case of a square, a circle boundary shape does not conform to the hexagonal disk arrangement. Unlike the case of a square or equilateral triangle, no obviously "natural'' packing of equal disks has been proposed. Perhaps, the curved hexagonal packings [2] can be taken as such. They are presented in the left column in Fig. 10. Three different equally dense curved hexagonal packings of 61 disks are known. Fig. 10 presents one of them. Similarly, 12 different equally dense curved hexagonal packings of 91 disks are known and Fig. 10 presents one of them.The curved hexagonal packings of 3k(k+1)+1 disks exist for arbitrary large k but they become not the densest ones for k = 6, 7,... (n = 127, 169, ...). The right column in Fig. 10 presents best found packings for the sequence 3k(k+3)+1. These demonstrate a different way of adjusting hexagonal packing to the circular-wall boundary.

Figure 8: Best found packings of equal disks in a square. Top row: first (k = 5) and last (k = 8) members of sequence k2-3. The packings consist of a heavier shaded (k-3)×(k-3) square packing in the bottom left corner and three lighter shaded alternating rows and columns and one unshaded rattler. Middle row: two out of four existing best packings of 34 disks, a member of sequence k2-2 for k = 6. Each packing consists of a (k-2)×(k-2) heavier shaded square pattern with two lighter shaded hexagonally arranged rows and two columns. Two pairs of insertion indices identifies a packing. Bottom row: two out of three existing best packing of 35 disks, a member of sequence k2-1 for k = 16. Each packing consists of a (k-1)×(k-1) heavier shaded square pattern with one  lighter shaded inserted row andone column. A pair of insertion indices identifies a packing.

Figure 9: Members of sequences k(k+1) (56 disks) and k2+[k/2] (5, 10, 18, and 52 disks). Optimal packing of 10 disks (middle row, left) does not follow the common pattern of the latter sequence. The inferior packing of 10 disks that follows the pattern is also shown (top row, right)

Figure 10: Left: best found packings of 3k(k+1)+1 equal disks in a circle for k = 1, 2, 3, 4, and 5. Right: best found packings of 3k(k+3)+1 equal disks in a circle for k = 2, 3, and 4 and the highest obtained density configuration of 121 disks (k = 5).



[1] R. L. Graham and B. D. Lubachevsky, Repeated patterns of dense packings of equal disks in a square, The Electronic Journ.of Combinatorics 3 (1996), #R16.
[2] R. L. Graham, B. D. Lubachevsky, K. J. Nurmela, and P.R. J. Östergård, Dense packings of congruent circles in a circle, Discrete Mathematics, 181 (1998), 139-154.
[3] B. D. Lubachevsky, How to simulate billiards and similar systems, J. Computational Physics 94 (1991), 255-283.
[4] B. D. Lubachevsky and F. H. Stillinger, Geometric properties of random disk packings, J. Statistical Physics60 (1990), 561-583.
[5] K. J. Nurmela and P. R. J. Östergård, Packing up to 50 equal circles in a square, Discrete & Computational Geometry,18 (1997), 111-120.
[6] N. Oler, A finite packing problem, Canad. Math. Bull.4 (1961), 153-155.
[7] G. E. Reis, Dense packings of equal circles within a circle, Math. Mag. 48 (1975), 33-37.


VisMath HOME