## Abstract

We answer the question of whether optimal packings of circles on a sphere are equilibrium solutions to the logarithmic particle interaction problem for values of *N*=3–12 and 24, the only values of *N* for which the optimal packing problem (also known as the Tammes problem) has rigorously known solutions. We also address the cases *N*=13–23 where optimal packing solutions have been computed, but not proven analytically. As in Jamaloodeen & Newton (Jamaloodeen & Newton 2006 *Proc. R. Soc. Lond. Ser. A* **462**, 3277–3299. (doi:10.1098/rspa.2006.1731)), a logarithmic, or point vortex equilibrium is determined by formulating the problem as the one in linear algebra, , where *A* is a *N*(*N*−1)/2×*N* non-normal configuration matrix obtained by requiring that all interparticle distances remain constant. If *A* has a kernel, the strength vector is then determined as a right-singular vector associated with the zero singular value, or a vector that lies in the nullspace of *A* where the kernel is multi-dimensional. First we determine if the known optimal packing solution for a given value of *N* has a configuration matrix *A* with a non-empty nullspace. The answer is yes for *N*=3–9,11–14,16 and no for *N*=10,15,17–24. We then determine a basis set for the nullspace of *A* associated with the optimally packed state, answer the question of whether *N*-equal strength particles, , form an equilibrium for this configuration, and describe what is special about the icosahedral configuration from this point of view. We also find new equilibria by implementing two versions of a random walk algorithm. First, we cluster sub-groups of particles into patterns during the packing process, and ‘grow’ a packed state using a version of the ‘yin-yang’ algorithm of Longuet-Higgins (Longuet-Higgins 2008 *Proc. R. Soc. A* (doi:10.1098/rspa.2008.0219)). We also implement a version of our ‘Brownian ratchet’ algorithm to find new equilibria near the optimally packed state for *N*=10,15,17–24.

## 1. Introduction

We consider solutions of the optimal packing problem of *N* equal circles on the surface of a sphere (Conway & Sloane 1998) and ask, if a particle of strength is placed at the centre of the *i*th circle (*i*=1,…,*N*), is this a relative equilibrium configuration for a system of interacting particles with logarithmic interactions (Bergersen *et al.* 1994; Jamaloodeen & Newton 2006)? In other words, we ask whether solutions of the Tammes (1930) are also point vortex equilibria on the sphere (Newton 2001; Newton & Sakajo 2009), and, if so, for what particle strengths. Surprisingly, solutions to the Tammes problem are only rigorously known for values of *N*=2–12 and 24, and for each of these values, except *N*=5, the solution is unique.

The notions of optimal distribution of points on a sphere and interacting particle equilibria can be unified by viewing each as an extremizer of the ‘Riesz-s interaction energy’ **∇***E*_{s}, where
1.1
The prime indicates that we exclude the terms *i*=*j*. In the limit *E*_{0} corresponds to the logarithmic (point vortex) interaction considered in this paper which arises naturally as a discrete model for incompressible fluid dynamics (Newton 2001), where each singularity corresponds to a delta-function representation of the vorticity field at that point. Minimizing *E*_{0} corresponds to maximizing the product of interparticle distances. This particle interaction problem has been considered extensively in the literature, and we refer simply to Aref *et al.* (2003), Newton (2001) or Newton & Chamoun (2009) for a relevant overview. When , the problem corresponds to the spherical packing problem of Tammes (1930) (see Clare & Kepert 1986, 1991) in which the minimum distance between circles is maximized. This can be seen by noting that as , the dominant term in equation (1.1) will be the one with minimum |*x*_{i}−*x*_{j}|. Thus, to minimize , one must maximize this minimum distance. In this sense, the logarithmic interaction problem and the spherical packing problem lie at opposite ends of the allowable values for *s*, with Coulomb interactions, *s*=1 (known as the Thomson 1904, or plum-pudding model of the atom), in between. The literature on the Thomson problem is large, and we refer simply to the recent work of Wales & McKay (2009) for an entry. We note that in all cases, the Euler constraint: *F*−*E*+*V* =2 must be enforced, where *F* is the number of faces of the converged polyhedron, *E* is the number of edges, and *V* =*N* is the number of vertices. See Saff & Kuijlaars (1997) for more discussions on extremizers of the Riesz-s energy, Bergersen *et al.* (1994) for logarithmic interactions, and Edmundson (1992), Erber & Hockney (1991), Glasser & Every (1992) for earlier works on this and related problems, the more recent relevant work of Cohn & Kumar (2007), and Alben (2008) for an interesting consideration of the equilibrium shapes of line singularities on the sphere.^{1} We also note that the energy, as written in equation (1.1), assumes that each of the particles has equal strength, which means . We will consider both this case and the more general case of whether there exist ANY particle strengths for which the optimally packed state is a point vortex equilibrium.

In §2, we review the analytically known solutions to the Tammes optimal packing problem of circles on a sphere for *N*=3–12,24. An entertaining introduction to this problem can be found in Bagchi (1997). We then review the criteria we use for identifying point vortex equilibria for the logarithmic interaction problem. Along with this comes a discussion of the singular spectral decomposition of the configuration matrix *A*. Then, in §3, we apply the methods of §2 to solutions of the optimal packing problem. In §4, we adapt a method introduced in Longuet-Higgins (2008) called the ‘yin-yang’ method, in which packed states are ‘grown’ using random walks and growing circles (shrinking spheres). During formation, preferential clustering into subgroups of particles is enforced. In some cases, the clusterings and growth lead directly to equilibria, whereas in others, a ‘Brownian ratchet’ adjustment of the positions of the particles is necessary to find a nearby equilibrium state. We demonstrate this and exhibit a catalogue of new examples of equilibria obtained using combinations of these random walk methods.

## 2. Solutions to Tammes problem

### (a) *N*=3–12,24

The analytically known solutions to the Tammes problem are listed in column 1 of table 1 and a (numerical) listing for larger values of *N* can be found in Clare & Kepert (1986, 1991), which is by no means the most comprehensive listing known. In fact, one can find excellent websites, such as http://www2.research.att.com/njas/ and http://www-wales.ch.cam.ac.uk/ which contain a wealth of computational information regarding these and related particle interaction problems. For our purposes, we focus first on the cases for which there are known analytical proofs—this is only for the values *N*=2–12,24 (we do not include the *N*=2 antipodal state in our table). One of the standard measures associated with optimal packing solutions is *l*(*r*), the diameter of the packed circles. In each case, this is listed in the third column of table 1. We take a sphere of unit radius, hence *r*=1. A second measure of the optimally packed state is the proportion *p* of the sphere’s surface that is covered by the spherical caps. A formula for this, given in Longuet-Higgins (2008), is , where *N* is the number of vertices and *α* is the angle subtended by an edge at the centre of the sphere. These values are also listed in table 1.

For *N*=4,6 and 12, the optimally packed state is the corresponding Platonic solid for that value of *N*, namely the tetrahedron (*N*=4), octahedron (*N*=6) and icosahedron (*N*=12). The Platonic solid for *N*=8 is the cube, which is not the optimally packed state. Instead, the staggered cube where the top square is rotated by *π*/4 is optimal (which can be viewed as a nice example of a single vortex street configuration). For *N*=5, there is actually a continuum of optimally packed states, we list the two at the extremes of the continuum range. One is the octahedron minus a vertex, while the other has a vertex at each of the poles, with the three others evenly distributed around the equator. For *N*=11, the optimal packing is the icosahedron minus a vertex, while for *N*=24, it is the snub cube. For values *N*=7,9 and 10, the optimal states are difficult to describe in simple geometric terms, but we refer the reader to Danzer (1986) for references and values of *l*(*r*) and *p* in table 1 for comparison purposes. In the first columns of figures 1 and 2 we show optimal configurations for each case *N*=3–11. The points shown in the Northern and Southern hemispheres are placed at the centres of the packed circles, which correspond to the vertices of embedded polyhedra. In figure 3, we show the icosahedral state *N*=12, a configuration that plays a special role in many biological packing settings (see Zandi *et al.* 2004; Johnston *et al.* 2009; Torquato & Jiao 2009) and as we will describe later, is also special with respect to our equilibrium considerations. Notions here of ideas ‘balanced symmetry’ and ‘balanced geometries’ (e.g. Wales 1990, 1998; Cohn *et al.* 2010) are potentially relevant, but they typically only apply to situations where the particle strengths are all equal. Our next goal is to see whether these optimally packed states are also point vortex equilibria.

### (b) Point vortex equilibrium criterion

Since our characterization of relative equilibria is based on the fact that they can be viewed as fixed points for the dynamical equations governing intervortical distances, convenient starting points are the dynamical equations of motion for *N*-point vortices moving on the surface of the unit sphere (see Newton 2001; chapter 4), which written in Cartesian coordinates are given by
2.1
The vector **x**_{α} denotes the position of the *α*th vortex the strength of which is given by *Γ*_{α}∈R. The prime on the summation indicates that the singular term *β*=*α* is omitted and initially, the vortices are located at the given positions , (*α*=1,…,*N*). The denominator in equation (2.1) is the intervortical distance, *l*_{αβ}, between vortex *Γ*_{α} and *Γ*_{β} since . As described in Newton & Sakajo (2009), the equations (2.1) can be written in Hamiltonian form, with logarithmic interaction energy *H*, where
2.2
By exponentiating both sides of this equation, one can see that minimizing *H* corresponds to maximizing the product of distances between pairs of particles (with powers *Γ*_{α}*Γ*_{β}).

The evolution equations for the relative distances are (Jamaloodeen & Newton 2006)
2.3
where . Here the ^{′′} means that the summation excludes *γ*=*α* and *γ*=*β*. *V*_{αβγ} is the volume of the parallelepiped formed by the vectors **x**_{α},**x**_{β} and **x**_{γ}
The sign of *V*_{αβγ} can be positive or negative depending on whether the vectors form a right- or left-handed coordinate system. Our criteria for a relative equilibrium is that the interparticle distances remain constant, namely we look for fixed points of the evolution equation (2.3). The relative equations of motion yield necessary conditions for relative equilibria,
2.4
or equivalently,
2.5
for each value of *α*,*β*=1,…,*N*. Based on the fact that equation (2.5) is linear in the vortex strengths, we write it as a linear matrix system
2.6
where is the vector of vortex strengths, and *A* is the *M*×*N* configuration matrix (*M*=*N*(*N*−1)/2) the entries of which, given by the terms *V*_{αβγ}*d*_{αβγ}, encode the geometry of the configuration. Thus, a necessary condition for a relative equilibrium is given by
2.7
where *A*^{T} is the transpose of *A*, and *A*^{T}*A* is the covariance matrix associated with the configuration. To render these conditions sufficient, one must check that the *l*_{αβ}’s are constant for all time which in practice is most easily checked numerically by simulating the converged configuration via the evolution equation (2.1).

Once the configuration is chosen, the entries of *A*, as given in equation (2.5), can be obtained and one must (i) determine if *A* has a kernel, and, if it does, (ii) obtain a basis set for its nullspace. Any that lies in this nullspace is then an allowable choice which makes that configuration a relative equilibrium. It is then a straightforward calculation to check whether the special choice lies in the nullspace.

## 3. Spectral decomposition of optimally packed states

### (a) Singular value decomposition of *A*

To carry out the steps outlined above in a comprehensive way, we spectrally decompose *A* via its singular value decomposition from which we obtain an optimal basis set for the nullspace of *A* (Trefethen & Bau 1997). We note that there is a redundancy in the full set of equations (2.4), which will be reflected in the singular value structure of *A*. We obtain the *N* singular values *σ*_{i} and corresponding left and right singular vectors , by solving the coupled linear system
3.1
where . The vector is called the left-singular vector associated with *σ*_{i}, while is the right-singular vector and it is clear from this equation that when *σ*_{i}=0, the right singular vector is an element of the nullspace of *A*. In terms of these, the matrix *A* has the factorization
3.2
Since *M*>*N*, the first *N* columns of *U* are the left-singular vectors , and the remaining *M*−*N* columns are chosen to be orthonormal so that *U* is orthogonal (*U*^{T}*U*=*I*). Then *U* is an *M*×*M* orthogonal matrix, *V* is an *N*×*N* orthogonal matrix the columns of which are the right-singular vectors , *Σ* is an *M*×*N* matrix with the *N* singular values on the diagonal and zeros off the diagonal. Equivalently, multiplying the first equation in (3.1) by *A*^{T}, the second by *A*, and uncoupling the two, we obtain
3.3
which expresses the fact that the singular values squared are the eigenvalues of the square covariance matrices *A*^{T}*A*, *AA*^{T}. We write these eigenvalues as *λ*_{i}≡(*σ*_{i})^{2}. The decomposition (3.2) expresses *A* as a linear superposition of the rank-1 matrices , (*i*=1,…,*N*) with weighting determined by the singular values *σ*_{i}.

The singular spectrum of each of the best packing configurations is shown in figures 1, 2, and 3. In these plots, we have normalized each of the eigenvalues of the covariance matrices so that they lie in the range from zero to one and can be interpreted either as probabilities, or as the percentage of energy contained in each mode. The normalized eigenvalues are given by
3.4
where *k* is the number of non-zero singular values, hence the rank of *A*. The sequence of numbers then forms a discrete stochastic sequence which uniquely characterizes each relative equilibrium. As discussed in Shannon (1948) in the context of communication theory, the *entropy*, *S*, is obtained from the sequence via the definition:
3.5
Equation (3.5) provides a convenient scalar measure of how the rank-1 matrices in equation (3.2) are distributed in forming the configuration matrix, and thus can be thought of as a measure of ‘disorder’ of the pattern. In particular, if all of the weighting is in one mode, then *A* has rank-1 and the Shannon entropy is minimized—its value is zero. The configuration matrix in this case can be viewed as the ‘least disordered’ in terms of how its rank-1 modes are distributed. On the other hand, if each mode has equal weighting, the entropy is maximum—its value is . In this case, the configuration matrix is the ‘least ordered’ in terms of how the rank-1 modes are distributed. In addition, this measure of maximum entropy increases monotonically with *k*. Interpreted slightly differently, the Shannon entropy of the configuration can be thought of as a *measure of how close the configuration matrix is to one of low rank*. The lower the entropy, the closer the matrix is to a rank-1 matrix, whereas the higher the entropy, the closer it is to one of full rank. Table 1 lists the nullspace dimension of *A* associated with each of the optimally packed states, and we show the Shannon entropy associated with the non-zero spectrum of *A* for each.

### (b) Results for *N*=3–12,24

The results can be summarized as follows. For values of *N* in the range 3≤*N*≤12, the optimally packed states give rise to a configuration matrix *A* with a non-empty nullspace except for the value *N*=10 where the nullspace is empty. For *N*=24, the optimally packed state also gives rise to an empty nullspace. Therefore, particle strengths can be chosen for each of the optimally packed configurations corresponding to *N*=3–9,11,12 to make the configuration a point vortex relative equilibrium. It is a straightforward exercise to show that lies in the nullspace, and thus that equal strength particles distributed at the vertices of an icosahedron forms a relative equilibrium. For *N*=3,4,6, the dimension of the nullspace of *A* is *N*, hence in each of those cases, the corresponding system with equal strength particles also forms an equilibrium. However, for *N*=5–9,11, the vector does not lie in the nullspace of *A*. In table 1, we list the Shannon entropy (3.5) for each configuration, as well as the Hamiltonian energy level (2.2) corresponding to the equal particle strength configuration if it is available, or if there is a one-dimensional nullspace which makes the energy uniquely defined.

### (c) The yin-yang algorithm and results for *N*=13,14,16

For values of *N* other than those listed in table 1, we need a numerical algorithm to find the optimally packed state. An interesting and useful algorithm is that introduced by Longuet-Higgins (2008). Motivated by an application of polyhedral theory to biological questions associated with the growth of the outer sheath of certain virus molecules (see Crick & Watson 1957; Klug & Finch 1965; Cohn & Kumar 2009; Johnston *et al.* 2009; Fejer *et al.* 2010), Longuet-Higgins introduced a random walk method which he dubbed the yin-yang method. First, *N* spherical caps are placed with centres at random points on the unit sphere. The radii of each of the caps are equal and initially small enough so that no two caps overlap. During the ‘yin’ phase, the centre of each cap, in turn, is given a small random step. After each step, if two or more caps overlap, the step is cancelled, if not, it is retained. After *N* such ‘yin’ steps (one for each cap), in the ‘yang’ phase, the minimum angular distance between the centres of any two of the caps is calculated. Then, the radii of each of the caps is allowed to increase by an amount proportional to this minimum distance, and the yin-yang process continues until two or more of the caps touch. The final ‘packed’ state is thus obtained by an iterative sequence of ‘growth’ and ‘random adjustment’, and compared with known maximally packed states for certain given values of *N*. In figures 4, 5, and 6 we show the known (Clare & Kepert 1986, 1991) optimally packed state grown using the yin-yang algorithm. These figures correspond to values of *N*=13,14,16. In each case we show (i) the configuration of points on the sphere, (ii) the singular spectrum associated with the equilibrium, and (iii) the particle strengths associated with the equilibrium as obtained by finding the right singular vector associated with the zero singular value of the configuration matrix *A*.

## 4. Equilibria via preferential clustering and random walks

### (a) The Brownian ratchet method

We now ask a fundamentally different question, one associated more with the *formation* of packed and equilibrium states. The assembly process of spherical macro-molecules, particularly virus molecules, has been the subject of intensive research since Crick & Watson’s (1957) paper in which they hypothesized that many viruses were made up of symmetrical arrays of subunits with overall icosahedral symmetry. Questions then immediately centred not only on the symmetry properties of these final configurations, but on how they arose. An excellent introduction to some of the fundamental issues can be found in ‘Virus Structure and Assembly’, edited by Casjens (1985) while discussions of packed states in the context of proteins and DNA molecules can be found in Liang & Dill (2001) and Purohit *et al.* (2003).

First, for the cases *N*=10,24, where the optimally packed state is NOT an equilibrium, we implement our Brownian ratchet algorithm, as detailed in Newton & Sakajo (2009) to find a ‘nearby’ equilibrium. Briefly stated, using the optimally packed state as an initial condition, we (i) give each point a small random perturbation, where the size is scaled with the smallest singular value; then (ii) if the smallest singular value decreases after the perturbation, we keep the new configuration, otherwise we go back to the previous configuration and repeat the previous step. The algorithm typically converges to a new equilibrium by ‘ratcheting’ the smallest singular value down to zero, thus, our method can be viewed as the one that finds an equilibrium by randomly ‘adjusting’ the optimally packed state with a sequence of smaller and smaller adjustments until an equilibrium is found. Shown in figures 7 and 8 are equilibria obtained near the optimally packed states for *N*=10,24. Each case shows the optimally packed state compared with the ‘nearby’ equilibrium configuration, along with the singular spectrum, particle strengths and Frobenius norm difference between the two configurations.

### (b) Preferential clustering using the yin-yang method

Next, we ask, if the *N*-particles are constrained during formation to lie in ‘preferentially clustered’ sub-units, how will this effect the structure of the converged states? This question was introduced to the authors in the previously mentioned recent work of Longuet-Higgins (2008). For our purposes, we are particularly interested in systems that form in the presence of icosahedral symmetries (see Zandi *et al.* 2004; Torquato & Jiao 2009). Thus, we first cluster sub-groups of *N*=12 particles during the yin-yang phase of formation and see if this leads to an equilibrium configuration different from the icosahedral configuration.

Let us first consider two groups of flower petals with five satellite particles, that is, two central vortices each surrounded by five equally spaced particles around a spherical cap. Putting point vortices at the centre of the spherical caps and aligning the two centres of the flower petals at the north and the south poles respectively, we have a configuration where each of the five satellite points are equally arranged along the lines of latitudes. Singular values corresponding to the converged configuration yield five singular values that are 0.2 and the other seven are zero. This indicates that the configuration, after the yin-yang phase, is the icosahedral state as shown in figure 3. Thus, preferential clustering into two subgroups of ‘flower-petal’ arrangements seems to yield the optimally packed icosahedral state.

By contrast, we apply the numerical algorithm to the groups of flower petals with three satellite particles, namely three groups of 4. Figure 9*a* shows spherical caps at the end of the yin-yang phase of the algorithm, which is not an equilibrium configuration since the smallest singular value is not zero. We then use this configuration as the initial state for the Brownian ratchet scheme detailed in Newton & Sakajo (2009), hence we randomly search for a nearby state, keeping the sub-unit structure intact. The converged equilibrium after this Brownian ratchet step is shown in figure 9*b*, which is not the icosahedral state, but some nearby state. Indeed, only one singular value corresponding to the equilibrium is zero as we see in figure 9*c*. The entropy of the equilibrium is 2.139352, which is higher than that of the icosahedron. The right eigenvectors to the zero singular value are shown in figure 9*d*, with which the energy of the configuration is 6.03151×10^{−2}.

Finally, we consider *m* groups of *k* particles that are equally spaced along a circle without a particle in the centre. For *N*=12, possible pairs of (*m*,*k*) are (3,4), (4,3) and (6,2). For (*m*,*k*)=(3,4) and (4,3), we obtain no convergent final state via the Brownian ratchet step, which means there is no equilibrium near the configuration that attains the maximum packing of the sphere. On the other hand, for (*m*,*k*)=(6,2), we have an equilibrium state as shown in figure 10*a* just by the ying-yang algorithm. The configuration consists of two staggered ring structures of six point vortices, the configuration matrix of which has one-dimensional null space (figure 10*c*). Figure 10*d* shows that the six point vortices of each ring cluster have the same strength, but the two clusters have strengths of opposite signs. Its entropy is 2.164777, which is lower than that of the icosahedron. The energy is 6.03151×10^{−2}. We can find the similar staggered ring states for 2 groups of *M* vortices equally spaced along the line of latitudes as in figure 10.

In figures 11, 12, and 13 we show a catalogue of equilibria obtained by the preferential clustering method using *N*=9,15,16,18,20,21 particles. The caption above each of the configurations indicates how the clusterings were enforced (+1 denotes a point vortex at a pole) while the points on the sphere are the final equilibrium states obtained. In each, we show the singular spectrum associated with the equilibria, each of which has a one-dimensional nullspace.

## 5. Discussion

We show in this paper that an optimally packed configuration of circles on a sphere is not necessarily an equilibrium configuration for interacting particles with logarithmic interaction energy (i.e. maximizing the minimum distance is not necessarily the same as maximizing the product of distances), nor does there seem to be a simple way of predicting if it will be for a given *N* value, and general particle strengths, based on symmetry considerations alone. The cube, for example, is NOT an optimally packed state, but it is an equilibrium configuration (see Jamaloodeen & Newton 2006), whereas the optimally packed states for *N*=10 and 24 (snub cube) are NOT equilibrium configurations. Furthermore, if we restrict ourselves to systems with equal strength particles, the optimally packed states for values of *N*=5,7,8,9,11 do not yield equilibria. The icosahedral configuration *N*=12 is special from this point of view—it is the only non-degenerate (rank greater than zero) analytically known optimally packed state that is also an equilibrium configuration for equal strength particles. There seems to be no simple symmetry based argument, such as those discussed in Wales (1990, 1998), or Cohn *et al.* (2010) leading to a clear correspondence rule which associates optimally packed states with equilibria, only a calculation involving the nullspace and singular spectrum of the corresponding configuration matrix *A* gives the answer.

Furthermore, we show that there are a whole host of interesting equilibria obtained by clustering particles into subgroups during the formation process that do not converge to the optimally packed state. This is demonstrated simply for the case *N*=12,24 where the icosahedron and snub cube are the optimally packed configurations, yet preferential clustering and Brownian ratchets lead to new equilibria—constraints imposed during equilibrium formation can determine which final state is obtained. As a practical matter, since a great deal is known (at least numerically) about optimally packed configurations of circles on spheres for impressively large values of *N* (in general, much more is known about these configurations than is known about interacting particle equilibria), a systematic investigation of these configurations and nearby configurations as candidate point vortex, Coulomb or Riesz-s equilibria seems like a potentially fruitful direction to pursue. The stability of these many configurations is, of course, a largely unexplored area.

## Acknowledgements

P.K.N. acknowledges financial support from the National Science Foundation through Grant NSF-DMS-0804629. T.S. is partially supported by JSPS grant 21340017 and by JST PRESTO Award.

## Footnotes

↵1 For a recent conference devoted to state-of-the-art developments in this field, see http://www.math.vanderbilt.edu/~optimal2010/.

- Received July 13, 2010.
- Accepted November 4, 2010.

- This journal is © 2010 The Royal Society