## Abstract

We describe a Brownian ratchet scheme that we use to calculate relative equilibrium configurations of *N* point vortices of mixed strength on the surface of a unit sphere. We formulate it as a problem in linear algebra, *A*** Γ**=0, where

*A*is a

*N*(

*N*−1)/2×

*N*non-normal configuration matrix obtained by requiring that all inter-vortical distances on the sphere remain constant and

**∈**

*Γ*^{N}is the (unit) vector of vortex strengths that must lie in the null space of

*A*. Existence of an equilibrium is expressed by the condition det(

*A*

^{T}

*A*)=0, while uniqueness follows if Rank(

*A*)=

*N*−1. The singular value decomposition of

*A*is used to calculate an optimal basis set for the null space, yielding all values of the vortex strengths for which the configuration is an equilibrium and allowing us to decompose the equilibrium configuration into basis components. To home in on an equilibrium, we allow the point vortices to undergo a random walk on the sphere and, after each step, we compute the smallest singular value of the configuration matrix, keeping the new arrangement only if it decreases. When the smallest singular value drops below a predetermined convergence threshold, the existence criterion is satisfied and an equilibrium configuration is achieved. We then find a basis set for the null space of

*A*, and hence the vortex strengths, by calculating the right singular vectors corresponding to the singular values that are zero. We show a gallery of examples of equilibria with one-dimensional null spaces obtained by this method. Then, using an unbiased ensemble of 1000 relative equilibria for each value

*N*=4→10, we discuss some general features of the statistically averaged quantities, such as the Shannon entropy (using all of the normalized singular values) and Frobenius norm, centre-of-vorticity vector and Hamiltonian energy.

## 1. Introduction

We introduce a ‘Brownian ratchet’ scheme to obtain relative equilibrium configurations of *N*-point vortices on the surface of a sphere. Using the criterion that all inter-vortical distances remain constant, we find relative equilibria as fixed points of the evolution equation of the inter-vortical distances. Since these equations are linear in the vortex strength vector ** Γ**∈

^{N}, we formulate the problem as one in linear algebra, namely

*A*

**=0, where**

*Γ**A*∈

^{M×N}is a non-normal configuration matrix with rows. Thus, we seek arrangements of particles on the sphere for which det(

*A*

^{T}

*A*)=0, or equivalently, for which

*A*has a non-trivial null space. When this condition is satisfied, we obtain the vortex strength vector

*Γ**a posteriori*by finding a basis set for the null space of

*A*.

The method used to produce an equilibrium is based on using the *k* smallest singular values of the configuration matrix as a ‘ratchet’, which we drive to zero by a random walk algorithm on the sphere. The number of singular values that are zero corresponds to the dimension of the null space (which we call the *degree of heterogeneity* of the configuration) and thus the number of basis vectors needed to span the subspace of ^{N} in which the vortex strength vector lies. The decomposition method based on the null space of the configuration matrix was introduced in Jamaloodeen & Newton (2006) and used to determine all vortex strengths for which the Platonic solid configurations with a point vortex at each vertex form an equilibrium. Subsequently, the Brownian ratchet scheme coupled with the use of the singular value decomposition of the configuration matrix was developed by Newton & Chamoun (2007) and used to study equilibrium configurations in the planar *N*-vortex problem. The singular value decomposition gives rise to the ‘optimal’ basis set in which to represent the vortex strength vector and also produces a characteristic ‘distribution’ of singular values that allows us to calculate other important quantities, such as the Shannon entropy, a quantity borrowed from the field of communication theory (Shannon 1948), and the size of the configuration, based on the Frobenius norm. The equilibria described in this paper all have configuration matrices with one-dimensional null spaces and hence a unique vector of vortex strengths and, typically, they have no discernible symmetries. Previous results on relative equilibria of point vortices on the sphere, such as those of Kidambi & Newton (1998), Lim *et al*. (2001), Laurent-Polz (2002), Cabral *et al*. (2003) or Newton & Shokraneh (2006), are much more restrictive. Typically, they assume the vortex strengths to be equal (hence without loss of generality unity), or occur in equal and opposite pairs in the case where *N* is even. This is also true of studies of equilibrium distributions of particles on a sphere with more general interaction laws, such as those of Erber & Hockney (1991), Edmundson (1992), Glasser & Every (1992), Bergersen *et al*. (1994), Altschuler *et al*. (1997), Saff & Kuijlaars (1997) and Altschuler & Perez-Garrido (2005, 2006). By allowing the vortex strengths to take on any value, we show that the set of relative equilibrium configurations is far richer than previously realized. Despite this fact, we also show that, surprisingly, many of the key statistical quantities associated with ensembles of these equilibria remain well behaved. For general background on the *N*-vortex problem, see Newton (2001).

Our paper is organized as follows. In §2, we describe the basic tool we use to construct relative equilibria on the sphere, namely the singular value decomposition of the configuration matrix associated with each equilibrium. The distribution of these singular values (properly normalized) gives rise to a scalar quantity that characterizes the equilibria—the Shannon entropy of the configuration matrix. In §3 ,we describe the Brownian ratchet algorithm that we use to calculate the collection of equilibria for each *N*. In particular, we describe our random walk algorithm on the sphere and how it is used to home in on configurations of particles that produce a configuration matrix with a non-trivial null space. We show examples of typical relative equilibria for *N*=4, 6, 8, 10. We also detail the convergence properties of the Brownian ratchet scheme. In §4, we discuss several statistical properties of the unbiased ensembles, including the statistically averaged Shannon entropy, Frobenius norm, Hamiltonian energy and centre of vorticity. Section 5 contains a discussion of our key findings.

## 2. Decomposing the pattern

The evolution equations for *N*-point vortices moving on the surface of the unit sphere, written in Cartesian coordinates, are given by(2.1)The vector *x*_{α} denotes the position of the *α*th vortex whose strength is given by *Γ*_{α}∈. The prime on the summation indicates that the singular term *β*=*α* is omitted and, initially, the vortices are located at the given positions *x*_{α}(0)∈^{3}, *α*=1, …, *N*. The denominator in (2.1) is the inter-vortical distance, *l*_{αβ}, between vortex *Γ*_{α} and *Γ*_{β} since . As described in Newton & Shokraneh (2006), equation (2.1) has two conserved quantities associated with it, the Hamiltonian energy(2.2)and the centre-of-vorticity vector(2.3)

The evolution equations for the relative distances are(2.4)where . Here, the ″ means the summation excludes *γ*=*α* and *γ*=*β*. *V*_{αβγ} is the volume of the parallelepiped formed by the vectors *x*_{α}, *x*_{β} and *x*_{γ},Notice that the sign of *V*_{αβγ} can be positive or negative, depending on whether the vectors form a right- or left-handed coordinate system. The relative equations of motion yield sufficient conditions for relative equilibria,(2.5)

### (a) The configuration matrix approach

Using condition (2.5) in (2.4) gives the equation for the relative equilibria(2.6)for each value of *α*, *β*=1, …, *N*. Based on the fact that (2.6) is linear in the vortex strengths, we write it as a linear matrix system(2.7)where is the vector of vortex strengths and *A* is the *M*×*N* configuration matrix (*M*=*N*(*N*−1)/2) whose entries, given by the terms *V*_{αβγ}*d*_{αβγ}, encode the geometry of the configuration. Without loss of generality, we normalize the vector of vortex strengths to have unit length, hence(2.8)Thus, we seek configurations so that(2.9)in which case *A* is rank deficient, and has a non-trivial null space. We seek a basis set for this subspace of ^{N}. In all the cases considered in this paper, Rank(*A*)=*N*−1, hence the vortex strength vector is unique up to the ± sign.

### (b) Singular value decomposition

The optimal basis set for the null space of *A* is obtained by using the singular value decomposition of the matrix (see Trefethen & Bau 1997). We obtain the *N* singular values *σ*_{i} and corresponding left and right singular vectors , by solving the coupled linear system(2.10)where . The vector *u*_{i} is called the left-singular vector associated with *σ*_{i}, while *v*_{i} is the right-singular vector. In terms of these, the matrix *A* has the factorization(2.11)Since *M*>*N*, the first *N* columns of *U* are the left-singular vectors *u*_{i}, 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 whose columns are the right-singular vectors *v*_{i} and *Σ* is an *M*×*N* matrix with the *N* singular values on the diagonal and zeros off the diagonal,(2.12)Equivalently, multiplying the first equation in (2.10) by *A*^{T}, the second by *A* and uncoupling the two, we obtain(2.13)which express the fact that the singular values squared are the eigenvalues of the square covariance matrices *A*^{T}*A* and *AA*^{T}. We write these eigenvalues as *λ*_{i}≡(*σ*_{i})^{2}. The decomposition (2.11) expresses *A* as a linear superposition of the rank-one matrices , *i*=1, …, *N*, with weighting determined by the singular values *σ*_{i}. Its optimality is seen by the fact that the *m*th partial sum, defined as(2.14)provides the best rank-*m* approximation to *A*, as measured by the Frobenius norm. In other words, any rank-*m* matrix has the property that , where ‖.‖_{F} denotes the Frobenius norm defined as .

### (c) Shannon entropy

To understand how the rank-one modes are distributed, it is useful to normalize 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(2.15)where *k* is the number of non-zero singular values, hence the rank of *A*. The sequence of numbers forms a discrete stochastic sequence that 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(2.16)Equation (2.16) provides a convenient scalar measure of how the rank-one matrices in (2.11) 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 one 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-one modes are distributed. On the other hand, if each mode has equal weighting, the entropy is maximum—its value is ln(*k*). In this case, the configuration matrix is the ‘least ordered’ in terms of how the rank-one 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-one matrix, whereas the higher the entropy, the closer it is to one of full rank. We mention also that low entropy distributions are less ‘robust’ to perturbations than high entropy ones. As discussed in Newton & Chamoun (in press), generic perturbations to a given configuration will tend to *increase* the entropy of a base configuration, i.e. spread out the distribution among the modes.1 If the distribution is already spread out in the base state (i.e. a high entropy base state), the perturbation has a smaller effect than if the modes are clustered among only a small number. See Newton & Chamoun (in press) for more comprehensive discussions of these ideas.

## 3. The Brownian ratchet idea

Our method of obtaining relative equilibria is based on a Brownian ratchet scheme that we implement by a diffusion process in the plane, which we then map to the unit sphere. The terminology we use is borrowed from the biological literature in which molecular motors are known to extract energy from their surrounding ‘heat bath’ and rectify it via a ratchet mechanism. See Riemann (2002) for a comprehensive recent review. For us, the heat bath is provided by a random walk algorithm on the sphere, while the ratchet that rectifies this motion is the smallest singular value of the configuration matrix that we drive to zero. The random walk problem on the sphere is interesting in its own right, and has been studied in the past by Brillinger (1997), who considered the motion of a particle on the unit sphere heading towards a specific destination but subject to random deviations, which he modelled as a diffusion process with drift. His motivation was to model the trajectories of certain marine mammals, and in so doing he obtained quantitative formulae for expected travel times to a spherical cap, as well as forms for limiting distributions. Indeed, before this work, Kendall (1974) was interested in modelling the navigation of birds and used a pole-seeking Brownian motion model to partially explain their behaviour. An early and quite general work on random walk models on the sphere and on more general Riemannian manifolds is that of Roberts & Ursell (1960).

### (a) Random walk on a sphere

The random walk algorithm on the sphere is the engine that drives our ratchet scheme, so we describe it first. As shown schematically in figure 1, we start with an initial ‘seed’ point (*θ*_{0}, *ϕ*_{0}) on the sphere. From this point, the random walk is computed in three simple steps:

First, we obtain a sample point in the plane from the two-dimensional Gaussian distribution, for which we compute the polar coordinate representation (

*r*,*φ*).Next, with a scale factor

*ϵ*(typically taken as*ϵ*=0.01), we rescale the point as (*ϵr*,*φ*) and then map it to a corresponding point on the surface of the unit sphere centred around the North Pole so that the point is represented by (*θ*_{ϵ},*ϕ*_{ϵ})=(*ϵr*,*φ*) in spherical coordinates.Finally, we rotate the point so that the North Pole maps to the original point (

*θ*_{0},*ϕ*_{0}), while (*θ*_{ϵ},*ϕ*_{ϵ}) maps to the new ‘diffused’ point (*θ*_{1},*ϕ*_{1}).

The process is then iterated to obtain each subsequent point (*θ*_{n+1}, *ϕ*_{n+1}), starting with (*θ*_{n}, *ϕ*_{n}) as a seed. Here, the procedure is implemented for a collection of particles initially clustered around the North Pole (those marked as circles), and South Pole (those marked as pluses), as shown in figure 2. As the particles evolve, they gradually diffuse over the surface of the sphere, eventually giving equal probability of finding a circle particle or a plus particle in any fixed two-dimensional spherical sector.

### (b) The ratchet scheme in practice

For each *N*, we seek configurations of particles on the unit sphere for which (2.9) is satisfied, hence Rank(*A*)<*N*. We find these configurations with the following ratchet algorithm.

First, we distribute

*N*points randomly on the surface of the unit sphere and calculate the configuration matrix*A*, finding its smallest singular value*σ*_{min}.We then allow each particle to execute one random step on the sphere in order to produce a new configuration matrix , along with its smallest singular value .

If , we keep the new configuration, otherwise we discard it.

The process is repeated until drops below a certain pre-determined threshold, which we typically choose to be

*O*(10^{−10}). This ‘converged’ configuration is what we call a relative equilibrium.We then compute a basis set for the null space in order to find the corresponding vortex strengths.

Typical convergence plots are shown in figure 3. Figure 3*a* shows the decay of the smallest singular value (squared) as a function of the step number for *N*=6, 8, 10, plotted on a log–log scale. In most cases, convergence is rapid. Figure 3*b* shows the actual path of one of the point vortices making up the configuration from its initial point to its final (converged) point (marked by a plus) on the sphere. Note that the vortex meanders initially before it homes in rather directly to its final location, which need not be nearby the initial location. As a general remark, we note that the singular values of a matrix are relatively insensitive to perturbations of the matrix (see Trefethen & Bau 1997); hence, we expect that the converged positions of the vortices are not far from the exact equilibrium positions when the smallest singular value is below *O*(10^{−10}).

### (c) Gallery of relative equilibria for *N*=4, 6, 8, 10

Typical examples of relative equilibria found this way are shown in the panels of figure 4 for *N*=4, figure 5 for *N*=6, figure 6 for *N*=8 and figure 7 for *N*=10. In each figure, we present a panel of six distinct relative equilibrium configurations showing both the vortex positions in the Northern and Southern hemispheres. In each case, the intersection of the centre-of-vorticity vector ** J** (as defined in (2.3)) with the unit sphere marked by an ‘×’. All of the cases treated in this paper have one-dimensional null spaces; hence, unique vortex strength vectors that we normalize to unity. Note that all of the configurations are manifestly asymmetric, a topic discussed in Newton & Chamoun (in press). Examples of asymmetric equilibria are indeed rare, the first discussion of this can be found in Aref & Vainchtein (1998).

In figure 8, we show histograms (for large collections of equilibria as discussed in §4) of the length of the ** J** vector for the cases

*N*=4, 6, 8, 10. In all cases, the peak is near the unit value, indicating that most of the states making up the ensemble can be described as not too different from a single dominant vortex of near unit strength resting near the tip of the centre-of-vorticity vector, with the remaining

*N*−1 weaker vortices distributed asymmetrically around the surface of the sphere. In all cases, the

*N*vortices have mixed signs and the spread around the most likely state tightens as

*N*increases, indicating that the limiting configuration (constrained to have Rank=

*N*−1) is a single vortex of unit strength resting at the tip of the centre-of-vorticity vector.

Likewise, histograms of the Hamiltonian energy (2.2) are shown in figure 9, and, in each case, the peak value is zero with a spread that tightens with increasing *N*. This limiting configuration suggests a relatively uniform distribution of points around the sphere with vortex strengths of mixed sign.

## 4. Statistical properties

In contrast to classes of equilibria obtained by other methods (see Aref *et al*. (2003) for a comprehensive overview), the approach described in this paper is capable of generating large *unbiased* ensembles of equilibria. This is both due to the random initial conditions used to start each Brownian based search and to the random search algorithm that is capable of finding all relative equilibria, not just those with prescribed symmetries or specific vortex strengths. Thus, it makes sense to use these ensembles to produce statistically averaged quantities that characterize the equilibria. We discuss some of these properties in this section.

### (a) Ensemble averages

For each value of *N*=4→10, we generate an ensemble of equilibrium configuration matrices, denoting each member of the ensemble *A*^{(j)}, with corresponding right null vector *Γ*^{(j)}. The initial sample size for each case is nominally *M*=1000, which we double to *M*=2000 by including both ±*Γ*^{(j)}. The singular values for the *j*th realization are denoted by and their corresponding left and right singular vectors are denoted by and , , respectively. We define the ensemble average of the collection of configuration matrices(4.1)as well as the ensemble averages of the singular components(4.2)(4.3)The standard deviation of each quantity is denoted with double brackets 〈〈.〉〉. We denote the averaged normalized values(4.4)(4.5)with standard deviations . We then define the Shannon entropy of the *j*th member of the ensemble to be(4.6)with ensemble average(4.7)and standard deviation 〈〈*S*〉〉.

### (b) Statistical properties

Here, we summarize the main results based on an analysis of the ensemble averages for the cases *N*=4, 5, 6, 7, 8, 9, 10. Table 1 shows the ensemble-averaged properties of the normalized singular values, listed in decreasing order, for the case *N*=10. For each of the 10 singular values, we show the maximum value in the ensemble , the minimum value , the sample mean and the sample standard deviation for *M*=1000. The smallest sample average is with a gap of 10 orders of magnitude between it and the next smallest value . The size of the smallest singular value, the gap between it and the next smallest and the steady decrease of the convergence curve shown in figure 3*a* give us confidence that we are in close proximity to an equilibrium configuration. Figure 10 shows the distribution of the normalized singular values for *N*=4, 6, 8, 10. A noteworthy feature is that the shape of the distribution for the final two cases *N*=8 and 10 is quite similar, indicating convergence to a fixed distribution as a function of *N*.

In table 2, we show the statistical properties of the averaged Shannon entropy and Frobenius norms for *N*=4, 5, 6, 7, 8, 9, 10. These quantities, shown as a function of the sample size *M*, are depicted in figures 11 and 13. It is interesting to note from figure 11 that the spacing of the converged values is quite regular, indicating an underlying scaling law. Indeed, in figure 12 we show the ensemble-averaged Shannon entropy values, shown in table 2, plotted as a function of *N* on a log–log scale. The data show power-law scaling of the form , with *α*∼0.305683 and *β*∼0.671424, as obtained via a least-squares fit to the data. In figure 14, we show histograms of the total vortex strength of each equilibrium. We note the tendency for to cluster at the extreme values ±1 in agreement with the observation that the histograms of ‖** J**‖ in figure 8 cluster around 1. The ‘pure translation’ case appears to be quite rare, although there are examples of pure translational equilibria in the samples.

## 5. Discussion

The Brownian search scheme described in this paper, based on a linear algebra formulation of the problem (in contrast to the classical variational approach used, for example, in Campbell & Ziff (1979)), offers an unbiased approach for finding all of the relative equilibrium configurations of point vortices on the sphere, regardless of their stability properties or symmetries. For the range of values of *N* used in this paper, the convergence properties of the algorithm were adequate—for larger values of *N*, we expect convergence to be more sluggish. The richness of the class of relative equilibria allowed us to use them as microstates from which to extract information on the macroscopic level via ensemble averages. There are two main findings:

The length of the centre-of-vorticity vector, ‖

‖, clusters near 1, as shown in the histograms of figure 8, while the total vorticity associated with each member of the ensemble, as expressed by , tends to cluster at the extreme values of ±1, as shown in the histograms in figure 14.*J*The averaged Shannon entropy scales very nearly as , with

*β*∼2/3. This quantity reflects the averaged distribution of the normalized singular values shown in figure 10 as a function of*N*and provides a scalar measure of the relative weighting of the rank-one components, , constituting the equilibrium ‘pattern’, as encoded in the configuration matrix and expressed in (2.11).

The first conclusion provides evidence that the macroscopic average vorticity can be thought of as one single vortex of unit strength, with either clockwise or anticlockwise circulation, discretized, in a sense, by the point vortices in their relative equilibrium configuration. Since this macroscopic state is in agreement with statistical results reported by mean-field theory using collections of equal strength vortices moving dynamically on the sphere or via Monte Carlo simulations (see the recent monograph of Lim & Nebus (2006)), it suggests that using the full family of relative equilibria (presumably most of them unstable) offers a useful and rich enough set of microscopic building blocks from which to extract meaningful macroscopic information. The second conclusion, we believe, is unexpected as there is no *a priori* reason for the averaged quantities to follow any clean scaling law. Indeed, as shown in figure 13, the ensemble-averaged Frobenius norms do not exhibit clear scaling features. As a final remark, we point out that the methods and conclusions reached in this paper are also relevant in treating the classical problem of optimally distributing *N*-charged electrons on the surface of a conducting sphere, an unsolved problem with a long history (e.g. Erber & Hockney 1991; Edmundson 1992; Glasser & Every 1992; Bergersen *et al*. 1994; Altschuler *et al*. 1997; Saff & Kuijlaars 1997; Altschuler & Perez-Garrido 2005, 2006), and listed by Smale (2000) as one of the outstanding mathematical problems for the next century.

## Footnotes

↵An alternative way of conveying this same idea is contained in Shannon's original 1948 publication. He proves that any

*averaging*operation on of the form , where , with all*a*_{ij}≥0, will*increase S*. Other properties of this logarithmic quantity are also discussed in Shannon's paper.- Received May 14, 2008.
- Accepted September 23, 2008.

- © 2008 The Royal Society