The motion of point vortices on closed surfaces

We develop a mathematical framework for the dynamics of a set of point vortices on a class of differentiable surfaces conformal to the unit sphere. When the sum of the vortex circulations is non-zero, a compensating uniform vorticity field is required to satisfy the Gauss condition (that the integral of the Laplace–Beltrami operator must vanish). On variable Gaussian curvature surfaces, this results in self-induced vortex motion, a feature entirely absent on the plane, the sphere or the hyperboloid. We derive explicit equations of motion for vortices on surfaces of revolution and compute their solutions for a variety of surfaces. We also apply these equations to study the linear stability of a ring of vortices on any surface of revolution. On an ellipsoid of revolution, as few as two vortices can be unstable on oblate surfaces or sufficiently prolate ones. This extends known results for the plane, where seven vortices are marginally unstable (Thomson 1883 A treatise on the motion of vortex rings, pp. 94–108; Dritschel 1985 J. Fluid Mech. 157, 95–134 (doi:10.1017/S0022112088003088)), and the sphere, where four vortices may be unstable if sufficiently close to the equator (Polvani & Dritschel 1993 J. Fluid Mech. 255, 35–64 (doi:10.1017/S0022112093002381)).


Introduction
The original development of the theory of vortex motion goes back to Helmholtz [1] and Kelvin [2], who formulated various theorems concerning vorticity, in particular the conservation of circulation in a 'perfect' fluid (e.g. inviscid, incompressible and subject to forces derived from a single-valued potential) [3,4] [5] first derived the equations governed by point vortices in the two-dimensional (2D) plane. Such vortices have finite, constant circulation but singular vorticity restricted to a point. Each vortex induces a singular velocity field, entirely in the azimuthal component (with no radial component owing to zero horizontal divergence), and as a result there is no self-induced motion. Instead, vortices move in the flow field generated by all others.
The resulting model-a nonlinear Hamiltonian dynamical system-has since been studied extensively and been shown to exhibit surprisingly rich dynamics, from steady motions to chaotic motions, strongly depending on the number of vortices (for a review, see [6][7][8]). The model also exhibits special statistical properties for large numbers of mixed-sign vortices, noted by Onsager [9], favouring the clustering of like-signed vortices at high energies.
Interest in 2D flow on a sphere has primarily stemmed from pioneering studies of weather forecasting around 1950 [10]. The so-called 'barotropic model' employed was founded on the conservative transport of the scalar 2D vorticity (the vorticity component normal to the surface) by a 2D non-divergent (incompressible) 'geostrophic' flow. This model has remained popular ever since in idealized studies of atmospheric dynamics, though it is too simplistic to be used predictively. In fact, the equations for 2D flow on the surface of a sphere go back much earlier, to Zermelo [11]. A Lagrangian formulation of the fluid particle motion in Cartesian coordinates x with origin at the centre of the sphere was later given by Dritschel [12]. These coordinates make use of the natural isotropy of the surface, and the constraint, |x| = constant, is automatically satisfied by the equations of motion (this also avoids apparent polar singularities in traditional latitude-longitude coordinates [13]). These equations were then applied to study point vortex dynamics on a (non-rotating) sphere by Polvani & Dritschel [14], who also extended Thomson's [15] stability analysis of a ring of n equal vortices [16] finding that as few as four vortices may be unstable if located within 35.264 . . . degrees of the equator. 1 Boatto & Cabral [17] further extended these analyses to prove nonlinear stability, while Boatto & Simó [18,19] included the stabilizing/destabilizing effect of additional polar vortices. In particular, it was found that, for suitably chosen circulations of the polar vortices, a latitudinal ring of n vortices is stable everywhere, for any latitude and for all n.
The motion of point vortices on more general surfaces has recently attracted mathematical interest, as the more general setting allows one to better understand the role played by the surface geometry. Kimura [20] analysed constant-curvature surfaces (the sphere and the hyperboloid, respectively, with positive and negative curvature). Boatto [21] extended Kimura's analysis to prove that positive and negative curvatures have, respectively, a destabilizing and a stabilizing effect on the linear and nonlinear stability of a ring of identical vortices. More recently, Kim [22] appeared to have further extended the analysis to surfaces of variable curvature, specifically the ellipsoid of revolution, adapting the analysis of the sphere carried out by Boatto and Simó in a preprint [18]. A general theory for compact surfaces of variable curvature, however, appears to be lacking. Significant advances have been put forth by Boatto & Koiller [23,24], who clarified the role of the surface geometry on the form of the Green function and the non-singular residual Green function called the 'Robin function'. In the absence of boundaries, the Robin function plays a role only on surfaces of variable curvature.
In this article, we develop a general theory for point vortex dynamics on a general compact, differentiable surface. This makes direct use of the Hamiltonian formalism, for which we state the general form of the vortex interaction energy for a system of n vortices. We apply the theory to surfaces of revolution and derive the equations of motion explicitly. We illustrate aspects of the vortex dynamics for several surfaces, and revisit the stability of a ring of vortices. We correct the previous analysis of Kim [22], who omitted the contribution due to variable surface curvature and used the incorrect area form (i.e. the incorrect Hamiltonian coordinates). This results in major differences in our conclusions.
The structure of the paper is as follows. In the next section, we derive the interaction energy for a system of n vortices on a closed surface. This is the Hamiltonian for the dynamical system. In §3, we show how to compute the Green function by conformal transformation, and how to properly account for the compensating vorticity field forced by the Gauss condition. In §4, we derive the equations of motion for surfaces of revolution; in §5, we illustrate the forms of various surface functions arising in the theory; and in §6, we analyse the linear stability and nonlinear evolution of a latitudinal ring of n identical vortices. A few conclusions are offered in §7.

Energy of a point vortex system on a compact surface (a) Incompressible flow induced by vorticity
We consider a general differentiable compact surface S. Such a surface can be viewed as embedded in R 3 , and we adopt the standard notation of r = (x 1 , x 2 , x 3 ) for the position of a particle on S. Nevertheless, only two parameters s = (s 1 , s 2 ) are needed to specify a particle on a 2D surface, i.e. r(s). 2 Suppose that a vorticity distribution ω(s, t) is given on S, 3 and that the corresponding flow u(s, t) is incompressible, i.e. ∇ · u = 0. Given that ω is conserved following a fluid particle and dr/dt = u, we seek to determine the motion of an arbitrary fluid particle induced by ω(s, t). Note that this motion generally redistributes ω(s, t) so the problem is fundamentally nonlinear.
As noted by Kirchhoff [5], for planar 2D incompressible flows, we can express the velocity field in terms of a suitable differentiable function ψ, called the 'streamfunction', via where n is the local unit normal to the surface, viewed as embedded in R 3 . Upon substituting this into the definition of the (scalar) vorticity, ω = n · ∇ × u, the streamfunction ψ is seen to satisfy Poisson's equation, where is the Laplace-Beltrami operator (Laplace's operator restricted to S). Then, given ω, the solution of this equation provides ψ, from which we obtain u from (2.1) and hence the particle motion and evolution of ω.

(b) Inversion of the Laplace-Beltrami operator
Central to the analysis is the solution of (2.2) for ψ: the 'inversion' problem. To begin, note that for a compact surface S an immediate consequence of the divergence theorem is that the vorticity ω must integrate to zero over S, i.e.
This is the so-called 'Gauss condition'. In the integral, dΩ s is the differential area element expressed in the s variables. Owing to the linearity of Poisson's equation (2.2), the general solution, to within a trivial constant, is where G(s, s o ) is the fundamental solution or Green function. This can be obtained from

(c) Point vortices
Consider henceforth the special vorticity field ω associated with n point vortices, of circulations Γ j and at positions s j (t), for j = 1, . . . , n. By linear superposition, we have where, for notational convenience, we define Then from (2.4), the streamfunction ψ(s, t) at any point s = s k , ∀k is given by Unfortunately, this cannot be used to derive the motion of the point vortices, since ψ is singular at each vortex and u is not defined. Instead, we use Hamilton's equations directly. The first step is to derive the interaction energy, the Hamiltonian H, for a system of n vortices on S. This energy is entirely kinetic, but, as we shall see, excludes a singular part which results in no motion.
The kinetic energy E is given by after integrating by parts. For point vortices, this reduces to where Γ = j Γ j will be used henceforth to denote the sum of the vortex circulations. In fluid mechanics, there is distinction to be made between a vortex, which generates the fluid dynamics, and a passive particle, which corresponds to an element of fluid that is passively advected by the velocity field produced by the vortices. For a passive particle, the equations of motion can be found by taking E to be the Hamiltonian H. The vortex motion, however, cannot be derived this way, since (2.8) is singular at each vortex position. In any case, the motion of a passive particle depends on the vortex motion, which we turn to next.

(d) Point vortex dynamics and the Hamiltonian
We make the standard assumption, proposed by Kirchhoff [5], that a vortex behaves as a particle in the velocity field of the other vortices. However, special care must be taken to remove the selfinteraction term that would result in singular behaviour. On the plane or the sphere, symmetry rspa.royalsocietypublishing.org Proc. R. Soc.
arguments alone lead to the conclusion that there can be no self-induced vortex motion, so that one merely has to exclude j = k from the sum in (2.8) when evaluating ψ at s = s j . This renders the energy E-now the 'excess energy'-finite, and we can take H = E for the Hamiltonian. Essentially, the circular streamlines around each vortex correspond to a flow with no radial component. For a general compact surface S, this motivates the following definition.
Definition 2.1. The streamfunction at the jth vortex is defined by the limiting process where d(s, s j ) is the geodesic distance between s and s j .
Remark. In the limit above, the (locally smooth) surface is equivalent to a portion of a planar surface, and on such a surface the streamfunction has circular streamlines. As a consequence there is no singular self-induced motion owing to the vorticity singularity.
Hence, using the expression for ψ in (2.8), it follows that is known as the Robin function [25][26][27]. 4 We can now state the form of the finite, excess energy.

Proposition 2.2. The excess energy of a system of n point vortices, equivalently the Hamiltonian H, is
14) The proof follows immediately from (2.10) using (2.12).

Remarks.
(a) This expression holds for any closed, differentiable, genus zero surface (i.e. any surface topologically equivalent to a sphere). (b) The last term does not contribute to the dynamics as the integral of ψ over the whole surface is a constant (see remarks after proposition 3.1 in §3 below). This enables one to simplify the vortex Hamiltonian to as in [23,24], where the Green function part describes the interaction between pairs of distinct vortices, while the Robin function part can be viewed as the Hamiltonian describing the interaction of a single vortex with its uniform compensating vorticity spread across the surface. It is through R that a single vortex can still move on S. Explicit forms for R are given in §5. See also appendix B. (c) Lin [25] derives a directly analogous expression for H, called the 'Kirchhoff-Routh' function (see (4.4) in that paper), for describing vortex motion in planar domains.

Calculation of the Green function by conformal transformation
We now consider a surface of revolution S (about the vertical z axis), for which much analytical progress can be made (a discussion of more general surfaces is deferred to §3d). The Cartesian coordinates x in R 3 of any point on S may be expressed as functions of two surface coordinates θ and φ, co-latitude and longitude, respectively (note: s = (θ, φ)). For surfaces of revolution, it is sufficient to take where ρ(θ ) and ζ (θ) are specified functions of θ . Without loss of generality, we may take 0 ≤ φ ≤ 2π and 0 ≤ θ ≤ π over S. Note, on a spherical surface ρ = sin θ and ζ = cos θ.

(a) Conformal transformations for surfaces of revolution
We first make a conformal transformation from the punctured surface S (i.e. the surface without a point, S p ) 5 to the plane R 2 . The first step is to compute the differential distance ds between two points on S, where primes denote differentiation w.r.t. θ. We then consider the map Φ : , the usual planar polar coordinates. The goal is to find the function r(θ) in terms of which ds 2 can be written as where λ 2 (θ ) is a conformal factor (distances are not preserved but angles are under a conformal transformation). By equating the above two expressions for ds 2 in (3.2) and (3.3), we find equations defining both the conformal factor λ 2 and the corresponding map r = r(θ): This implies Starting from the natural boundary condition r(0) = 0, in principle these equations can be solved given ρ(θ ) and ζ (θ ). Some examples are given in §5.

Remarks.
(a) Note that r → ∞ when θ → π . That is, the point θ = π is mapped to infinity, since the conformal transformation applies only to the open surface θ < π. (b) For surfaces conformal to the sphere, as considered here, we can alternatively map directly from the sphere [23,24]. The planar map however is simpler. is given by the product of distances in the perpendicular coordinate directions, i.e.
This may also be written −dμ ∧ dφ where μ(θ) is determined from The coordinates μ and φ are said to be the 'natural' coordinates-also called the Darboux coordinates-for the surface S. For a more formal approach, see Do Carmo [29].
(b) The Green function for the punctured surface S p We have just shown that every punctured surface can be conformally mapped into the plane R 2 . That is, there exists a conformal transformation Φ : S p → R 2 , with conformal factor λ 2 (s), which maps S p into R 2 . On the plane, the vortex motion depends simply on the planar Green function, where in this subsection x and x o denote 2D planar coordinates. G P is the solution of Poisson's equation for a very special source term, a Dirac delta distribution, Here a subscript x is attached to to distinguish the Laplace operator in R 2 from the Laplace-Beltrami operator on S. We continue to denote the latter by throughout.
To relate G P to the Green function G S p for the punctured surface S p , it proves convenient to write G P explicitly in polar coordinates (r, φ), As shown in appendix C (and generally taken for granted), the conformal transformation Φ : S p → R 2 maps (3.8), the equation satisfied by G P , into We can therefore conclude that the Green function for the punctured surface is and moreover G S p satisfies with a bare Dirac distribution as its source.
Remark. This Green function does not apply to the closed surface S, since G S p does not enforce zero mean vorticity, as pointed out in §2, see (2.3). That is, the integral of G S p over S does not vanish (in fact, it is equal to 1).
(c) Extending the Green function to the whole surface S We next show how to recover the Green function G S for the whole closed surface S, once the Green function G S p for the punctured surface S p is known.
Here, recall dΩ is the total area of the surface.
Proof. Applying the Laplace-Beltrami operator to G S , we obtain which is exactly the condition we require for the Green function on S. As verified in appendix A, the extra terms added to G S p in (3.13) cancel the singularities in G S p at θ = π and θ o = π , so that G S is regular at these points. That is, G S is finite at θ = π and θ o = π , justifying taking the domain of integration in (3.16) as the whole compact surface S. Hence, G S is the Green function for S, to within an unimportant constant.
Remarks. Properties of the Green function.
is just a constant. In the second line, we have used (3.14) to evaluate the integral over G S p . Strictly speaking,Ḡ is found by integrating over the punctured surface S p , but the additional point in S contributes nothing to the integral. A simple proof is given in appendix A.
Remark. For a spherical surface, discussed in §5a below, one can show that this construction gives the correct form for G S , proportional to the chord distance between s and s o .

(d) General surfaces
For more general surfaces (of genus zero) without rotational symmetry such as a tri-axial ellipsoid, we can in principle construct the Green function in the same way as outlined for surfaces of revolution. The procedure is as follows: (1) find a conformal transformation between the surface S and the plane R 2 ; (2) make use of the known Green function G S p for the punctured surface; and (3) correct for the condition of zero total vorticity, i.e. set G S (s, The recipe is clear, though the devil is in the detail!

(e) The Hamiltonian
The Hamiltonian H derived in §2 simplifies when we make use of the Green function derived just above. To simplify the notation, from now on we shall suppress the subscript S and simply use G(s, s o ) for the Green function on the surface S.
As previously remarked, the last term in the Hamiltonian H in (2.14) is just a constant, independent of the vortex positions, so does not contribute to the vortex dynamics. (This is explicitly shown using the expression for G above in the remarks after proposition 3.1, in (3.17).) It is sufficient therefore to consider the simplified Hamiltonian given in (2.15).
For a surface of revolution, the Robin function R(s j ) defined in (2.13) also simplifies considerably. In the expression for R, in the limit s → s j , the geodesic distance d(s, s j ) has the same form as the infinitesimal distance, ds, introduced at the beginning of §3a. Explicitly, where r = r(θ) and r j = r(θ j ), and the second identical expression is introduced to draw a parallel to G S p : using (3.10) and (3.11), with x = (r cos φ, r sin φ), x j = (r j cos φ j , r j sin φ j ), as before. Now it is evident that the log singularities in R cancel, leaving where λ(s j ) = ρ(θ j )/r(θ j ) depends only on θ j for a surface of revolution. Furthermore, as shown in (A 5),Ḡ(s j ) depends only on θ j , so that the Robin function R also depends only on θ j .

Remark.
For a bounded 2D domain, the Robin function R has the exact same dependence on the conformal factor λ arising in the conformal transformation mapping the boundary to the unit circle (e.g. §3.3 in [8] and also appendix C). However, there is no analogue of the functionḠ, which arises here because the surface S is closed, as shown in proposition 3.1.

Equations of motion
Given the Hamiltonian H(s 1 , s 2 , . . . s n ), the equations of motion may be generally written where a dot indicates a time derivative, and is the standard skew-symmetric 2 × 2 matrix.

(a) Surfaces of revolution
Surfaces of revolution have a natural 'position-momentum' canonical pair (q k , Here μ and φ are the 'natural' coordinates for which the area element dΩ s = −dμ dφ (see remark (c) at the end of §3a). Moreover, incompressibility in these coordinates reduces to (in planar polar coordinates, for example, μ = r 2 /2; see [30]). Hence, there exists a streamfunction ψ such thatφ For the point vortex dynamics, the very same equations apply for each vortex k, and the streamfunction is just that given in (2.12), i.e. the regular part of the streamfunction evaluated at the vortex position. Equivalently, and consistently,φ k andμ k can be obtained from Hamilton's equations (4.1), identifying φ k as the 'position' and Γ k μ k as the 'momentum': In fact, the identification of Γ k μ k with 'momentum' is appropriate since it is immediate to show that is conserved owing to the rotational symmetry of H (see below). This is just the angular impulse. We defer the explicit calculation of the partial derivatives in (4.6) and first introduce a transformation of the equations of motion that is beneficial for their numerical implementation. We begin by rewriting the equations of motion (4.6) in the coordinates (φ, θ), for whicḣ where μ = −ρ (ρ ) 2 + (ζ ) 2 ; see (3.6). Note that μ (θ) < 0 everywhere except at the poles θ = 0 and π where μ = 0. There,φ k is singular because arbitrarily rapid angular variations occur for vortex trajectories crossing through either pole (θ k , however, is finite). The singularity inφ k can be avoided using instead three coordinates, analogous to Cartesian coordinates, constrained to the surface [12], X = sin θ cos φ, Y = sin θ sin φ and Z = cos θ, (4.9) satisfying X 2 + Y 2 + Z 2 = 1. Their time derivatives are given bẏ where u k = sin θ kφk is the 'azimuthal' velocity and v k = −θ k is the 'meridional' velocity-both of which are finite at all points, including the poles. While these equations are redundant, they are regular everywhere and convenient for numerical implementation. We next turn to the calculation of the partial derivatives ∂H/∂φ k and ∂H/∂θ k in (4.6). To this end, note that, for a surface of revolution, where r j = ρ(θ j )/λ(θ j ) and r k = ρ(θ k )/λ(θ k ), while (4.12) To simplify notation, we suppress the subscript k and change j to o. Functions evaluated at θ are subscript free, while those evaluated at θ o acquire a subscript o.
It proves useful to eliminateḠ using the definition of R above in (4.12): Then, after using r = ρ/λ and r o = ρ o /λ o (here, as before, α = φ − φ o ). We next introduce the function ξ = ρ 2 λ , (4.14) which is finite for all θ . Then, since R depends only on θ, the partial derivative of G w.r.t. φ is given by The partial derivative of G w.r.t. θ is Here again the numerator and denominator involve only finite functions of θ and θ o . The derivatives ξ and λ are given explicitly by and, notably, they vanish at θ = 0 and π and are finite elsewhere. In ∂G/∂θ , only the derivative of the Robin function R is needed, and after some manipulation we find where we have taken μ(0) = A/(4π ) without loss of generality (then μ(π ) = −μ(0)). Note that R , like ξ and λ , vanishes at θ = 0 and π . Putting everything together (and noting that distinct pairs of integers appear twice in the double sum involved in H; see (2.15)), we can now write down the explicit form of the azimuthal and meridional velocity components, restoring the j and k subscripts: where we have introduced three additional functions and which are finite for all θ, and moreover β(0) = β(π ) = 0.

Examples
(a) The unit sphere S 2 In (3.1), the unit sphere, see figure 1a, is described by the functions ρ = sin θ and ζ = cos θ. The planar radius r(θ ) satisfies which can be integrated starting from r(0) = 0 to give in terms of which the conformal factor is given by 3) The meridional coordinate satisfies  Next, we turn to the calculation ofḠ(θ ). From §3a, this is given by which is singular only as θ → π , as anticipated from the discussion in §3a. The Green function, therefore, takes the form where |x − x o | is the chord distance (usual distance in R 3 ) between x and x o . This is the classical result [12]. The Robin function given by is just a constant, and therefore plays no role in the dynamics. Finally, the equations of motion can be expressed simply in terms of the Cartesian coordinates x k of each vortex (one needs only to derive the equation forż k and use symmetry to obtain the others forẋ k andẏ k ; see [12]):ẋ The forms of the various surface functions such as μ(θ) and λ(θ) are markedly more complex than in the case of the sphere. Defining q = τ −1 = (ρ ) 2 + (ζ ) 2 = cos 2 θ + b 2 sin 2 θ, the meridional coordinate μ takes the form (5.10) Note that the surface area A = 4πμ(0). The conformal factor λ is more complicated still, and we find where we have taken λ(0) = 2 without loss of generality (any value could be taken as only ratios of λ appear in the dynamical equations (4.19)). These forms of μ and λ, for both b < 1 and b > 1, have been verified by direct numerical integration to machine precision accuracy. Figure 2 shows the form of μ(θ) and λ(θ ), as well as ν(θ), the rotation rate of a single vortex of circulation Γ = 2π about the axis of symmetry. The rotation rate is given in terms of β(θ) from (4.21), specifically The meridional coordinate μ(θ) is antisymmetric about the equator θ = π/2. In general, μ varies more strongly with θ as b increases (as the shape becomes more prolate), which is to be expected since 2π (μ(0) − μ) measures the surface area from the north pole θ = 0 to any 'co-latitude' θ. Taller shapes have greater surface area than shorter ones. The conformal factor λ decreases to zero more rapidly for prolate shapes than for oblate ones. As b → 0, λ develops an extensive flat region where λ ≈ 2. The rotation rate of a single vortex ν, like μ, is antisymmetric about the equator θ = π/2. Positive-signed vortices on prolate shapes b > 1 rotate cyclonically (counter-clockwise) about the north pole if they are located in the northern hemisphere. In the southern hemisphere, they rotate the opposite way. The fastest rotation rates occur at the poles. The situation is reversed for oblate shapes b < 1, while, for a sphere b = 0, a single vortex is stationary. Notably, the motion of a vortex dipole (with Γ = Γ 1 + Γ 2 = 0) appears to be quasi-periodic in general. The average position of the dipole does not generally form closed orbits, unlike in the case of a spherical surface on which all orbits are great circles.
(5.13) When a = 0 and b = 1, we recover the sphere, while for general b the surface is a spheroid (an ellipsoid with semi-axis lengths 1, 1 and b). When a > 0, the upper part of the surface is dented inwards, producing a region of negative curvature-a 'bean' shape. An example is provided in figure 1c for the parameters a = 0.6 and b = 0.4. By 'curvature', we mean the Gaussian curvature K = κ 1 κ 2 , the product of the largest and smallest normal curvatures of curves tangent to S at a specific point. The normal curvature is measured in the plane spanned by the tangent vector and the normal to the surface. The Gaussian curvature can be computed in terms of the first and second fundamental forms, and, for a surface of revolution, reduces to For the bean-shaped surface, we find showing that-if 2a > b-then K = 0 at the 'rim' θ = θ r = cos −1 (b/2a) (where the surface reaches its maximum height a + b 2 /(4a)) and also at the smaller angle θ = θ m = cos −1 ( 3 b/2a). For θ m < θ < θ r , the Gaussian curvature is negative.
Remark. For a = 0, it does not appear possible to obtain closed-form expressions for the various functions λ, μ, etc. Nevertheless, they can be obtained numerically to high precision. Special care is needed, however, to find the conformal factor λ(θ ). From (4.17), we have The r.h.s. vanishes as θ → 0 but is singular for θ → π . To obtain an accurate expression for λ, it is necessary to add and subtract (cos θ − 1)/ sin θ from the r.h.s., as this term can be integrated exactly and compensates the singularity at θ → π : Note that λ(0) = 2 while λ(π ) = 0. The integral above and that to find μ(θ) from (where we are free to choose C so that μ(π ) = −μ(0)) are obtained using two-point Gaussian quadrature with 1800 equal divisions of θ between 0 and π , giving results accurate to one part in 10 14 for moderate values of a and b. Once these functions are computed, all other functions can be found directly, without integration.   Figure 3 shows μ(θ ), λ(θ) and ν(θ) for various values of a and for b = 0.5. For a > 0, the meridional coordinate μ(θ) is no longer symmetric about the 'equator' θ = π/2. An additional inflection occurs for a > 0.25; this inflection arises from a region of negative Gaussian curvature K < 0 located in a belt around the north pole; in general, this occurs when 2a > b (see discussion above). This inflection is also seen in the conformal factor λ(θ), which develops a short plateau before steeply dropping off for θ > π/2. The rotation rate ν(θ) like μ(θ) loses symmetry about the equator when a > 0. Two zones of strong anti-cyclonic rotation develop for the larger values of a. The first zone is centred on the rim at θ = cos −1 (b/2a), where the Gaussian curvature K changes sign. The second zone is centred on the south pole and is evidently associated with the increasing curvature there as a increases-note the similarity to the b = 2 case in figure 2. Essentially, the region near the south pole is becoming locally prolate, inducing anti-cyclonic rotation about the axis of symmetry (or cyclonic rotation about the south pole, considering this pole as 'up') (figure 1c).

Stability of a ring of vortices
To explore the effects of variable curvature on the dynamics of point vortices, we examine next the stability of a ring of vortices located at 'co-latitude' θ on a general surface of revolution (figure 4). We first consider linear stability, then illustrate several instabilities in the fully nonlinear dynamics. Specific results are presented for the ellipsoid and compared with the known results for a sphere [14].

(a) Linear stability
A system of n identical vortices (Γ k = Γ 0 , Γ = nΓ 0 ) initially lying on the co-latitude ring θ k = θ at equally spaced longitudes φ k = 2π k/n, k = 1, 2, . . . , n rotates steadily about the north pole at a fixed angular velocity ν (it is said to be a 'relative equilibrium'). The vortex motion is simply given by φ k (t) = 2π k/n + νt, θ k (t) = θ . The angular velocity ν =φ k can be calculated from u k / sin θ k (for any k) using (4.19). Taking Γ 0 = 2π without loss of generality, we obtain Above, all surface functions (β, τ , . . .) are evaluated at θ. However, ξ = ρ 2 /λ so λξ = ρ 2 . Also, λξ + λ ξ = (λξ ) = 2ρρ , so the sum reduces simply to 2(n − 1)ρ /ρ. Therefore, Using τ = 1/ (ρ ) 2 + (ζ ) 2 and the definition of β in (4.21), we obtain the final expression When n = 1, this reduces to the single-vortex rotation rate in (5.12). For a spherical surface, ν = (n − 1) cos θ/2 sin 2 θ , in agreement with that found in [14,18]. Linear stability is determined by displacing the vortices infinitesimally, taking φ k = 2π k/n + δφ k and θ k = θ + δθ k , and solving the resultant linearized equations for δφ k (t) and δθ k (t). Seeking eigen-solutions of the form δφ k (t) = {φ k e σ t } and δθ k (t) = {θ k e σ t }, the objective is to find the growth rate σ (or frequency if σ is purely imaginary) by solving the resultant linear algebraic equations. As in [14], this procedure can be simplified considerably by exploiting the rotational symmetry of the basic state. This symmetry is apparent in the general form of the matrix eigen-problem, (6.4) where v = (φ 1 , . . . ,φ n ,θ 1 , . . . ,θ n ) T and A, B, C and D are cyclic n × n matrices of the form a 1 a 2 a 3 · · · a n−1 a n a n a 1 a 2 · · · a n−2 a n−1 a n−1 a n a 1 · · · a n−3 a n−2 · · · · · · · · · · · · · · · · · · a 2 a 3 a 4 · · · a n a 1 As a result, the eigenvectors v correspond to special vortex displacements of the form φ k =φ e 2πikm/n andθ k =θ e 2πikm/n (m = 1, 2, . . . n). (6.6) Note, there are m distinct modes of displacement, and any general displacement can be represented as a linear combination of these modes. Substituting these forms forφ k andθ k into the matrix eigen-problem, we obtain a 2 × 2 problem for each mode m: The calculation of the a j , b j , c j and d j is straightforward but laborious, so details are omitted. It turns out that all a j = d j = 0, while, for the other coefficients, we obtain , j > 1 (6.9) and , j > 1, (6.10) where in b 1 and K is the Gaussian curvature defined in (5.14). Then, using the coefficientsb m andc m for each mode m simplify tõ Putting this altogether, the growth rate of each mode m is just σ = σ m = ± b mcm , or, using (6.11) and (6.13), where (m = 1, 2, . . . , n). In this expression, we have taken all vortex circulations Γ k to be 2π . If we instead take Γ k = Γ 0 , the above expression must be multiplied by Γ 0 /2π . Note that the rotation rate ν, defined in (6.3) and appearing in (6.14), is dimensionless. Like σ m the dimensional rotation rate is given by ν times Γ 0 /2π . For the sphere, this reduces to the form derived in [14], namely 2 sin 2 θ . (6.15) Instability occurs if m(n − m) > (n − 1)(1 + cos 2 θ ) or, equivalently, if This never occurs for m = 1, n − 1 or n. Hence, ring configurations of three or fewer vortices are stable on a sphere. For n = 4, instability occurs for cos 2 θ < 1/3, i.e. in an equatorial belt extending between θ = 54.7 . . . to 125.2 . . . degrees, and, as n increases, the instability region widens. For n ≥ 7, instability occurs at all θ. We now examine the ellipsoid. This surface is described by the functions ρ = sin θ and ζ = b cos θ (the sphere corresponds to the special case b = 1). Figure 5 shows the domain of instability as a function of b for n = 2-6 vortices (a single vortex is always stable while seven or more vortices are always unstable). Configurations with θ n (b) < θ < π − θ n (b) are unstable. Notably, as few as two vortices can be unstable, if sufficiently close to the equator, when b < 1 (oblate surfaces) or when b > 1.684541781 approximately (prolate surfaces). Three vortices always have a domain of instability near the equator except on the sphere b = 1. Four to six vortices are unstable sufficiently near the equator for all b > 0.
As b → 0, the domain of instability appears to shrink to zero. However, dθ does not represent distance on an ellipsoid. An alternative is to use arc-length s measured from the north pole to the co-latitudes θ = θ n (b), normalized by the arc-length to the equator s e (note: arc-length is equal to , where E is the incomplete elliptic integral of the second kind). Figure 6 Figure 6. Same as figure 5 except that the arc-length ratio s(θ n (b))/s e is shown. Here, s e is the arc-length from the pole to the equator. The instability domain lies above the curves shown for each n, wherever s/s e < 1. Note that s = 0 corresponds to the north pole. ratio s/s e as a function of b for n = 2-6 vortices. Comparing with figure 5, the domain of instability is seen to be larger for b > 1 and smaller for b < 1. The domain of instability is still seen to shrink to zero as b → 0. Overall, this implies that vortices on oblate surfaces are generally more stable than vortices on prolate surfaces.
These stability results differ markedly from those reported in [22]. In that work, the effect of the Robin function-arising from the compensating uniform vorticity field associated with the Gauss condition (2.3)-was overlooked. The Robin function plays a key role in the dynamics when the sum of the vortex circulations is non-zero, as here.

(b) Nonlinear evolution
We illustrate next several instabilities using the full nonlinear dynamics in (4.19). First, consider two equal vortices (Γ 1 = Γ 2 = 2π ) lying on the equator of a slightly oblate ellipsoid, b = 0.99. According to the linear analysis, this state is unstable ( figure 5). As an initial perturbation, we displace the vortex co-latitudes by ±0.001 (radians), keeping their longitudes at φ = 0 and π . As shown in figure 7, the vortices move far from their equilibrium positions, passing close to the poles. In fact, vortex 1 remains in the northern hemisphere, while vortex 2 remains in the southern hemisphere. Initially, vortex 1 moves south then east, while vortex 2 moves north then west, in such a way that the two vortices remain on a diagonal passing through the equator at φ = π/2 (the trajectory resembles a great circle, though the surface is not a sphere). By t = 400, the vortices have exchanged positions, then they repeat their polar excursions but with a displacement of π in φ. The whole evolution repeats periodically nearly every 800 time units. Notably, two vortices on a sphere are stable [14]. Only a slight oblateness therefore changes the stability properties markedly.
Instability for two vortices also occurs on sufficiently prolate ellipsoids, b > 1.6845 . . . . Here, we illustrate a case with b = 1.69, with the same initial conditions as before. This time the vortices do not move far from their initial positions, so it is more useful to plot φ and η = μ/μ(0) directly as a function of time. This is done in figure 8 for φ 2 (t) (left) and for η 2 (t) (right); vortex 1 exhibits a   symmetric motion, with φ 1 = π − φ 2 and η 1 = −η 2 . The strong, relative growth in η 2 at early times is consistent with linear instability. However, the instability saturates at a small amplitude, then repeats itself periodically.
The periodic behaviour in the two cases illustrated is generic, given the number of degrees of freedom available. For two vortices, there are four coordinates. But for a general surface of revolution, there are two conserved quantities: the Hamiltonian H and the angular impulse M, (4.7). Moreover, H only depends on the longitudinal difference α = φ 2 − φ 1 . This leaves just one degree of freedom, and all that is possible is steady or periodic behaviour.
The motion of three vortices can be decidedly richer. Now there are three degrees of freedom, and therefore the potential for chaotic motion. On a sphere, though, three vortices are stable [14], even to finite-amplitude displacements [17]. Notably, the sphere has two additional conserved quantities (the other two components of the angular impulse), which constrain motions to be steady or periodic. However, a slight departure from a spherical surface leads to instability, both for b < 1 and for b > 1. This is illustrated in figures 9 and 10, respectively, for b = 0.99 and for b = 1.01-both slight deformations of a sphere. The vortex trajectories now exhibit a more irregular motion, potentially chaotic (this has not yet been investigated).

Conclusion
In this paper, we have developed general equations of motion for point vortices, singular concentrations of vorticity, moving on compact surfaces. We have sought to reveal the effects of surface geometry on vortex motion, in particular the self-induced motion associated with variable surface curvature. This self-induced motion stems from the mathematical requirement that the integral of the vorticity must vanish on a (differentiable) compact surface. As a result, a single point vortex must coexist with a background sea of uniform, non-zero vorticity (in fact, it could coexist with any background vorticity which satisfies this mathematical requirement, but a uniform background is the simplest as it reduces the dynamics to the motion of a single vortex). On surfaces of constant curvature like the sphere, the effect of this background sea is not directly apparent, but on surfaces of variable curvature, its effect is to induce vortex motion, at least when the sum of the vortex circulations is non-zero.
As an application, we investigated the linear and nonlinear stability of 'vortex rings', configurations of n equal vortices located in equilibrium on the same co-latitude and equally spaced in longitude, a problem first studied by Thomson [15] in the plane and later generalized to the sphere [14,19]. Here, we presented results for the ellipsoid of revolution, characterized by a vertical : horizontal aspect ratio b. We have found that vortex rings enjoy a wider range of stability on oblate surfaces (b < 1) than on prolate ones (b > 1). As few as two vortices may be unstable on oblate surfaces or sufficiently prolate ones, compared with four vortices for the sphere [14]. Our results correct the analysis in [22], in which the contribution by the Robin function was overlooked.
The formulation of vortex motion presented here appears to readily extend to vortex patches, uniform regions of vorticity bounded by material contours [12,31]. Vortex patch dynamicscontour dynamics [32]-requires only the Green function G, which we have derived for any differentiable compact surface of revolution. We hope to report on this extension in the near future. the hyperboloid). The only other surface with constant curvature is the plane R 2 , and again there is no self-induced motion.
However, compact surfaces are special in that they require the Gauss condition (2.3) for the solution of Poisson's equation (2.2). The simplest way of enforcing this condition is to add a compensating uniform vorticity distribution, as this does not add any degrees of freedom to the dynamics (material conservation of vorticity simply keeps the vorticity uniform). This background vorticity induces a flow, an axisymmetric one on a surface of revolution, which any point vortex must respond to. The only exception is when the sum of the vortex circulations is zero: then there is no background vorticity and no self-induced motion, regardless of the shape of S.
We note in passing that point vortex motion within a domain having exterior and/or interior boundaries exhibits analogous features. The Robin function R, defined as in (2.13), still induces the dynamics of a single vortex through the Hamiltonian −Γ 2 1 R(s 1 )/2. For example, for flow in the upper half-plane, it is simple to show that R(s 1 ) = R(y 1 ) = − 1 2π log(2y 1 ), (B 1) in terms of the usual planar coordinates (x 1 , y 1 ). As a result, the flow is entirely in the x direction, as is well known.