Dense packings of spheres in cylinders I. Simulations
Abstract
We study the optimal packing of hard spheres in an infinitely long cylinder, using simulated annealing, and compare our results with the analogous problem of packing disks on the unrolled surface of a cylinder. The densest structures are described and tabulated in detail up to (ratio of cylinder and sphere diameters). This extends previous computations into the range of structures which include internal spheres that are not in contact with the cylinder.
pacs:
I Introduction
The dense packing of monodisperse (equalsized) hard spheres in a cylinder has been found to produce a remarkable sequence of interesting structures as the ratio of the cylinder diameter to the sphere diameter is varied. We have explored computationally the densest of these packings in detail using simulated annealing, following the work of Pickett et al Pickett:2000 (), and the first four packings in this sequence are shown in Fig 1. We have reported extensive results, together with analytic approximations which serve to elucidate our findings Mughal:2011 (). The present paper and its future continuations amplify the previous reports by Mughal et al. Mughal:2011 () and by Chan (second author of the present paper) chan:2011 (), and extend them in various directions. In particular, new results are provided for larger cylinder diameters and we include full details of all structures in a form suitable for future reference.
We shall refer to these quasi 1D packings as columnar crystals since they are periodic. Each structure can be assembled by stacking unit cells ad infinitum along the length of the cylinder with each subsequent unit cell rotated by the same twist angle with respect to the previous one.
For smaller cylinder diameters, up to and beyond the range investigated by Pickett et al. Pickett:2000 (), the densest packings only consist of spheres that are in contact with the surface of the cylinder. Eventually at larger cylinder diameters dense packings are found for which there are internal spheres that are not in contact with the surface. Nevertheless we extend the study of structures with only cylindertouching spheres to larger diameters, for its own sake. As already explained in Mughal:2011 (), these structures are readily understood by recourse to a yet simpler problem, in which circular disks are placed on a cylinder. This can be fully worked out in analytical terms.
In parallel with our computational study using simulated annealing and a separate investigation by Chan chan:2011 () using sequential deposition, a number of experiments have been undertaken on the packing of solid spheres and small bubbles under gravity. Many of the simulated structures are observed. While we will have cause to mention these experiments, their presentation will be reserved for a further paper Meagher:2012 ().
Some of the structures that we describe (particularly the more elementary ones) have been observed in biological microstructures Erickson:1973 (). Columnar crystals have also been realised in numerous experimental contexts including foams (both dry Pittet:1996 (); Weaire:1992 (); Pittet:1995 (); Boltenhagen:1998 (); Hutzler:2002 () and wet Meagher:2012 ()) in tubes, colloids in micro channels Moon:2004 (); Moon:2005 (); Li:2005 (); Tymczenko:2008 (); Lohr:2010 (), and fullerenes in nanotubes Khlobystov:2004 (); Yamazaki:2008 (); Warner:2010 ().
We have not attempted to rigorously prove that the densest packings identified by simulations are indeed the densest possible. This ought to be achievable for some of the simplest ones.
Ii Motivation and Background
The identification of dense packings has always played an important part in condensed matter physics and physical chemistry Aste:2008 (). The entities which are to be packed are often spherical, or well approximated as such, and they may be monodisperse. In the idealised case of hard spheres (with no interaction when not in contact) the problem of finding the maximum density with given boundary conditions constitutes a classic mathematical challenge. If posed for an unbounded packing, it is long associated with the name of Kepler. It is still debated, in terms of full mathematical rigour, although there is no room for practical doubt about the nature of densest packings (fcc etc) Hales ().
Here we are concerned with the case in which spheres of diameter are to be contained in an unbounded cylinder of diameter . We shall pursue it up to which represents the limit of the capability of our present computational resources and methods.
Our primary objective is to identify the structure with the largest value of , defined as the fraction of volume occupied by the spheres.
There are close connections between this topic and that of phyllotaxis, a subject arising out of biology and having to do with the dense arrangement of similar units on the surface of a cylinder, exemplified by pine cones, pineapples or corn cobs Airy:1872 (). This connection will be explored in detail here.
Iii Some simple columnar crystals
For certain specific values of the optimal packing structure can be easily surmised. These structures are directly related to the analogous twodimensional problem of finding the smallest diameter circle into which nonoverlapping circles, each of diameter , can be packed. A description of these special cases will serve to illustrate the general problem that we will address.
The first five of these special solutions, labelled , are shown schematically in Fig 2 and (for ) consist of disks placed at the vertices of a regular polygon. The diameter of the enclosing circle is given by,
(1) 
where , , , .
Replacing the circles with spheres we define a unit cell which can be repeated along the length of the tube as follows. Each successive layer, indicated by alternating red and green spheres, can be generated by translating the previous layer through a distance
(2) 
along the tube and rotating it by a twist angle . Using the above formula the volume fraction of these simple columnar crystals can easily be computed, , , , , .
We shall label these simple columnar crystals with the notation to indicate that their structure can be derived in a simple way from the circle packing problem.
Iv Simulation
For general values of , the optimal packing structure cannot be guessed easily and we must turn to heuristic methods. Our search method is confined to structures that are periodic in the following sense. There is a primitive cell, of length , containing spheres, the structure being generated from this by the screw operation of (i) translation along the cylinder axis by (where is any integer) combined with (ii) rotation about the axis by an angle . This screw operation represents the underlying symmetry of columnar crystals (which are not to be confused with columnar phases, originating in the study of liquid crystals and related instead to the packing of columns).
Our primary method of simulation is based on the well known approach of simulated annealing. This provides a reasonably exhaustive and unbiased search for maximum density. Appendix A summarises the technical details of the present application (such as annealing schedules).The search procedure described in Appendix A looks for the minimum possible value of , for a given , treating and the sphere positions as variables. It does this for increasing values of until the (likely) optimal structure is apparent, or further computation is not practical.
In the previous simulations of Pickett et al. Pickett:2000 (), a periodic boundary condition without any twist was used. If for example the densest structure has a twist angle equal to , then ten cells would combine to satisfy the simpler boundary condition.
V Numerical Results
In this section we present the full range of results for several significant properties as a function of . These are the volume fraction and chirality, for the densest structures that we have found. The following figures and Table 1 summarise the main results.
Table 1 enumerates all the densest packings observed, from our numerical work, in the range . The table lists the range of over which each structure is observed, the number of spheres in the primitive unit cell (from which the extended structure can be constructed) and the average number of contacts per sphere. The table also classifies the structures into two groups: those which are chiral and those which are not (achiral)  see below for a more detailed discussion of chirality.
The break in Table 1, between structures 32 and 33, highlights the fact that beyond the densest structures are those which contain spheres not in contact with the cylindrical boundary. The unit cell for packings with internal spheres are described using the notation , where the first integer is the number of internal spheres and the second is the number of spheres which are in contact with the cylindrical boundary. So for example in the case of structure 30 the unit cell contains 6 spheres of which 5 are in contact with the cylindrical boundary and 1 is not.
Structure  Range  Number of spheres in Unit Cell  Notation  Average Contact Number  Description  chirality 

1 ()  1  straight chain  2    achiral  
2  1  zigzag  2    achiral  
3  1  twisted zigzag  4    chiral  
4  2  (2,1,1)  4  lineslip  chiral  
5 ()  2  (2,2,0)  5  maximal contact  achiral  
6  2  (2,2,0)  5  lineslip  chiral  
7  1  (3,2,1)  6  maximal contact  chiral  
8  2  (3,2,1)  5  lineslip  chiral  
9  3  (3,2,1)  16/3  lineslip  chiral  
10 ()  3  (3,3,0)  6  maximal contact  achiral  
11  3  (3,3,0)  16/3  lineslip  chiral  
12  2  (3,2,1)  5  lineslip  chiral  
13  2  (4,2,2)  6  maximal contact  achiral  
14  2  (4,2,2)\(4,2,2)  5  lineslip  chiral  
15  3  (3,3,0)  16/3  lineslip  chiral  
16  1  (4,3,1)  6  maximal contact  chiral  
17  3  (4,3,1)  16/3  lineslip  chiral  
18  4  (4,3,1)  11/2  lineslip  chiral  
19 ()  4  (4,4,0)  6  maximal contact  achiral  
20  4  (4,4,0)  11/2  lineslip  chiral  
21  3  (4,3,1)  16/3  line slip  chiral  
22  1  (5,3,2)  6  maximal contact  chiral  
23  3  (5,3,2)  16/3  lineslip  chiral  
24  4  (4,4,0)  11/2  lineslip  chiral  
25  1  (5,4,1)  6  maximal contact  chiral  
26  4  (5,4,1)  11/2  lineslip  chiral  
27  5  (5,4,1)  28/5  lineslip  chiral  
28 ()  5  (5,5,0)  6  maximal contact  achiral  
29  5  (5,5,0)  28/5  lineslip  chiral  
30  6 (1,5)  26/6  (1,5)  chiral  
31  6 (1,5)  40/6  maximal contact  chiral  
32  6 (1,5)  26/6  +(1,5)  chiral  
33  11 (1,10)  40/11  (1,10)  see Sect. VIII  
34  11 (1,10)  60/11  maximal contact  achiral  
35  11 (1,10)  58/11  +(1,10)  chiral  
36  7 (2,5)  32/7  (2,5)  chiral  
37  7 (2,5)  34/7  maximal contact  chiral  
38  7 (2,5)  30/7  +(2,5)  chiral  
39  15 (2,13)  72/15  (2,13)  chiral  
40  15 (2,13)  90/15  maximal contact  chiral  

v.1 Volume fraction
Fig 3 presents the maximum density found by our prescribed procedure, for the full range of that we have explored. By simple finite difference we approximate the derivative of the volume fraction as a function of , as shown in Fig 4. This is to clearly identify the singular behaviour discussed below.
The vertical lines indicate a discontinuity in the derivative of the volume fraction; there is either (i) only a sudden changes in the derivative, or (ii) the simultaneous presence of such a sudden change and a squareroot singularity in the derivative, the two cases of which are denoted by vertical bluedashed and reddotted lines, respectively. Note the squareroot singularities in the derivative coincide with the (circle packing) structures.
Let us discuss this behaviour by reference to Fig (5) which is a simplified cartoon of the observed changes in the volume fraction as a function of . Qualitatively, the variation of the maximum density takes the following form: its dependence on is everywhere continuous while the derivative is piecewise continuous. We distinguish two types of singular points, as in Fig (3), at which the derivative changes.
The first occur when a maximum number of contacts is reached. These points correspond to highly symmetric structures, such as the structures described above. At such points  represented by the vertices labelled and , in the inset  the previous trend of structural change cannot be continued; in effect, the structure has jammed due to the formation of new contacts. In other words: if the separations of contacting sphere centres is to be maintained, the system is overconstrained and no such solution exists beyond this point. Instead, some existing contacts are released and a new structural trend proceeds. The structure itself varies continuously through the singular point, with a downward change in the derivative of the density, corresponding to the line segments and in the cartoon.
At the second kind of singular point the structure is simply overtaken by another more dense packing. Here the structure itself changes discontinuously, and the derivative obviously must change in a positive sense. In the inset this transition corresponds to the points labelled and where, for example, the line segment illustrates the increasing trend in density immediately following such a transition. The dashed lines indicate the continuation of the previous structure which now has a lower density compared with the optimal packing.
These remarks are based on observation of the results. We have for example no proof that the change of derivative at points of the first type is always negative, although we can offer a proof that the dependence of (maximum) density on cannot have any discontinuities at which the density undergoes a finite upward or downward change. This is discussed in detail in appendix C.
v.2 Chirality
A chiral object is one which cannot be superimposed on its mirror image (or inverse)by translations and (proper) rotations. The question of chirality is interesting in the present context, for example in relation to designing chiral molecular filters that will discriminate between enantiomers.
According to this definition an object is either chiral or it is not (it is achiral). The reported structures are classified accordingly in Table 1.
But chirality may be manifested in physical properties to a greater or lesser extent: the degree of chirality of the structure itself is a tempting concept. While it may be useful to think of this in practice, there is no unique quantitative definition that will have general relevance. Any appropriate property could be used as an index of chirality. Nevertheless one may offer a working definition for use in relation to simulations, as Pickett et al did Pickett:2000 (). Here we will employ another definition that is perhaps more transparent.
We define a chirality index in terms of the degree to which an object can be superimposed upon its mirror image. For a given columnar packing, our method is to start with the sphere centres and generate a mirror image of the structure by reflecting the coordinates of the centres in the  plane, i.e. the crosssectional plane of the cylinder. For each sphere in the packing we compute its distance to the nearest sphere in the mirror image. The overlap function is defined as the sum of these distances divided by the number of spheres in the original packing. Clearly when a packing can be completely superimposed upon its mirror image the overlap function vanishes.
The computational challenge then is to find the arrangement of the mirror image (by rotation and translation) that minimises the overlap with respect to the original structure. This is done in a straightforward manner by simulated annealing. Clearly for chiral structures the overlap cannot vanish and by plotting the overlap function we have a measure of the chirality as a function of , as in Fig(6).
Vi Details of Structures for up to 2.71486
Up to , all the structures that are found are of a special character, helpful for their interpretation. Every sphere is in contact with the cylinder surface. Therefore all of their centres lie on an inner cylinder, of diameter Dd, and they touch another cylinder of diameter D2d. They may therefore be considered as the densest packings of spheres on a cylinder.
Clearly inner spheres must appear in the optimal packing when D2d is greater than d, that is , in fact they are first found at , and we reserve them for a separate discussion in the next section. We also treat structures for low separately in section VI.3.
vi.1 Maximal contact structures
We have already noted the special points at which the structure changes continuously, with the formation of new contacts and the breaking of old ones. We have previously called them symmetric structures, but this term is really only meaningful in the diskpacking analysis below. Here we shall apply the term maximal contact instead. For , all of these structures are depicted in Fig (7) and labelled using the phyllotactic notation explained in Appendix B.
They include the simple packings which we have already described and are labelled in the phyllotactic notation, see Appendix B. This relates to the pattern of sphere centres and contacts, “rolled out” on to a plane, and identified with rhombic or triangular lattices. This does not apply to the first of the structures, the straight chain, and is rather confusing if applied to the next two, hence we use it only for the subsequent structures.
We begin with the case of : this is the closepacked, achiral, structure previously described. Each sphere has five contacts (not including contact between a given sphere and the cylinder); these consist of a contact with the neighbouring sphere in the same unit cell and four contacts with spheres in adjacent cells. All of the rest of the maximal contact structures so far considered are composed of spheres with six contacts. Packings of the type are achiral while the rest are chiral.
vi.2 Line slip structures
We now turn to the intervening structures, which are modified versions of the maximal contact packings, adjusted to fit around the cylinder. Such staggered helices were first observed by Pickett et al in their simulations. We label these line slip structures to indicate that the modification consists of the release of half of the contacts along one line separating two of the spiral chains of the symmetric structure, and a consequent relative slip of the two sides. Such line slip arrangements are a feature of hard sphere packings and are evident in the rolled out patterns of Fig 8.
There are in general three possible choices for the direction of the line slip and the results of Fig 8 correspond to those that maintain the highest density as in Fig 3.
We continue to higher values of in section VIII.
vi.3 Structures below , and the square root singularity
We return to examine the structures for low D/d, which are somewhat different. The reason for such a difference will become clearer when we explore the relationship to disk packing, in a later section. These first four structures are shown in Fig 1.
Structure 1 is the elementary case of a straight chain of spheres, like peas in a pod, as shown in Fig 1. This trivial, achiral, structure has a volume fraction of . Similarly structure 2 is obvious and consists of a zigzag planar arrangement of spheres such that the change of the azimuthal angle from one ball to the next is equal to ; each sphere makes contact with one sphere from above and one from below. When the zigzag packing, i.e structure 2, encounters additional contacts between its second neighbours at (so that each sphere now has four contacts), it is forced to take the form of an increasingly twisted spiral (structure 3). By direct numerical calculation one may follow this to its end point at structure 7, but structure 6 intervenes with higher density, so that there must be a transition below D/d = 2.0. This is found at D/d=1.990, strikingly close to 2.0, beyond which the densest packing is a lineslip modification of structure 6. The smallness of the interval in which this is found is largely due to the existence of a squareroot singularity.
The variation of volume fraction, as any of the or structures is approached with increasing , exhibits a squareroot form, whereas it is linear in other cases. That is
(3) 
where the quantities with a subscript 0 belong to a or structure. The derivative plot shown in Fig 2 is intended to make this singular behaviour more evident. The squareroot arises from the special symmetry of these achiral structures.
Vii Relationship to circular disk packings on a cylinder
Many aspects of these results are made more understandable by recourse to the packing of circular disks on a cylinder. We have already noted that the same sequence of ÒsymmetricÓ structures is found. The similarity extends further to the lineslip structures and some of the details of the analytic form of the curve in Fig (3), not yet noted.
This qualitative correspondence attracts our attention to the cylindrical disk packing problem, which we pursue below.
vii.1 Disk Packings
Densest packing of circular disks in a plane places their centres on a triangular lattice where each disk has six contacts; this was rigorously proved about a century ago Aste:2008 (). Cylindrical surface packing of disks with the same density is generated by rolling this pattern on to a cylinder, when possible.
To see how this can be achieved seamlessly let us define a periodicity vector between any pair of disk centres of the triangular lattice packing as shown in Fig (9)a. The region between the start and end of , which is bounded by lines perpendicular to , can then be excised and wrapped onto a cylinder whose circumference is equal to the length of . The resulting structure is a dense, homogenous packing  in the sense that all disk sites are equivalent  which we call a symmetric packing. Any such symmetric packing is characterised by this periodicity vector. In the phyllotactic scheme can be defined by a set of three ordered positive integers , as discussed in appendix B (a simple operational method to assign indices is to count the number of lattice rows that cross the periodicity vector in the three directions, as shown in Fig (9)a and Fig (9)c).
For other values of the cylinder circumference, the symmetric packing may be distorted in some way to wrap around the cylinder. The most obvious adjustment is an affine transformation of the triangular lattice, as shown in Fig (9)b. We illustrate the effect of the transformation by reference to the shaded triangle in Fig (9)a. The affine transformation distorts the equilateral triangle into an isosceles triangle by keeping the length of two adjacent sides fixed and varying the third; consequently, the disk centres form a rhombic lattice which has a lower density (area fraction) compared with the triangular lattice since each disk now has only four contacts (the lowest area fraction corresponds to a square lattice).The transformation may be adjusted to make the length of equal to . The resulting pattern can be wrapped onto the surface of the cylinder. We may call these asymmetric or affine lattice packings.
If is varied the disk centres of such a structure are eventually brought back into coincidence with the sites of a triangular lattice, see Fig (9)c. In this manner the strained structure proceeds from one symmetric packing to another. Since there are three possible choices for the affine transformation, the rules for this process  as reported previously Mughal:2011 ()  are, when applied to the second and third phyllotactic indices, as follows,
where the new phyllotactic indices may have to be rearranged into descending order. Intermediate asymmetric packings may be labelled using rhombic notation , where the ordered indices are the indices common to both the initial and final symmetric states.
We illustrate the use of these rules with an example. In the case of the symmetric packing (illustrated in Fig (9)a) the above rules yield,
Thus the symmetric packing is connected to via the rhombic structure , to via  which are shown in Fig (9)c and Fig (9)b, respectively, and to via .
Although such asymmetric packings are the simplest type of dense structure, intermediate between two symmetric packings, they are not the densest. As reported previously Mughal:2011 (), asymmetric packings are in general superseded by another type of packing involving an inhomogeneous shear of the symmetric lattice in which there is a localised strain or slip along a line (and its periodic replicas).
We can describe such lineslip packings with reference to the symmetric packing . Again there are three possibilities and we describe each of these in turn. In the first case, four lattice rows cross the periodicity vector in the horizontal direction, as shown in Fig (9)a. We allow the final, fourth, row to “slideover” the previous row until the disk centres are once again arranged into a triangular lattice. Thus as illustrated in Fig (9), the line slip in question is intermediate between the the symmetric packing and . By allowing the disks in the final row to slide in the opposite manner we find that the symmetric packing is connected to . In the second case five lattice rows cross the periodicity vector. Four rows are held fixed while the final layer of disks slides over the penultimate row  thus is connected to and . In the third case only a single lattice row crosses the periodicity vector and we find that is connected to and
This localised shear allows the length of to vary continuously. The disks involved in the lineslip have five contacts while the rest of the structure remains symmetric and closepacked, with six contacts for every disk. Clearly the area fraction has a maximum, corresponding to the symmetric packings, while the intermediate lineslip packings have a lower value. As reported previously Mughal:2011 (), there is a simple rule for the closepacked structures that are the end points of lineslip solutions
where the leading numbers denote the direction , or of the lineslip. Again the above rules apply to the second and third phyllotactic indices of a given close packed structure and keep either , or constant. For example, using the above rules, the symmetric packing is connected by a line slip along to or , a line slip along yields or and a line slip along yields or (note in the third case  along  a rearrangement of the phyllotactic indices into descending order was necessary and the new structure is the same as the initial structure).
In Table 1 we use bold numerals to denote the direction of the line slip. In the table we only have cause to mention the first column of lineslip structures, that is: we denote
A full derivation of such connection rules for both the asymmetric and lineslip packings will be given in a future publication.
Fig (10) presents analytical results for the density of these disk packings. Here we see for the first time some of the structures of lower density, as well as those of maximum density. Even at this stage, we find a qualitative resemblance between the Fig (10), the disk packing problem, and the corresponding sphere packings (as shown in Fig (3)).
A simple method for achieving a semiquantitative description of the sphere packing problem relies on the fact that for the range all the spheres in the packing are in contact with the cylindrical boundary. Thus there exists an inner cylinder on which the centres of all the spheres are to be found. The pattern formed by the intersection of the spheres with the surface of the inner cylinder bears a close resemblance to the disk packing problem, which we now consider in detail.
vii.2 Relation to sphere packing
We now explore more closely the quantitative correspondence between disk and sphere packings, and find a way of bringing them as close to each other as is possible without undue complication.
Consider again a packing of spheres, all touching the confining cylinder of diameter . Their centres lie on an inner cylinder of diameter , so that a given centre is located at in cylindrical polar coordinates. The separation of contacting spheres in 3D is but if the inner cylinder is rolled out onto the plane, i.e. a given sphere centre is mapped to the 2D cartesian coordinates , then the distance between the centres of contacting spheres in 2D depends on their mutual orientation.
In general, the required contact rule, in 3D, for a pair of spheres (each having diameter ) whose centres, in cylindrical polar coordinates, are located at and is chan:2011 (),
(4) 
We plot the parametric equation Eq. (4) in Fig. 11, for various values of (where the distance is referred to as the ‘arc length’). We see that, as increases above 2, the ‘surface packing’ resembles the packing of oval disks, which may be approximated by ellipses. For contacting spheres that lie directly in a line parallel to the direction of the cylinder axis, the separation of their centres is still , while for those that lie in the perpendicular plane, it is,
(5) 
For simplicity we take the major and minor axes of the ellipse to have the and , so that the distance on the plane for a pair of contacting spheres is approximated by,
(6) 
Comparing the circumference of the inner cylinder with the length of the periodicity vector (measured in units of the disk diameter ), gives a stretch factor,
(7) 
where the last term is the ratio of the arc distance between contacting spheres, which lie perpendicular to the direction of the cylinder axis. In this way the sphere packing is related to the planar packing of elliptical disks, which is merely a stretched version of a circular disk packing. Given any packing of circular disks in the plane, it may simply be stretched in the direction of the vector , by a factor , and wrapped onto a cylinder of diameter
(8) 
(obtained using Eq. ( 7)), to create a good approximation to a sphere packing.
Thus a correspondence that at first seems distant may be brought much closer. Fig (12) presents the resulting transformed curves for sphere packing density; asymmetric and line slip structures are shown by red (dashed) and black (continuous) lines, respectively. The peaks in the density correspond to symmetric packings. The lower (dashed) heavy black line accentuates the curve of maximum density, as predicted by our analytical approach. This is to be compared with the results of simulated annealing shown by the upper curve. We have extended our simulated annealing study of single layer packings to , as shown by the grey region in Fig (12), by searching for the densest structures in which all the spheres have centres on the surface of a cylinder of diameter . From this we see that the maximum possible density, for single layer packings, is found at and corresponds to the symmetric packing . Beyond this the packing density steadily diminishes since spheres are to be found only on the surface of the cylinder and the interior is empty.
The analytical method presented here gives the correct lineslip solution in the majority of cases. There are, nevertheless, notable discrepancies between the analytical and numerical results. These are due to the simplicity of the transformation used to map the disk packing problem into the sphere packing problem; as a consequence, the analytical results always predict a type (3) lineslip solutions leading up to packing while simulations in fact find a region split between type (2) lineslip followed immediately by a type (3) lineslip. A more accurate (and more complicated) transformation ought to account for this by pushing the type (2) curve above the type (3) solution, for part of this region. This is borne out by the fact that as the diameter of the cylinder increases, and higher order corrections to the transformation diminish in importance, the type (2) solution in the simulations are observed over an increasingly smaller range in the split region.
Viii Structures beyond D/d = 2.71486, with internal spheres.
We are able to pursue the same simulations to higher cylinder diameters, finding the anticipated occurrence of internal spheres, initially at . The structures involved are tabulated in the final section of Table 1. In some cases phyllotactic notation may still be useful, but much of the simplicity of the previous sections is lost, and increasingly so. The extension of the curve of volume fraction, as far as present computing resources allow, is shown on the right hand side of Fig (3).
It is of the same general character as before; that is we observe peaks corresponding to maximal contact structures  which are labelled using the notation (m,n) in table I  and decreasing or increasing yields structures with fewer contacts which can in some way be deformed steadily into the maximal contact packing  these are labelled and respectively.
Packings with internal spheres are depicted in Fig (13) where in each case the middle diagram shows the entire structure viewed sideon. Although the introduction of inner spheres produces more complicated structures there remains an outer shell of spheres contacting the confining cylinder. The rolledout pattern made by these outer spheres is shown to the left of each packing. The arrangement of the inner spheres is depicted in the diagram to the right of each packing, whereby the outer layer is transparent and the inner spheres are in yellow.
The first packing with internal spheres is structure 30, or , which consists of a basic unit cell containing six spheres. Of these one sphere is not in contact with the cylindrical surface and is found at the centre of the cylinder. The remaining five spheres, which are in contact with the cylinder, form a tilted pentagonal ring around the central sphere (i.e. tilted at an angle with respect to the plane perpendicular to the cylindrical axis). The extended structure can be made by translating the unit cell along the cylinder and rotating it by . The result is a chiral packing with an internal chain of noncontacting spheres along the cylindrical axis and an outer layer of touching pentagonal rings. As can be seen, in the corresponding rolled out diagram of the outer layer in Fig (13), there is only a single point of contact between outer rings corresponding to a sphere with four surface neighbours. The location of this sphere alternates by an azimuthal angle of between successive layers.
Some of the spheres in the outer layer have only two contacts with other spheres in the outer layer. However, mechanical stability is assured by their contact with the internal chain of spheres and the average number of contacts per sphere for the structure as a whole is above four.
Increasing reduces the tilt of the pentagonal rings until at we have the maximal contact packing . Thus in structure 31 the spheres in the internal chain are in contact with each other and, compared with structure 30, there are now a greater number of contacts between successive pentagonal rings. Increasing breaks some of these surface contacts and forces the chain of internal spheres to form a twisted zigzag structure. Thus structure 32, or , is chiral and is superseded at by a new type of packing.
Structure 33, or consists of a basic unit cell containing eleven spheres. Of these ten touch the cylindrical surface and are arranged as a pair of pentagonal rings stacked on top of each other. The eleventh, internal, sphere is found above the ten surface spheres and is located at the centre of the cylinder. The extended structure can be described as an alternating sequence of surface spheres followed by an internal sphere. However, on average each sphere is only in contact with other spheres, this number is too low for mechanical stability, which we explain as follows.
As shown by the rolled out diagram of the outer sphere centres, surface spheres from adjacent unit cells are not in contact. Thus the pairs of pentagonal rings are free to rotate, by a certain amount about the cylindrical axis while internal spheres remain fixed in position. Whether structure 33 is achiral or chiral depends on the relative orientation of the surface spheres with respect to each other. No other optimal cylindrical packing yet discovered has this property.
With increasing the pentagonal rings from adjacent unit cells are eventually locked into position at to give . Thus packing 34 is a maximal contact achiral structure, corresponding to a peak in Fig (3). The spheres in the outer layer form a perfect rhombic lattice when rolled out onto the plane and the internal spheres lie along the cylindrical axis. An increase in results in the loss of two contacts (on average) between the surface spheres and the internal spheres, and proceeds downwards in density. As the trend continues we observe a decrease in the separation of neighbouring internal spheres and a modulation in the rhombic surface pattern; this latter symmetry breaking is responsible for structure 35 being a chiral structure.
Packing 36 includes an internal linear chain of noncontacting spheres lying along the central axis of the confining cylinder. These are surrounded by an outer layer of spheres which form a complex chiral structure. As is increased the spheres in the internal chain are brought ever closer to each other until they make contact at , which corresponds to the maximal contact, chiral, structure 37.
A further increase in forces the internal chain of spheres to form a twisted zigzag structure. As a result there is a loss of two types of contacts: the breaking of contacts between spheres in the outer layer (as seen from the corresponding rolled out pattern in Fig (13)) and a break in contact between spheres in the outer and inner layer. Thus there is a decrease in density as we proceed from structure 37 along 38, i.e. as we increase for .
Structure 39 is remarkable in that we find an internal chain of spheres (along the central axis of the cylinder) but the chain is composed of a pair of touching spheres followed by a gap, this is the structure . Increasing produces structure 40, a maximal contact chiral packing , which compared with structure 39 has an increased number of contacts between internalsurface and surfacesurface spheres.
This remarkable sequence of structures would have been difficult to imagine in advance. It is likely to be followed by an equally rich scenario, whenever it becomes possible to pursue higher values of . The structures reported are new, to our knowledge, except in so far as they correspond to dry foam structures mentioned below in section .
Ix Related Experiments
Columnar sphere packings appear in a variety of different experimental contexts, such as the packing of buckyballs inside carbon nanotubes Khlobystov:2004 () and polystyrene spheres inside the pores of silicone membranes Tymczenko:2008 (). The structures that are commonly identified are the straight chain (structure 1 in our table I), zigzag (structure 2 or ), twisted zigzag (structure 3) and structure 5, although more complicated structure have also been observed, e.g. Li:2005 ().
However, care needs to be taken when directly comparing these experimental observations with the simple hard sphere simulations described here, which might well only serve as a first guide of what structures to expect. The minimization of interaction energy may replace the maximization of density as the guiding principle.
For example the polystyrene particles in the experiments of Tymczenko:2008 () have been chargestabilised and thus repel each other. Unlike the situation in our simulations, this results in a preference to sit near the wall of the pores (i.e. the cylinder wall).
Also the comparison of our simulation results with packings of buckyballs might be limited. Experiments by Khlobystov:2004 () show that while it is possible to fill into doublewalled nanotubes down to the ratio (resulting in a simple linear chain), it is not possible to fill singlewalled carbon nanotubes below ( is the inner diameter of the tube in both cases). This reflects the role of the van der Waals forces between and the confining nanotube walls.
ix.1 Present experiments with ball bearings and bubbles
However, we have identified two other experimental systems that provide very good comparison with our simulation data (details to be provided in a followup paper).
The first are metal spheres of a few millimeter in diameter (the type used in “ball bearings”) packed into perspex tubes. The filled tubes were mechanically agitated over an extended period. At the end of this “annealing” procedure this resulted in clearly identifiable ordered sphere packings. These experiments were carried out for a large range of values of D/d (up to 3.15), and resulted in about 16 different structures, of both the maximal contact and the lineslip type.
We have also carried out more extensive experiments with equalvolume gas bubbles (with diameter of a few hundred microns). These are produced by bubbling gas into surfactant solution with the bubble diameter easily variable by adjusting the gas flow.
The bubbles are then gathered into vertically placed cylinders. They maintain a (nearly) spherical shape, even when in contact, up to a column height of a few millimetres, corresponding to the capillary length for the surfactant solution in use.
We were able to produce straight chain, zigzag and a large number of maximal contact structures, and also many structures with internal bubbles, up to and well beyond the point reached by the present simulations. There is also some evidence of the twisted zigzag structure and lineslip structures. The subsequent paper presenting these experimentally obtained structures will include the use of Xray tomography to establish their internal configurations.
ix.2 Ordered dry foam structures
Extensive and detailed results have been published for the various structures that form when equal volume bubbles with diameter exceeding a few millimeters (i.e. larger than the capillary length) are collected in vertical tubes Pittet:1995 (); Tobin:2011 (). In this case, the liquid drains and a dry foam is created  a packing of polyhedral bubbles, with its volume fraction approaching unity.
The surface pattern displayed by these ordered foams is hexagonal (apart from the socalled bamboo structure which is simply an array of parallel soap films). The phyllotactic notation is thus the obvious choice for their classification, at least as long as there are no internal bubbles (see for example Fig. 14).
Restricting ourselves for the moment to this case it would appear that all ordered dry foam structures have corresponding sphere packings (wet foams), i.e. they can be classified by the same phyllotactic notation. These include the straight chain (bamboo), and zigzag structure (called 211 or staircase structure in the foam literature), and all of the 10 maximal contact structures, except structure 25 (541) which has not yet been observed. None of the ordered dry foam structures is of the lineslip type.
Note that since the bubbles are deformable, the average contact number in the case of dry foams is generally higher than in hard sphere packings.
Unlike their sphere packing counterparts, the various dry foam structures are found over ranges of values of (in the foams literature refers to the equivalent sphere diameter of a bubble). Hysteresis plays a large role in the standard experimental procedure, leaving to overlapping ranges of stability for the various structures.
The structure of minimal energy for a given value of may be determined from computer simulations WeairePhelan1996 (); WeaireBradleyPhelan2000 (); SaadatfarEtal08 (); Hutzler:2009 (); Tobin:2011 (). We find these ranges always to be lower than the value of D/d for the corresponding sphere packing, reflecting their much smaller volume fractions. It remains for future research to establish the phase diagram (which may be rather complex) that connects these two limiting cases.
There are also observations of over 20 different ordered dry foam structures with internal bubbles Pittet:1995 (), but the precise internal arrangements have only been identified for the three simplest cases, with the aid of Surface Evolver simulations Brakke (); SaadatfarEtal08 ().
The dry foam with a surface pattern equivalent to structure 28 () has six bubbles in the periodic unit cell, i.e. one more than . A pentagonal dodecahedron in the centre of the cylinder is surrounded by a ring of five bubbles in contact with the cylinder surface. The dodecahedra of neighbouring unit cells are in contact.
There is also a foam structure with a socalled Kelvin cell (tetracaidecahedron) in the centre, surrounded by 6 bubbles in contact with the surface. Again the internal bubbles of neighbouring unit cells are in contact with each other.
The third identified structure with internal bubbles consists of a a total of 13 bubbles in the unit cell, 12 touching the surface, and a bubble with 12 pentagons and 3 hexagons (Goldberg3) in the centre SaadatfarEtal08 (). Interestingly in this case the internal bubbles of neighbouring unit cells are not in contact, similar to the also achiral structure 34. The surface pattern of this foam structure consists of an arrangement of bidisperse hexagons.
X Conclusions
Our simulations have provided detailed results for 40 distinct columnar crystals, of which many are new. They fall into three categories: the very simplest cases, where the sphere diameter is of the same order as that of the cylinder, a wide range of further ”phyllotactic” structures, which may be understood as related to surface packings of disks, and structures that incorporate internal spheres.
Corresponding experimental observations are now available and will be reported in a second paper. Experiments can already be pushed to much larger cylinder diameters, providing further insights and challenges to simulation and interpretation.
Such a further development should motivate a reassessment of the simulation methods to be used; we make no strong claims for the efficiency of that used here. We have necessarily been cautious in using it, stipulating high degrees of convergence and undertaking many repeated runs in search of the optimum, to increase confidence in our conclusions. In parallel studies, Chan chan:2011 () used a sequential deposition method that was highly expeditious in determining optimal structures within the restriction of that procedure, whose effect could not be known a priori. Suitably adapted, it was able to reproduce the structures reported here, up to at least . Future work will include an extension of this deposition approach to higher values of .
Finally we note some possible extensions to this work. An example is the study of disordered hard sphere packings in a cylinder; for such problems the definition of the volume fraction recently given by Chan chan:2011 () may prove useful. A simple problem that may be of significant practical interest is to find the densest packing of hard spheres in a cylinder which is capped on either side by hard walls (i.e. packing in a cylindrical box). Even richer possibilities are offered by packing hard spheres in a cylinder that is capped on both ends by spheres held in a fixed position. Indeed, hard sphere packings that are bounded by templated surfaces on all sides are of significant current interest and are capable of realising a range of different morphologies (including the WeairePhelan structure) Gabbrielli:2012 ().
Xi Acknowledgments
This research was supported by Science Foundation Ireland (08RFPMTR1083) and European Space Agency (MAP AO99108:C14914/02/NL/SH and AO99075:C14308/00/NL/SH). H. K. C. would like to acknowledge support from the Irish Research Council for Science, Engineering and Technology (EMPOWER Fellowship).
Xii Appendix A: Simulation technique
xii.1 Energy function
The simulation is addressed to a cylindrically shaped cell of length and diameter . Contained within this space are points which represent the centres of spheres, each of diameter . If a pair of spheres is sufficiently close that they overlap we account for this using a pairwise potential as described below. A similar overlap potential is used to prevent the spheres from escaping the simulation cell in the radial direction. The final structures will be the densest that we find which have zero overlap energy.
In addition, we impose twisted periodic boundary conditions on the top and bottom of the cylindrical cell as described in the main text, with a twist angle .
We model the overlap potential between spheres using a Hookean, or “springlike”, pairwise interaction between the th and th spheres, which have their centres at and , the interaction energy between spheres is then given by,
(9) 
where is the distance between the centres of the spheres. Note the interaction energy falls to zero when there is no overlap between the spheres.
The interaction energy between the th sphere and the boundary is given by,
(10) 
where .
The twisted periodic boundary conditions on the ends of the cylinder are incorporated as follows: The th particle in the simulation cell has an image at the top and bottom of the simulation cell, the coordinates of these images are given by and , respectively, where is the length of the cylinder and is the twist angle.
Thus the total energy of the system is given by the sum of the spheresphere, sphereboundary and sphereimage interactions.
xii.2 Numerical Method
For a tube of diameter we wish to find the unit cell, composed of spheres, which when rotated and stacked along the tube has the highest volume fraction. In this section we describe the simulation protocol used to achieve this.
For a given and we assign initial starting positions to the spheres and an initial value to the twist angle by using a random number generator. A small initial value for the cylinder length is chosen to insure overlap between the spheres.
Keeping and fixed, we search for the lowest energy arrangement for the spheres by varying their coordinates and the twist angle. This is done using the standard Metropolis simulated annealing algorithm, where for a cluster of spheres the algorithm was run with typically Monte Carlos steps. The temperature of the simulation was decreased linearly. The average displacement of the spheres at each temperature step was chosen by an automatic process to give an acceptance probability of . The results of the simulated annealing are then put through a conjugate gradient routine to ensure that a local minimum has been reached.
The whole process is repeated many times, using a new randomly generated initial configuration, to give confidence that the lowest energy state has been found. From this ensemble we take the lowest energy configuration as the final state for that particular run.
After this first run we perform a subsequent run with a slightly longer cylinder. Since the spheres have more room the energy of the final state is lower compared to the final state of the previous run. Using a divide and conquer approach we are able to establish the cell length at which the energy per sphere in the system falls to zero, in practice this means the value of the energy is within a critical bound which is set to be E/N=. At this point the spheres have a small overlap corresponding to a five decimal place accuracy in the volume fraction (this is deduced by comparing simulation with the analytically derived volume fraction for the circle packing structures).
From this final result we compute the volume fraction which is defined as,
(11) 
where is the volume of one of the hard spheres and is the volume of the simulation cell. Note that the volume fraction depends on both and since different values of yield unit cells with different structures. Thus the whole procedure is repeated for series of different values of , typically up for the structures without internal spheres, before accepting the one which gives the largest volume fraction, as shown in Table 1.
Xiii Appendix B: Continuity of the volume fraction
We demonstrate that there can be no sudden finite discontinuities in the curve describing the maximum volume fraction (density) as is varied. This follows from lower and upper bounds, as described below, which limit the variation in volume fraction in the neighbourhood of any point . We will see that square root singularities observed in the numerical results for the volume fraction are due to the arguments given below.
xiii.1 Lower bound
For increasing an obvious variational argument bounds the density below. As is increased from , the cylinder expands radially. We may take the structure which has maximum density at as a trial structure for . The volume fraction of the structure, at is,
(12) 
and that of the trial structure is,
(13) 
where is the volume of a sphere of diameter and is the average separation between spheres in the direction. Since it follows that,
(14) 
providing a lower bound for .
xiii.2 Upper bound
For decreasing , a more subtle argument bounds the variation of density below by a squareroot function. Decreasing forces the spheres to move radially inwards to avoid contact with the cylindrical boundary. The resulting overlap between spheres can be eliminated by displacement of their centres parallel to the cylinder axis.
Let us index the sphere centres in ascending height using the index  so that for . Given the structure at which has a maximum density, then an extreme case is one where successive spheres all have the same height. Consider in this case a pair of contacting spheres with separation . Reducing the cylinder diameter by a factor of , so that the new diameter is , will force an overlap so that the separation between sphere centres is now . The overlap can be removed by moving one of the spheres vertically a distance , so that their separation is once again . It follows that, . Thus a constant can be chosen so that the following choice is sufficient,
(15) 
In order to eliminate the overlap, the sphere centres are shifted to a new height . The resulting trial structure has the volume fraction,
(16) 
which when combined with Eq. ( 12) results in the following bound
(17) 
Hence we again arrive at a (lower) bound for which goes continuously to as but in this case with a square root form.
References
 (1) G.T. Pickett, M. Gross, and H. Okuyama, Phys. Rev. Lett. 85, 3652 (2000).
 (2) A. Mughal, H.K. Chan, and D. Weaire, Phys. Rev. Lett. 106, 115704 (2011).
 (3) H.K. Chan, Phys. Rev. E 84, 050302(R) (2011).
 (4) A. Meagher, in preperation .
 (5) R.O. Erickson, Science 181, 705 (1973).
 (6) N. Pittet, P. Boltenhagen, N. Rivier, and D. Weaire, Europhys. Lett. 35, 547 (1996).
 (7) D. Weaire, S. Hutzler, and N. Pittet, Forma 7, 259 (1992).
 (8) N. Pittet, N. Rivier, and D. Weaire, Forma 10, 65 (1995).
 (9) P. Boltenhagen, N. Pittet, and N. Rivier, Europhys. Lett. 43, 690 (1998).
 (10) S. Hutzler, D. Weaire, F. Elias, and E. Janiaud, Phil. Mag. Lett. 82, 297 (2002).
 (11) J.H. Moon, S. Kim, G.R. Yi, Y.H. Lee, and S.M. Yang, Langmuir 20, 2033 (2004).
 (12) J.H. Moon, G.R. Yi, and S.M. Yang, J. Coll. Int. Sci. 287, 173 (2005).
 (13) F. Li, X. Badel, J. Linnros, and J. Wiley, J. Am. Chem. Soc. 27, 7262 (2005).
 (14) M. Tymczenko, L.F. Marsal, T. Trifonov, I. Rodriguez, F. RamiroManzano, J. Pallares, A. Rodriguez, R. Alcubilla, and F. Meseguer, Adv. Mater. 20, 2315 (2008).
 (15) M.A. Lohr, A.M. Alsayed, B.G. Chen, Z. Zhang, R.D. Kamien, and A.G. Yodh, Phys. Rev. E 81, 040401(R) (2010).
 (16) A. Khlobystov, D. Britz, A. Ardavan, and G. Briggs, Phys. Rev. Lett. 92, 245507 (2004).
 (17) T. Yamazaki, K. Kuramochi, D. Takagi, Y. Homma, F. Nishimura, N. Hori, K. Watanabe, S. Suzuki, and Y. Kobayashi, Nanotech. 19, 045702 (2008).
 (18) J. Warner and M. Wilson, ACS Nano 4, 4011 (2010).
 (19) T. Aste and D. Weaire, The pursuit of perfect packing, 2nd Ed. (CRC Press, 2008).
 (20) T. Hales, Annals of Math., 2nd Ser. 162, 1065â1185 (2005).
 (21) H. Airy, Proc. R. Soc. 21, 176 (1872).
 (22) S. Tobin, J. Barry, A. Meagher, B. Bulfin, C. O’Rathaille, and S. Hutzler, Coll. Surf. A 382, 24 (2011).
 (23) K.A. Brakke, Exp. Math. 1, 141 (1992).
 (24) D. Weaire and R. Phelan, Phil. Trans.: Math. 354, 1989 (1996).
 (25) D. Weaire, G. Bradley, and R. Phelan, Soft Condensed Matter: Configurations, Dynamics and Functionality (Kluwer Academic Publishers, 2000).
 (26) M. Saadatfar, J. Barry, D. Weaire, and S. Hutzler, Phil. Mag. Lett. 88, 661 (2008).
 (27) S. Hutzler, J. Barry, P. GraslandMongrain, and D. Weaire, Coll. Surf. A. 344, 37 (2009).
 (28) R. Gabbrielli, A. Meagher, D. Weaire, K. Brakke and S. Hutzler, Phil. Mag. Lett. 92, 1 (2012).