## Abstract

A theory capable of producing equilibrium configurations of point vortices in the plane, along with a numerical scheme to compute them, is described. The theory is formulated as a problem in linear algebra where one must find solutions to the matrix equation , where *A* is the (1/2)*N*(*N*−1)×*N* non-normal configuration matrix obtained by requiring that all intervortical distances remain fixed, and are the *N*-vortex strengths. For existence of an equilibrium, *A* must have a non-trivial nullspace. We consider the singular values of *A*; when this has one or more zero singular values, the nullspace of *A* is non-empty and an equilibrium exists for some choice of ** Γ**. New equilibrium configurations are found numerically by randomly depositing

*N*points in the plane, which generically gives rise to a configuration matrix

*A*with empty nullspace. Using the sum of squares of the

*k*smallest singular values of

*A*as a ‘ratchet’, we ‘thermally fluctuate’ the configuration, allowing each point to execute a random walk in the plane, retaining only those configurations which reduce this quantity at the next step. The configuration is thus driven to one with nullspace (

*A*)=

*k*>0. These converged states are not necessarily nearby their initial configurations, typically they are asymmetric, and often we can drive the same initial state to several different equilibria. A reverse-ratchet method is also described, which can produce initial conditions that would evolve to a specified equilibrium state. Once a converged final state is achieved, the full singular value decomposition of

*A*is used to calculate an optimal basis set for the nullspace of

*A*and thus all allowable

**. The**

*Γ**distribution*of the singular values gives important information on the size of each equilibrium state (as measured by Frobenius norm), their distance from each other (spacing and density) and how far a randomly chosen system of

*N*points in the plane is from the nearest equilibrium configuration with a specified rank, as well as its Shannon entropy.

## 1. Introduction

The study of equilibrium patterns of point vortices in the plane has a long and interesting history, dating back to the 1800s with the work of Mayer (1878), Thomson (1878), Warder & Shipley (1888) and Wood (1898) on configurations of floating magnetic devices which self-assemble into symmetric geometric configurations reminiscent of those that would form by a collection of straight, rectilinear thin vortex filaments, i.e. point vortices. These objects were once thought of as physical models of ‘vortex atoms’ by Thomson (soon to become Lord Kelvin). In an elegant reincarnation of these early experiments, Grzybowski *et al*. (2000) produced new and more robust symmetric patterns, some of which had first been seen in the famous Los Alamos catalogue of two-dimensional vortex patterns produced by Campbell & Ziff (1978). For type II superconductors, vortex lattice configurations were first predicted by Abrikosov (1957). The first direct image of such a lattice can be found in the paper by Essmann & Trauble (1967), which shows a very clear ‘perfect’ triangular lattice configuration of point vortices, as well as some which are not so perfect. More recent experiments on Bose–Einstein condensate vortex lattices have produced still more elaborate, precise and smaller-scale configurations (see Abo-Shaer *et al*. 2001), captured some of their interesting dynamics in the formation process (Tsubota *et al*. 2002), under forcing (Adhikari & Muruganandam 2002), uncovered states with multiply quantized vortex sites (Lundh 2002) and lattices with ‘broken symmetries’, i.e. which exhibit defects or dislocations in the field. Particularly interesting and relevant to this work are the experiments reported by Engels *et al*. (2002, 2003), in which they show a vortex lattice with no discernible symmetries that reforms after being perturbed. This work has reinvigorated a classical subject and brought to light new issues that remain poorly understood. Why, for example, are symmetric states far more prevalent than asymmetric ones and what causes the dislocations, that sometimes form, which break the high degree of symmetry associated with a lattice? Two-dimensional vortex lattices can also be produced using magnetically confined non-neutral plasmas (Durkin & Fajans 2000) and superfluid helium (Yarmchuk *et al*. 1979), and although each of these continuum systems is governed by a different mean-field equation, when the vortices at the lattice sites are sufficiently compact, their azimuthal structure becomes irrelevant and they can be modelled more directly as a linear superposition of velocity fields, which are radially symmetric around each site and decay with distance *Γ*_{β}.*r*^{−a}, with *Γ*_{β}∈ being the strength of the *β*th particle and *a*≥1. The velocity field at the lattice site of the *β*th particle is thus a superposition of the velocities induced by the other *N*−1 particles. When *a*=1, the system has been rigorously derived, under suitable assumptions, as a particle approximation to the two-dimensional Euler equations of incompressible fluid flow (Marchioro & Pulvirenti 1994), and it is generally accepted that such a decay law also generically models these other continuum systems, at least to leading order and when all the particles are stationary. We refer the reader to the fairly comprehensive recent review on these so-called ‘vortex crystals’ (Aref *et al*. 2003).

Notably lacking in all of these studies are examples of completely asymmetric patterns and *heterogeneous* states made up of vortices whose strengths are not equal or do not sum to zero. Though it was long suspected that patterns exhibiting no particular symmetries existed, a commonly held view was that they were elusive owing to their lack of stability (Tkachenko 1966*a*,*b*). Indeed, most of the numerical methods used to find equilibria, such as those of Campbell & Ziff (1978), relied heavily on driving the free energy to a local minimum that was close to an initial symmetric configuration. A common assumption is that the vortices are spaced evenly on nested rings from a central point (Aref 1982); hence, they can be thought of as a family of nested polygonal patterns of various shapes and sizes (Lewis & Ratiu 1996). The only paper that focuses on producing asymmetric states is that of Aref & Vainchtein (1998) who found such states of *equal strength* point vortices using a clever parameter continuation strategy.

Our goal in this paper is to formulate an interacting particle theory as a problem in linear algebra and derived from the *N*-point vortex equations governing a wide class of vortex lattice systems. Equally important to the development of such a theory is a practical numerical scheme capable of producing these states; hence, we describe what we call a ‘Brownian ratchet scheme’ (in analogy to molecular and biological motors that rectify Brownian motion) to home in on new *N*-vortex equilibrium patterns, showing examples for *N*=6, 8 and 10. Our method gives insights into how prevalent an equilibrium configuration is with respect to all possible *N*-vortex configurations, why asymmetric patterns are far less prevalent than symmetric ones (not relying on stability arguments) and why equilibrium states with pre-assigned vortex strengths (such as all being equal) are even less prevalent than more complex heterogeneous states. Since the search algorithm we use is driven by unbiased Brownian motion, the method has an equally likely chance of finding any existing equilibrium configuration and does so relatively quickly.

## 2. The equilibrium problem

First, we pose the equilibrium problem as a constrained variational problem based on an augmented Hamiltonian and then we formulate it as a classical problem in linear algebra.

### (a) Variational formulation

The equations of motion for *N*-point vortices in complex form, which assume radial symmetry at each site and have a velocity field *Γ*_{β}.*r*^{−1}, are given by Newton (2001) and Aref *et al*. (2003) as(2.1)Here, *z*_{α}=*x*_{α}+i*y*_{α} denotes the position of the *α*th vortex in the complex plane, *Γ*_{β}∈ denotes its strength and the prime on the summation indicates that the singular term *β*=*α* is omitted. Equation (2.1) arises from a Hamiltonian(2.2)as is well known (Newton 2001). If we assume that the entire configuration of point vortices moves as a rigid body, then each vortex can, in principle, translate and rotate, which gives rise to the ansatz(2.3)where *V* is the (complex) translational velocity, *ω* is the (real) rotational frequency, and these are the same for each of the *N* vortices making up the configuration. By substituting this in the equations of motion (2.1), we obtain the system of equations that must be solved to obtain an equilibrium(2.4)When the vortex strengths are first given, by solving this algebraic system, we simultaneously find the vortex positions *z*_{α}, the translational velocity *V* and the rotational frequency *ω* of the rigid structure.

This algebraic system can also be obtained by extremizing the augmented Hamiltonian(2.5)where *V*=*u*+i*v* and *ω* play the role of Lagrange multipliers and the constraints are the conservation of linear impulse and angular impulse (Aref *et al*. (2003) for discussions on this variational principal due to Kelvin). To see this, we differentiate equation (2.5) with respect to *x*_{α} and *y*_{α,}(2.6)(2.7)We also know that the equations of motion for the system in the complex plane are(2.8)Using equation (2.8) in equation (2.6), equation (2.7) yields(2.9)Combining these relations gives(2.10)or(2.11)where *V*≡*u*+i*v*. These are precisely the equations that guarantee that the system of *N* vortices moves as a rigid collection of points in the plane, with translational velocity *V* and angular velocity *ω*. Thus, equilibrium configurations of *N*-point vortices in the plane arise as extremizers of the interparticle Hamiltonian energy (2.2), subject to the constraints that the linear and angular impulse be conserved. The stability of the configuration thus depends on whether it sits at a local minimum (stable), maximum (unstable) or saddle point (unstable) on the appropriate energy landscape.

### (b) Linear algebra approach

Alternatively, if we form the equations for the intervortical distances from (2.1), i.e. the equations governing , we obtain (Aref *et al*. 2003)(2.12)Here, *A*_{αλβ} is the *signed* area subtended by the three vortices located at points *z*_{α}, *z*_{λ} and *z*_{β}. If they appear in clockwise order, then the area is considered positive, otherwise it is negative. The double prime indicates that the singular terms *β*=*α* and *β*=*λ* are omitted. Any configuration of *N*-point vortices such that(2.13)is an equilibrium state; hence, a necessary and sufficient condition for the *N*-point vortices to form a relative equilibrium is(2.14)A key observation is that this set of equations is *linear* in the vortex strengths *Γ*_{β} (in contrast to the interparticle Hamiltonian energy (2.2) which is quadratic in the vortex strengths) and hence can be written as a matrix equation(2.15)where is the vector of vortex strengths and *A* is an (1/2)*N*(*N*−1)×*N* non-normal (*A*^{T}*A*≠*AA*^{T}) matrix, before any reduction which takes advantage of the conservation of the Hamiltonian and linear and angular impulse. The dimension of the nullspace of the matrix *A* is denoted by nullity(*A*), the rank by rank(*A*) and the range of *A* by *R*(*A*). We call *A* the *configuration matrix*, and the dimension of the nullspace of *A* represents the configuration's degeneracy. Since rank(*A*)+nullity(*A*)=*N*, non-trivial equilibria exist iff rank(*A*)<*N* (in this case, *A* is said to be rank deficient), or equivalently nullity(*A*)>0. Multiplying equation (2.15) by *A*^{T} yields the matrix equation(2.16)where *A*^{T}*A* is now a square, symmetric *N*×*N* matrix, called the configuration's *covariance* matrix. The condition for equation (2.16) to have a non-trivial solution is det(*A*^{T}*A*)≡0, which is equivalent to saying nullity(*A*)>0. Equivalently, if *A*^{T}*A* has one or more zero eigenvalues, then there is a non-trivial nullspace of *A* and ** Γ** can be chosen so that an equilibrium exists for that configuration. If the dimension of this nullspace is one, then there is a

*unique*choice for the vector of vortex strengths

**(up to multiplicative constant) for which the configuration is an equilibrium. Thus, each equilibrium configuration can be naturally categorized according to the dimension of the nullspace of its corresponding configuration matrix (Jamaloodeen & Newton 2006; Newton & Chamoun 2006) as well as its Hamiltonian energy level. This formulation gives immediate insights into several interesting issues, such as the ‘likelihood’ of finding an equilibrium state and thus their density in the space of all possible configurations.**

*Γ*Figure 1*a* shows the smallest eigenvalue of *A*^{T}*A* corresponding to 10^{6} random arrangements of *N*=10 points in the plane. Since none of the eigenvalues are zero, none of the configurations are equilibria, indicating that a generic arrangement of *N* points in the plane will not be an equilibrium for any choice of ** Γ**1. Alternatively, if we first choose

**and ask which configurations lead to equilibria (as is the common practice), these will be a subset of those produced by random arrangements as shown in figure 1, making them even more rare. Figure 1**

*Γ**b*shows the square root of sum of squares of the smallest two eigenvalues, and our conclusions are the same for the case with two- and higher-dimensional nullspaces—the problem of finding an equilibrium state is akin to that of finding a needle in a haystack (figure 2).

## 3. The singular value decomposition

### (a) Summary of important properties of the singular value decomposition

The appropriate tool for understanding and characterizing the nullspace, range and rank of a non-normal matrix is the singular value decomposition (SVD; Golub & Van Loan 1996). The *N* singular values, *σ*^{(i)} (*i*=1, …, *N*), of the *M*×*N* real matrix *A*, are non-negative real numbers that satisfy(3.1)where *u*^{(i)}∈^{M} and *v*^{(i)}∈^{N}. 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 (spectral decomposition)(3.2)where *U* is an *M*×*M* unitary (i.e. *U*^{T}*U*=1) matrix whose columns are the vectors *u*^{(i)}; *V* is an *N*×*N* unitary matrix with columns given by *v*^{(i)}; and *Σ* is an *M*×*N* matrix with non-negative numbers on the diagonal and zeros off the diagonal(3.3)The singular values can be ordered so that *σ*^{(1)}≡*σ*^{(max)}≥*σ*^{(2)}≥⋯≥*σ*^{(min)}≥0 and of course one or more may be zero. As is evident from multiplying the first equation in (3.1) by *A*^{T} and the second by *A*,(3.4)i.e. the singular values squared are the eigenvalues of the matrices *A*^{T}*A* or *AA*^{T}, which have the same eigenvalue structure, while the left-singular vectors *u*^{(i)} are the eigenvectors of *AA*^{T} and the right-singular vectors *v*^{(i)} are the eigenvectors of *A*^{T}*A*. From equation (3.1), we also note that the right-singular vectors *v*^{(i)} corresponding to *σ*^{(i)}=0 form a basis for the nullspace of *A*. From a practical point of view, it is well known (Golub & Van Loan 1996) that numerical algorithms which calculate the singular values of *A* are much more stable and accurate than those that calculate the eigenvalues of *A*^{T}*A*, particularly when calculating small or zero values.

The main use here is that we seek configuration matrices with one or more singular values that are zero, and the SVD provides an explicit and optimal representation of the range and nullspace of the matrix *A*. In particular, the right-singular vectors *v*^{(i)} corresponding to the singular values that are zero span the nullspace of *A*, while the left-singular vectors *u*^{(i)} corresponding to the non-zero singular values span the range of *A*. The rank of *A* equals the number of non-zero singular values and since rank(*A*)+nullity(*A*)=*N*, we know that the number of zero singular values equals the dimension of the nullspace of *A*. In addition, the ranks of *A*, *A*^{T}*A* and *AA*^{T} are the same, and *A*^{T}*A* and *AA*^{T} have the same non-zero eigenvalues and nullspaces.

If the rank of *A* is *r*, then there are *r* non-zero singular values ordered so that *σ*^{(1)}≥⋯≥*σ*^{(r)}≥*σ*^{(r+1)}=*σ*^{(r+2)}=⋯*σ*^{(N)}=0 and one obtains the following representation of the matrix *A*:(3.5)In this way, the singular value decomposition expresses *A* as a sum of *r* rank-one matrices and also provides an optimal method of approximating *A* by another matrix of reduced rank. In particular, we can define a rank *k*<*r* approximation to *A* by keeping only the first *k* terms of equation (3.5),(3.6)This matrix is sometimes referred to as the *k*-truncated SVD. It can be proved that *A*_{k} is the optimal reduced rank approximation to the matrix *A*, meaning that any other rank *k* matrix approximation to *A* has greater error, as measured by their Frobenius-norm difference.

### (b) Variational aspects of the SVD

It is perhaps worth pointing out, in view of Kelvin's variational characterization of equilibria based on the augmented Hamiltonian, that the singular values can also be characterized variationally. Consider the quadratic function formed from the configuration matrix(3.7)Here, ** u** and

**are unit vectors in**

*v*^{M}and

^{N}, respectively. It can be proved that

*σ*attains its maximum when

**≡**

*u*

*u*^{(1)}and

**≡**

*v*

*v*^{(1)}, i.e. when they correspond to the left- and right-singular vectors of

*A*. The maximum value achieved by the function

*σ*corresponds to the largest singular value of

*A*,

*σ*(

*u*^{(1)},

*v*^{(1)})≡

*σ*

^{(1)}. Likewise, the next largest singular value maximizes (3.7) when

**and**

*u***are required to be orthogonal to**

*v*

*u*^{(1)}and

*v*^{(1)}. The

*k*th largest singular value maximizes (3.7) when

**and**

*u***are required to be orthogonal to**

*v*

*u*^{(1)}and

*v*^{(1)},

*u*^{(2)}and

*v*^{(2)}, …,

*u*^{(k−1)}and

*v*^{(k−1)}.

In addition, it is also interesting to note that if we define the Shannon entropy, *E*, of the system in terms of normalized eigenvalues of the covariance matrix(3.8)where the probabilities *P*_{i} are defined in terms of the normalized eigenvalues(3.9)then it is a standard fact that the basis set obtained using the SVD *minimizes* the Shannon entropy, which, to some extent, dictates the distribution of the singular values.

## 4. The Brownian ratchet scheme

Effective tools to find a needle in a haystack are algorithms driven by Brownian motion, so that all portions of the configuration landscape are sampled in an unbiased way. Here, we use what we call a Brownian ‘ratchet’ scheme. The idea behind a Brownian ratchet is simple—it is a device that sits in a heat bath and rectifies the non-equilibrium thermal fluctuations to generate motion in a preferred direction. One can read surveys by Doering (1995, 1998) and Reimann (2002). In analogy with these studies, we ratchet the Brownian motion of the particles, allowing it only in the direction that decreases the smallest singular value (or sum of the smallest *k* singular values) of *A* until it drops below a pre-assigned threshold. In this way, we drive the configuration towards an equilibrium (i.e. one or more zero singular values), but we do not know which values to assign ** Γ** until we derive a basis set for the nullspace of

*A*corresponding to the converged state. The number of singular values that are simultaneously driven to zero correspond to the dimension of the nullspace and thus determine whether or not the equilibrium configuration is unique (up to a multiplicative constant) with respect to the choice of

**.**

*Γ*### (a) The Brownian ratchet algorithm

The algorithm consists of four main steps.

First, we randomly deposit

*N*points in the plane in an unbiased way and compute the*N*singular values of the configuration matrix*A*. These can be ordered and denoted*σ*_{1}≡*σ*_{(max)}≥*σ*_{2}≥⋯≥*σ*_{N}≡*σ*_{(min)}≥0. The minimum singular value,*σ*_{N}, is positive with probability one. In practice, we cycle through 10^{6}random configurations and choose the one with the smallest*σ*_{(min)}as our initial state in order to speed up the convergence in later steps.After picking an initial state, we allow each of the

*N*points to execute an unbiased random walk in^{2}, and we compute the singular values of*A*at each random step. The singular values, after the*n*th step, are denoted , and the step size is scaled so that it is linearly proportional to the smallest, hence decreases in size as the configuration gets closer to an equilibrium.To find a configuration with a one-dimensional nullspace, at the (

*n*+1) step, we only keep the new arrangement if the minimal singular value decreases from that of the previous step, i.e. if . Otherwise, we discard the configuration, produce a new random arrangement and repeat this step. For equilibria with*k*-dimensional nullspaces, we drive the scalar quantity to zero in a similar fashion.When (or equivalently

*δ*_{k}) is below a certain predetermined threshold, i.e.*δ*_{k}<*δ*_{threshold}, the algorithm has converged. Typically, we take*δ*_{threshold}∼10^{−10k}.

In practice, the algorithm converges fairly quickly to equilibrium configurations, particularly if *δ*_{k} associated with the initial guess is small. Figure 3 shows some convergence properties of the algorithm. Figure 3*a* shows the evolution of the smallest singular value squared of *A* as the configuration homes in on a converged state for *N*=6, 8 and 10. Figure 3*b* shows a magnified view of the convergence path of one of the vortices making up the configuration. We also note that there is nothing that prevents us from reversing the ratchet, as described in step (iii), and only retaining new configurations at each step that *increase* the smallest singular value. This reverse-ratchet procedure can be used when starting from a *specific* equilibrium configuration in order to find a path leading to a generic starting point with empty nullspace. The convergence path identified in this manner is also a convergence path to the equilibrium from that starting point using the forward-ratchet procedure.

For us, the truncated SVD plays an important practical role. As stated in step (iv) of the numerical procedure, we are searching for configurations which lead to matrices *A* which have one or more zero singular values. In practice, we say the configuration has converged if one or more of the singular values is smaller than 10^{−5}. We then approximate this ‘converged’ configuration matrix (say of rank *r*) with the optimal rank *k*<*r* matrix *A*_{k} as in (3.6) (i.e. by setting *σ*^{(k+1)}=*σ*^{(k+2)}=⋯*σ*^{(r)}=0) to compute a basis set for the nullspace of *A*_{k} and hence the vortex strength vectors ** Γ** for which the converged configuration is an equilibrium. Figure 2 shows examples of equilibria for

*N*=6, 8 and 10 corresponding to configuration matrices with one-dimensional nullspaces. Note that these converged states show no underlying symmetries (figure 3).

Of course, the method can be used to find symmetric states as well, as shown in figure 4. For this, we use the reverse-ratchet procedure, as described earlier, increasing the smallest singular value at each step instead of decreasing it. The case shown in figure 4*a* starts as an *N*=6 triangular relative equilibrium state with a one-dimensional nullspace (Newton & Chamoun 2006). As the configuration evolves to its final asymmetric state (filled circles), the smallest singular value increases to its final value *O*(1), showing that the configuration made up of filled circles is not an equilibrium. However, if this asymmetric configuration were used as an initial state for the forward-ratchet scheme, one realization of an ensemble of random runs would converge to the symmetric equilateral triangle. Likewise, figure 4*b* starts with an *N*=7 polygonal relative equilibrium with a two-dimensional nullspace (Newton & Chamoun 2006) and uses the reverse-ratchet method to evolve to a final asymmetric state (filled circles) that is not an equilibrium, but if it were used as the initial state for the forward-ratchet scheme, then it would converge to the polygonal state in one of its random runs.

### (b) The Frobenius norm and distribution of singular values

A key characterization which allows us to quantify the size of each equilibrium configuration as well as their distances from each other is the Frobenius norm of a matrix, which is the sum of squares of its singular values. Thus, the Frobenius norm of a rank *r* matrix *A* is given by(4.1)The distance between this matrix and its *k*-truncated SVD is(4.2)We use this characterization to calculate the size of each of the equilibrium configurations, how far each converged state in figure 2 is from its initial state, as well as how far each of the respective equilibria are from each other. Table 1 shows the Frobenius-norm size for each of the equilibrium states (column 1), the distances between the initial and converged states (column 2), the Hamiltonian energy value (2.2) for each of the asymmetric equilibria (column 3) and the distances between each of the two converged states for each *N* (column 4). Table 2 shows all of the singular values for each of the equilibrium states from figure 2, with one singular value in each case nearly zero. The size of the next smallest singular value gives a measure of how far these equilibria with one-dimensional nullspaces are to ones with two-dimensional nullspaces, as measured by the matrix 2-norm (Golub & Van Loan 1996). The distribution of all the singular values, in turn, contains important information regarding the spacing and density of the equilibrium states. Golub & Van Loan (1996, p. 73) have proved that the set of full rank matrices in ^{mxn} is both open and dense. Whether or not this is true of the set of configuration matrices is not clear, although figure 1 certainly provides numerical evidence that the set is dense. In addition, since the size of the smallest singular value can be thought of as the distance between a full rank matrix and the nearest singular one (Trefethen & Embree (2005) for interesting discussions of this topic related to robust control theory), the numerics also shows that the set of equilibrium configurations is not dense.

It is interesting to note that the typical distribution of singular values associated with the equilibrium configurations is definitely *not* what a random matrix would produce. As an example, we show in figure 5 the plots of the six singular values associated with matrices that are *N*(*N*−1)/2×*N* for *N*=6. The singular values are normalized by dividing by their sum, i.e. plotted are for two cases. The unfilled circles, connected by short dashed lines, are the averaged normalized singular values that a random matrix produces. We independently choose each of the entries of the matrix in this case using a random number generator, where each entry is a real number between 0 and 1, chosen from a uniform distribution. For this data, 10^{4} matrices are independently produced and their normalized singular values are averaged. By contrast, the filled circles, joined by longer dashed lines, are the averaged normalized singular values produced by actual configuration matrices whose entries are chosen to satisfy the relations (2.14). The vortex positions are chosen by a random number generator, but because the entries must satisfy (2.14) to be a configuration matrix, the entries are correlated. For this, the distribution has more of a Gaussian shape, decreasing from the peak much more slowly than that of the random matrix ensemble. Whether or not the actual distribution of singular values corresponding to equilibrium configurations, for large *N*, are, in fact, Gaussian distributed, or perhaps follow some other probability distribution, is at this point an open question. In fact, characterizing the distributions of classes of matrices, random or otherwise, is an active area of research (Sengupta & Mitra 1999).

## 5. Discussion

The method presented in this paper characterizes an equilibrium configuration as a point in a vector space whose non-normal configuration matrix has at least one zero singular value (i.e. is rank deficient) and whose strength vector lies in the spanning set of the right-singular vectors of *A* corresponding to singular values which are zero. This viewpoint complements the classical way of viewing each as an extremizer of an augmented Hamiltonian via Kelvin's variational principle. We believe this characterization adds new insights into the problem and allows us to conveniently classify equilibria according to the SVD of its configuration matrix as well as its Hamiltonian energy (2.2) and gives us a practical way of measuring quantities such as the size of each of the equilibria, distances between them and the distance between a generic non-equilibrium configuration and the nearest equilibrium configuration to it. The Brownian ratchet scheme, in practice, converges fairly quickly to configurations with non-trivial nullspaces if the nullspace is low dimensional, but from a practical point of view, for large *N* and for configurations with high-dimensional nullspaces, more efficient and stable random walk and SVD schemes will need to be implemented. It is a curious fact that, for the most part, the equilibrium states identified by the ratchet scheme have turned out to be asymmetric, despite the fact that there is nothing, in principle, that prevents the scheme from converging to symmetric states. One explanation for this could be that our initial configurations are generated randomly, hence typically have no inherent symmetries. It therefore seems unlikely for the scheme to converge towards an equilibrium configuration which is highly symmetric, unless there are no asymmetric states around. However, another more intriguing possibility might be that asymmetric states are far more prevalent than symmetric ones, despite the fact that much more is known about symmetric equilibria.

## Footnotes

↵The smallest eigenvalue can also be thought of as a 2-norm measure of the distance between this non-equilibrium configuration and the closest configuration which is an equilibrium for some choice of

.*Γ*- Received November 2, 2006.
- Accepted February 15, 2007.

- © 2007 The Royal Society

## References

## Notice of correction

Equations (2.12) and (2.14) are now presented in their correct forms.

A detailed erratum will appear at the end of the volume. 26 April 2007