## Abstract

We formulate the problem of finding equilibrium configurations of *N*-point vortices in the plane in terms of a gradient flow on the smallest singular value of a skew-symmetric matrix *M* whose nullspace structure determines the (real) strengths, rotational frequency and translational velocity of the configuration. A generic configuration gives rise to a matrix with empty nullspace, and hence is not a relative equilibrium for any choice of vortex strengths. We formulate the problem as a gradient flow in the space of square covariance matrices *M*^{T}*M*. The evolution equation for drives the configuration to one with a real nullspace, establishing the existence of an equilibrium for vortex strengths that are elements of the nullspace of the matrix. We formulate both the unconstrained gradient flow problem where the point vortex strengths are determined *a posteriori* by the nullspace of *M* and the constrained problem where the point vortex strengths are chosen *a priori* and one seeks configurations for which those strengths are elements of the nullspace.

## 1. Introduction

The problem of finding and classifying all fixed and relative equilibria of *N*-point vortices in the plane or on the sphere is of long-standing interest. In the plane, attacks on the problem date back to the 1800s with the works of Kirchhoff and von Helmholtz, who first formulated the problem as a Hamiltonian system of interacting particles—a complete bibliography of papers on the subject has been compiled recently by Meleshko & Aref (2007). When the particles are restricted to lie along a line, equilibrium positions are known to be related to the zeros of classical functions (Aref 2007). More generally, one finds deep connections with equilibria in the plane and rational solutions of KdV and the Adler–Moser polynomials—see Aref *et al.* (2003) and Aref (2007) for entries into this literature. On the sphere, the problem is related to Smale’s 7th problem (Smale 2000) of how to arrange *N* points on the surface so as to minimize the Rieszs-energy (Saff & Kuijlaars 1997) of the system in the limit . These points are generally known as Fekete points, or elliptic Fekete points when the interaction is logarithmic. Although much is known for special values and ranges of *N*, the general solution to this problem has not been found and is of interest, for example, to gain a better understanding of the ways in which complex virus molecules arrange themselves (Longuet-Higgins 2009), or how to distribute the nodes of a computational mesh on a general surface (Szeliski *et al.* 1993). It has emerged as a canonical ‘benchmark’ problem in constrained optimization theory that is routinely used as a challenge for numerical codes of all sorts, ranging from genetic algorithms, simulated annealing and gradient methods (Thomson 1904; Bergersen *et al.* 1994; Altschuler & Perez-Garrido 2005). The article of Aref *et al.* (2003) contains a relatively comprehensive survey of what was known at the time of writing, and the recent article of Newton & Chamoun (2009) highlights questions that can be addressed with the matrix-based approach used in this paper. The closely related problem of classifying relative equilibria for *N*-body gravitational interactions is also not yet solved and continues to present challenges (Hampton & Moeckel 2006). For a general introduction, see Newton (2001).

By formulating the problem as one of identifying fixed points of the system written in terms of interparticle distances, recent progress has been made in identifying and classifying equilibria of all types. As the system is linear in the particle strengths, the relative equilibrium problem can be reduced to a problem in linear algebra (Jamaloodeen & Newton 2006; Newton & Chamoun 2007, 2009). In order to set the stage for the gradient flow method detailed in this paper, we first give the necessary background.

The *N*-vortex problem in the plane can be written in terms of the *N* position vectors and real particle strengths *Γ*_{α}≠0 as
1.1
where
1.2
is a 2×2 skew-symmetric matrix and ∇_{α} is the two-dimensional partial gradient vector with respect to . The Hamiltonian ℋ is given by
1.3
with being the intervortical distances. The condition for the configuration to form a relative equilibrium for a given set of particle strengths is that all intervortical distances remain constant, i.e. *r*_{αβ}=const. for *α*≠*β*.

As the dynamics takes place in , one can formulate the problem as a dynamical system in the complex plane, with representing the position coordinates of the *α*th vortex. This *N*-dimensional complex dynamical system takes the form
1.4
with Hamiltonian
1.5
where we treat *z* and *z** as independent variables, adopting the notation
1.6and
1.7
The dynamical equations obtained from equations (1.4) and (1.6) are then
1.8
where the prime on the sum indicates that the term *β*=*α* is omitted in order to avoid the singular term. In addition to the Hamiltonian equation (1.5), the interacting particle system (1.8) has three other conserved quantities. The first,
1.9
is due to invariance of the system with respect to rotations in . The other two arise from translational invariance of the system
1.10
As long as the total vorticity of the system is not zero, it is traditional to normalize this by the total vorticity and to define the centre-of-vorticity as
1.11
Hence, a rigidly rotating configuration must rotate about the centre of vorticity *z*_{0}, which we can always choose as our origin *z*_{0}=0. For the time being, we will focus on configurations which rotate, so we assume that and thus *z*_{0}=0. Rotating patterns have been studied in Campbell & Ziff (1979), Gueron & Shafrir (1999), Aref & van Buren (2005) and Lewis & Ratiu (1996), for example. A more comprehensive treatment can be found in Aref *et al.* (2003) and Newton & Chamoun (2009). Later, we will consider the case of translation. Note that the above argument implies that a collection of vortices can only translate if
1.12
so configurations which translate uniformly should be of co-dimension one in the set of all configurations which are rigid.

The equations for rigid rotation about the origin are the following *N*+1 complex equations
1.13
where *ω* is the rotational frequency. This gives *N*+1 homogeneous equations in *N*+1 real unknowns *Γ*_{α},*ω*, with *N* complex parameters *z*_{1},*z*_{2},…,*z*_{N}. This same algebraic system can be obtained by extremizing the augmented Hamiltonian
1.14
where *u*+*i**v*=*V* and *ω* play the role of Lagrange multipliers and the constraints are the conservation of linear impulse and angular impulse . Thus, relative equilibria are obtained as critical points of the augmented Hamiltonian equation (1.14), a variational principal owing originally to Kelvin, as discussed in Aref *et al.* (2003) and used by Campbell & Ziff (1979) for computational purposes, and Gueron & Shafrir (1999) for analysis.

## 2. The linear system

As the system (1.13) is linear in the particle strengths and rotational frequency, one can write the equations as a linear system
2.1
where is an (*N*+1)×(*N*+1) matrix with complex entries and is a vector in . As such, we would like to characterize the nullspace structure of *M*, noting that *M* is skew-symmetric, i.e. *M*^{T}=−*M*.

### (a) Singular value spectrum

One may obtain the nullspace of a general matrix by performing a singular value decomposition and using the right singular vectors corresponding to the zero singular values as a basis set. In our case, we obtain the *N*+1 singular values *σ*_{i} (called the singular spectrum) of *M*, along with their corresponding left and right singular vectors , by solving the coupled linear system
2.2
where and *M*^{†} denotes the hermitian conjugate (adjoint) of *M*. The left and right singular vectors are used as columns to construct the unitary matrices *U* and *V*
2.3
which produces the singular value decomposition of *M*
2.4
is a diagonal matrix with (real) singular values down the diagonal, ordered from smallest (top left) to largest (bottom right)
2.5
The decomposition (2.4) expresses *M* as a linear superposition of the rank-one matrices , (*i*=1,…,*N*+1) with weighting determined by the singular values *σ*_{i}. The dimension of the nullspace of *M* corresponds to the number of zero singular values, while the right singular vectors corresponding to those zero singular values form a basis set for the nullspace of *M*. A comprehensive discussion of the many advantages and ramifications of adopting this approach is described in Newton & Chamoun (2009).

Note that *N**u**l**l*(*M*)≠∅ is a necessary but not sufficient condition for the existence of a vortex equilibrium. As while the vorticities are assumed to be real, one needs the nullspace of *M* to be real. It is easy to see that the condition that *M* have a real nullspace is equivalent to the condition that *M*_{R} and *M*_{I}, the real and imaginary parts of the matrix *M*=*M*_{R}+*i**M*_{I}, have nullspaces with a non-empty intersection. It is not hard to see that if and only if
2.6
The one direction is obvious. To see that implies the existence of a real non-zero vector which is in the nullspaces of *M*_{R} and *M*_{I}, note that, given , we have
For a generic collection of position coordinates *z*_{α}, both the matrices *M*_{R} and *M*_{I} will have full rank, hence an empty nullspace. However, if we view the smallest singular value associated with this matrix as a scalar which we would like to ‘drive’ to zero, we can formulate a gradient flow method based on the singular spectrum of *M*, which we describe next.

## 3. Gradient flow on the singular spectrum

### (a) A simple example

We first introduce some notation. Let **M**_{n,k} denote the space of *n*×*k* matrices, and be a function from *n*×*k* matrices to the reals. The gradient of *f* with respect to *M* is an *n*×*k* matrix with entries given by
3.1
For example, from the minor expansion, one has the well-known relation that the derivative with respect to the argument of the determinant is given by
3.2
where the last equality holds if . Recall that the cofactor matrix is the matrix whose (*i*,*j*) entry is (−1)^{i+j} times the (*i*,*j*) minor determinant of *M*. More generally, if *M* is a real non-square matrix, then the following identity holds:
3.3

To start, let us suppose that, given a matrix with empty nullspace, we wish to find a nearby matrix with a non-trivial nullspace. One tempting strategy is to employ a gradient flow, with a Lypanuov function given by . More specifically, consider the flow
3.4
where *s* denotes a generic ‘flow’ parameter. This flow has the following nice properties. Assuming that , it is easy to check that this flow leaves the singular value basis invariant, changing only the singular values, as
which is manifestly diagonal. As *M* and *M*_{s} share a singular value basis, the dynamics takes place entirely in the singular values. If *σ*_{i} are the singular values of *M*, then the evolution of *σ*_{i} is given by
3.5
This flow possesses a number of conserved quantities. It is clear that the difference of the squares of any two singular values is conserved, , and thus the flow is reducible to quadratures:
3.6
Also, it is easy to see that the flow for *σ*_{1} has an unstable fixed point at *σ*_{1}=*σ*_{2} and a stable fixed point at *σ*_{1}=0, so for generic initial data it will flow to the configuration where the smallest singular value is zero and the rest are non-zero. We note the similarities to some previous work on completely integrable gradient flows by Bloch *et al.* (1992).

### (b) The vortex system

We now specialize to the vortex system (2.1). Using the identities derived earlier, it is straightforward to find the gradient flow associated with this Lyapunov function. If we define , then using the chain rule, we have that
This gives the following expression for the gradient flow:
3.7
and similarly for the *y* coordinates. This system admits a clear physical interpretation. Originally, we considered matrices , but we are in fact only interested in matrices of the form (2.1), which forms an (complex) *n*-dimensional subset. As *t**r*(*M*^{†}*M*) defines an inner product on the space of square matrices, and form a basis for the tangent space to the set of matrices derived from a vortex configuration, the above flow simply represents the original gradient flow projected down onto the constraint set.

It is worth noting the following that for the unconstrained gradient flow, it is clear that the vector field vanishes only when (as the cofactor matrix of a non-zero matrix cannot be identically zero). For the above constrained flow, it is no longer clear that this is the case. It is possible that the gradient of the Lyapunov function is orthogonal to the tangent space of the constraint set (, and similarly *y*_{i}). Checking this analytically seems prohibitively difficult, and we have not done so. Nevertheless, this does not seem to be a problem numerically as we now show.

## 4. Numerical implementation

We now present some numerical results based on the above formulation. Figure 1 represents a direct numerical implementation of the gradient flow algorithm. The initial locations of the vortices were chosen randomly in the unit square. The vortices were then allowed to flow according to the constrained gradient flow (3.7). In the figure, the dot size representing the individual vortices are chosen to be proportional to the absolute magnitude of the vorticity *Γ*_{i}.

While the gradient descent algorithm worked well for relatively small configurations of vortices, of the order of 10 or so, it was found to be too slow for large ensembles of vortices. For this reason, we compute rigid equilibria by applying standard optimization algorithms to find local minima of the Lyapunov functions just described. We use a commercial implementation, the Matlab Optimization Toolbox, subroutine ‘fminunc’. The algorithms used by this software are suitable to a large-dimensional problem with a dense, unstructured Hessian; a subspace trust region interior reflective method (Branch *et al.* 1999). A basic outline of the methods is as follows.

—

*Step 1.*Determine the two-dimensional subspace such that and solves , where*H*(*f*) is the Hessian of*f*. That is, the subspace contains the gradient direction and the ‘Newton’ direction.—

*Step 2.*Take a trial step by approximately solving the minimization problem: That is, model the function locally by its Taylor series and seek a minimum of this approximation within a bounded region of the current point. Detailed algorithms to perform this predictive step can be found in Byrd*et al.*(1988) and Shultz*et al.*(1985).—

*Step 3.*If in fact , keep the new value and return to step 1. If not, reduce the size*C*of the ‘trust region’ and repeat step 2.

The implementation of this algorithm typically converges quickly for *N*≈50 vortices. Because this implementation is Matlab-scripted, a compiled version would have the potential for a significant (by a factor of 10–100) speedup. We have used this methodology to numerically compute rigid vortex configurations which either rotate, translate or are stationary. We have done this in both the cases where the vortex strengths are fixed (typically to ±1) or where, as above, they are determined as part of the nullspace of the underlying linear algebra problem (which we refer to as the case of free vortex strength).

### (a) Rotating solutions ( free vortex strengths)

Next, we applied the quasi-Newton root-finding algorithm to the function as outlined in the previous section. Using the quasi-Newton algorithm, we were able to handle a larger number of vortices than was possible with the gradient algorithm. Some configurations that were computed using this algorithm are shown in figure 2. Note that there is a pattern—while still quite irregular, the vortices seem to form a rough triangular lattice, with the vortices in the centre tending to have smaller vorticity and the vortices near the edges tending to have larger vorticity.

### (b) Rotating solutions (fixed vortex strengths)

One can also use similar ideas to compute rigidly rotating solutions where the circulations of the vortices are all the same. We can proceed in a similar way with a somewhat different choice of Lyapunov function. Having all circulations equal is equivalent to having the vector in the nullspace of *M*_{R} and *M*_{I}. The total vorticity *ω* is positive in this case and can, by the scaling invariance, be chosen to be 1. The simplest way to construct such solutions is to choose the Lyapunov function
The gradient flow can be written simply in terms of *M*_{R}, *M*_{I} and the derivatives of these matrices with respect to *x*_{j} and *y*_{j}. The results of one of these experiments are shown in figure 3. The pictures produced in this case are rather strikingly different from those produced in the case where the circulations are not fixed. In the case where the circulations are free, the equilibrium configurations tend to resemble a somewhat irregular triangular lattice, with the vortices near the centre tending to have smaller circulation and those near the boundary tending to have larger circulation. In the case where the circulations are fixed to have the same magnitude, however, the pictures look very different, and have a much more one-dimensional feel. The vortices tend to lie along curves in the plane with large fairly large ‘faces’ containing few or no vortices.

### (c) Translating solutions (free vortex strengths)

Another problem that one can consider is to find vortex configurations which translate uniformly. The equations governing uniform translation are slightly different. One needs that
Here, we have chosen the translation to occur in the *x* direction at unit speed. Obviously, this can be rotated/translated to occur in any direction and at any velocity.

The above system can be written in the form
4.1and
4.2
with and the real and imaginary parts of the matrix and the matrix of all ones. Note that the matrix is not quite the same matrix as we used to find uniformly rotating solutions, and is in fact *n*×*n* rather than (*n*+1)×(*n*+1). From the above equations, we need to require that have a nullspace, and that We seek solutions by minimizing the Lyapunov function with the associated gradient flow on vortex coordinates. Once we have a solution, we recover the circulations via .

### (d) Translating solutions (fixed vortex strengths)

We can also seek uniformly translating solutions with specified circulation strengths. If is the desired vector of circulation strengths, we can use the Lyapunov function 4.3 A necessary but not sufficient condition for translating solutions is that the circulation strengths sum to zero. Examples of solutions found by both approaches are shown in figure 4.

### (e) Stationary solutions (fixed vortex strengths)

Stationary solutions satisfy
that is, and (stated as in §4*c*) share a real nullspace. For the case of free vortex strengths, we can seek solutions by minimizing while for the case of fixed vortex strengths we can minimize where is a constant vector of vortex circulations. It is known that for vortices constrained to have the same absolute circulations (but allowing either sign), there is a strong constraint on the vortex circulations that admit a stationary solution. The number of positive vortices and negative vortices must be consecutive triangular numbers, *N*^{−}=(1/2)*n*(*n*−1) and *N*^{+}=(1/2)*n*(*n*+1), in order for a stationary solution to exist (Aref 2007 and references therein), and thus the total number of vortices must be a square.

In figure 5, we see several examples of stationary equilibria for constrained vortex circulations. Figure 5*a*,*b* show symmetric configurations in which the nine vortices are relatively equally spaced. Figure 5*c*,*d* shows configurations in which vortices are grouped in two scales: small equilibrium configurations are nested inside a large-scale vortex arrangement. Insets show the nested clusters which are reminiscent of the (1,3) Adler–Moser polynomial solutions.

As the solutions to the problem of stationary vortices with circulations of equal magnitude are known to be solved in terms of Adler–Moser polynomials, this makes an excellent test problem. After running the quasi-Newton algorithm, we then checked the polynomials numerically to verify that they were, to good approximation, Adler–Moser polynomials. Here is a brief synopsis of the calculation for figure 5*c*, with nine vortices (six of the majority species, three of the minority species).

The locations of the vortices of the majority species (six vortices) are given by the roots of the sixth-order polynomial
while the locations of the vortices of the minority species are given by the roots of the cubic
The general Adler–Moser polynomials of degrees three and six, respectively, are given by
4.4and
4.5
Note that frequently in the literature, the translation invariance is removed by choosing *z*_{0}=0. We take *z*_{0}=(5.7221−1.9185*i*)/6 which exactly removes the fifth-order term from the sextic polynomial. We can read off the coefficients *τ*_{2},*τ*_{3} from the cubic and linear terms of the polynomial, giving
4.6and
4.7
Taking the difference between the Adler–Moser polynomial with the above coefficients and the polynomial found numerically gives
4.8and
4.9
Note that no effort has been made to make any kind of best fit—we simply read off the coefficients. Nevertheless, the residual is quite small, especially when compared with the coefficients of the original polynomial. We have done similar calculations with the other configurations—to a good approximation, they appear to be Adler–Moser polynomials.

## 5. Discussion

The methods described in this paper are based on a formulation of gradient flow for the singular spectrum, i.e. equation (3.4), which for the vortex equations become the more complicated system (3.7). The power of the method is that it is capable of finding all equilibria associated with equation (1.8), irrespective of their stability properties. This has advantages as well as disadvantages over gradient methods based on the augmented Hamiltonian equation (1.14). It is clear from our numerical explorations that as *N* increases, so do the number of equilibria (typically, this increase is cited in the literature as exponential); hence a categorization of the equilibria based on the full singular spectrum of the configuration matrix, as detailed in Newton & Chamoun (2009), becomes valuable. In the future, it would, of course, be interesting to distinguish between those that are stable and those that are not. This could, in principle, be achieved by combining the spectral flow method described in this paper with an algorithm that simultaneously minimizes the augmented Hamiltonian equation (1.14). An open question in the field is whether there exist any stable asymmetric equilibria (Aref & Vainchtein 1998; Cabral & Schmidt 1999).

## Footnotes

- Received August 11, 2009.
- Accepted November 23, 2009.

- © 2009 The Royal Society