## Abstract

Stationary configurations of identical point vortices on the sphere are investigated using a simple numerical scheme. Configurations in which the vortices are arrayed along curves on the sphere are exhibited, which approximate equilibrium configurations of vortex sheets on the sphere. Other configurations (found after starting from random initial conditions) exhibit net-like distributions of vorticity, dividing the sphere into many cells that contain no vorticity or diffuse vorticity and forming a stationary ‘vortex foam’ on the sphere. They may be viewed as intermediate-energy elements in the set of all identical point vortex equilibria on the sphere. In the continuum limit, these foam states may correspond to stationary states of multiple intersecting vortex sheets. Stationary configurations of point vortices are not found to have this character when vortices of opposite circulations are included.

## 1. Introduction

The dynamical system consisting of point vortices interacting on the sphere is a convenient (although singular) discretization of vorticity that captures many aspects of ideal fluid motion on the sphere in a finite-dimensional setting. Stationary configurations of vortices are interesting as critical points of the interaction energy, and have been the object of much investigation; many references are found in Newton [1] and Meleshko *et al.* [2].

Known stationary configurations of identical vortices fall into two categories, which may be termed as low energy and high energy. ‘Energy’ refers to the interaction energy of the vortices, which is also the Hamiltonian for the system. Low-energy equilibria have point vortices spread out quite uniformly over the sphere [3–5], resulting in low interaction energy. The uniform distribution of the vortices roughly cancels the uniform opposing background vorticity and produces low fluid velocities across the sphere. High-energy solutions have all the vorticity confined to a few curves on the sphere, such as latitudinal circles; symmetry plays a fundamental role in finding these configurations [6]. Physically, the point vortices lined up along curves act as discrete approximations to continuous distributions of vorticity known as vortex sheets [7,8]. Fluid velocity is high over relatively large areas on the sphere, with correspondingly larger interaction energies. These configurations may be compared with relative equilibrium states of vortex sheets on the plane [9,10].

In this article, stationary states of identical point vortices on the sphere are computed using a simple but efficient numerical scheme that converges quickly for small and moderately sized vortex systems. Examples are presented for a system with 100–400 point vortices. Through the use of this scheme, a novel class of stationary configurations is identified—differing from the two types earlier noted—that has most or all of the vorticity located on one-dimensional sets forming a net-like arrangement of vorticity on the sphere. These sets, discrete analogues of continuous vortex sheets, divide the sphere up into numerous cells and could aptly be described as a stationary ‘vortex foam’. Some of these cells can contain diffuse arrays of vortices. The configurations appear to be intermediate to the high- and low-energy equilibria previously known, in two senses. First, the interaction energies lie between the high and low values previously observed. Second, creating a sequence of stationary vortex foams by consolidating the curves and decreasing the number of cells leads eventually to the high-energy states, while increasing the number of cells, and filling some with diffuse vorticity, leads eventually to the low-energy equilibria. Because of the nature of the state space of the vortex system, it will be argued below that these stationary vortex foams may be the ‘most common’ equilibrium configurations.

The paper begins with a transformation of the dynamical equations of the point vortex system into a form more convenient for numerical work. Then, the algorithm for generating the stationary configurations is described, and a number of representative examples are presented. Some evidence is shown that point vortex equilibria, viewed as discrete approximations of continuous vorticity distributions, can indeed converge to vortex sheet equilibria. Others are stationary vortex foams. A brief discussion of the characteristics of these configurations concludes the paper.

## 2. Stationary configurations of point vortices and vortex sheets on the sphere

The equations of motion for a system of point vortices on a sphere have been known for over a century; many references can be found in Newton [1] and Meleshko *et al.* [2]. The derivations of the following paragraphs are presented for the convenience of the reader. Consider the unit sphere centred at the origin in **R**^{3}. The velocity field owing to a point vortex on the sphere can be described as follows: a point vortex of circulation *Γ* at the point (*X*,*Y*,*Z*) induces a velocity at the pole (0,0,1) given by the tangent vector
2.1The velocity field may be expressed in terms of spherical coordinates (*θ*,*ϕ*) or (more conveniently for this work) using complex numbers as in Newton [1]. Stereographic projection of the unit sphere onto the complex plane identified with the equatorial plane {*Z*=0}, taking (0,0,1) to the origin and (0,0,−1) to the point at infinity, projects the point (*X*,*Y*,*Z*) to the complex number *z*=(*X*+i*Y*)/(1+*Z*). Using the complex coordinate *z*, equation (2.1) takes the form
2.2This velocity vector can be identified in the natural way with the complex number . Given a collection of point vortices on the sphere with circulations *Γ*_{1},…,*Γ*_{n} and positions corresponding to the complex numbers *z*_{1},…,*z*_{n} (assuming that all the *z*_{j} are finite), the fluid velocity induced by these vortices at (0,0,1) is zero precisely when the sum is zero.

Now suppose that a configuration of *n* vortices is stationary, i.e. the velocity at each due to the presence of all the others is zero. The sphere, together with its vortex configuration, may be rotated around an axis passing through the equator in a way that takes the *j*th vortex to (0,0,1). This rotation of the sphere corresponds (via stereographic projection) to a linear fractional transformation on the complex plane taking the *k*th vortex from *z*_{k} to . The rotated configuration is still stationary, so that the velocity of the vortex now at (0,0,1) still vanishes. Consequently, we have
2.3Every vortex velocity is zero, so this relation holds for each *j*. If *z*_{j} is infinite (i.e. the corresponding vortex is at the pole *Z*=−1), the limiting form of this relation may be used, namely . Thus, the *n* equations (2.3) form a system of equations in *n* complex variables with solutions that correspond to stationary point vortex configurations on the sphere. This system appears in Newton [1] in a slightly different form.

A suitable criterion for a vortex sheet on the sphere to be stationary can be obtained by adapting the Birkhoff–Rott equation [7,9]. Let *z*(*s*) be a differentiable complex function of the real parameter *s*∈(0,1) that traces out a vortex sheet in the complex plane, with linear vorticity density *Γ*/|*z*′(*s*)| and total circulation *Γ*. The function *z* may be periodic, so that the vortex sheet forms a closed curve. Another possibility is that derivative *z*′ may diverge to infinity as the parameter *s* approaches either end of the unit interval, causing the vorticity density to drop to zero at the ends of the vortex sheet. An example is the uniformly rotating vortex sheet in the plane formed from a uniform elliptical vortex patch in the limit as thickness goes to zero [7].

This sheet is motionless in the *plane* if the following relation holds for all *s*:
(The principal value (p.v.) of the integral is taken.) In the light of equations (2.3), however, the corresponding vortex sheet on the *sphere* will be stationary if, for all *s*,
2.4If several vortex sheets *z*_{j}(*s*) are specified, the condition that all sheets be stationary holds when the following system of equations is satisfied for all *j* and all *s*:
2.5The equations may be modified to accomodate sheets that pass through the point at infinity, but it is no loss of generality to constrain *z*(*s*) to have compact support. Numerically, one can approximate the vortex sheets by arrays of point vortices and simply use equations (2.3).

If none of the vortex sheets intersect, then the integrals in the summation in (2.5) are not singular, and principal values are not needed. If two or more sheets intersect at a point *p*, then we may reparametrize so that they meet at their endpoints with parameter value *s*=0. The requirement of convergence of the integrals at the point of intersection imposes a form of Kirchhoff's point law on the sheets. Suppose that the sheets labelled 1 through to *m* intersect at *p*, and let so that *Γ*_{j}|Δ_{j}| is the (possibly zero) terminal vorticity density of the *j*th sheet. Substitution of a linearization of each sheet in an *ϵ*-neighbourhood of *p*, *z*_{k}(*t*)≈*p*+*t*/Δ_{k}, into equation (2.5) with *s*=0 leads to the relation
Thus, we obtain a necessary condition at intersections
2.6The complex conjugate of equation (2.6) constrains the terminal linear vortex densities *Γ*_{j}|Δ_{j}| and tangent directions arg(*Γ*_{j}Δ_{j}) so that the ‘vorticity currents’ at the intersection point sum to zero. The fact that the vortex sheets are stationary imposes further relations on the intersecting vortex sheets; the planar case is discussed in the study of Bergkvist & Rullgåard [11].

## 3. Numerical method

The solution set to the system of equations (2.3) is known completely in three- and four-vortex systems [12,13], and many symmetrical stationary configurations of identical vortices are also known [6,14]. An algorithm for finding vortex configurations on the sphere that are equilibria with respect to *some* choice of circulations (that is, the circulations are not initially specified) can be found in the study of Newton & Sakajo [15]. Since *relative* equilibria (rotating configurations) are critical points of the energy function, these states may be found numerically by seeking energy-minimizing configurations either through a gradient search (the planar case treated by Campbell & Ziff [16]) or by selective-acceptance random walks [3]. On the sphere, these methods will not distinguish rotating from non-rotating configurations, and convergence is exponentially slow as the minimum is approached. Alternatively, stationary configurations of a large number of identical point vortices can be found by solving system (2.3) numerically through the use of the algorithm discussed below. This algorithm applies more generally to any vortex system with specified circulations, including circulations of differing sign. It gives very rapid (quadratic) convergence at the solutions and allows stationary configurations to be computed to arbitrary precision.

Denote the real and imaginary parts of the expressions on the right-hand sides of equations (2.3) by *f*_{i}, *g*_{i}, respectively, and let **F**=(*f*_{1},*g*_{1},…,*f*_{n},*g*_{n}), **z**=(*x*_{1},*y*_{1},…,*x*_{n},*y*_{n}). The Jacobian matrix of partial derivatives is **F**′. Starting with some initial value **z**=**z**_{0}, define the sequence **z**_{i} as follows:
3.1The positive parameter *λ* controls the rate and region of convergence of the scheme. The value *λ*=1 is Newton's method, while the limit *λ*→0 corresponds to Davidenko's method, i.e. following the gradient flow of −∥**F**∥^{2} to the origin. In place of the inverse matrix, it may be more practical to use a truncated pseudoinverse matrix that is computed by rounding to zero all singular values that lie below some threshold. Under suitable conditions, this sequence converges fairly rapidly.

The only complication encountered when following this scheme is a consequence of the degeneracy of the system (2.3) with respect to rotations on the sphere. Any vortex configuration on the sphere satisfying (2.3) may be rotated so as to put the first vortex at (0,0,1) and the second vortex on a given meridian of the sphere, making *z*_{1} vanish and giving *z*_{2} a prescribed argument. In effect, the system (2.3) of 2*n* real equations has only 2*n*−3 independent real variables. Thus, the matrix **F**′ is singular at every solution to the system (2.3), with nullspace dimension at least 3. So long as only three singular values of **F**′ are small, however, Newton's method (implemented with a truncated pseudoinverse) will compute the solution configuration with quadratic convergence, i.e. the error is squared at each step, and to any desired degree of accuracy.

In the computations below, the sequence (3.1) was computed with values of *λ* between 0.2 and 0.8, switching to *λ*=1 for the last few terms of the sequence to speed convergence. Singular values of **F**′ less than 10^{−5} times the largest singular value were rounded to zero before computing the pseudoinverse matrix. The sequence was deemed to have converged when the residual ∥**F**(**z**_{i})∥ (i.e. the sum of the magnitudes of all the vortex velocities) was reduced to less than 10^{−12} and precisely three of the singular values of **F**′ are less than 10^{−6}. This criterion gives assurance that the computed configuration is indeed stationary, not rotating. As an additional verification of the validity of these solutions to the system (2.3), several of the solutions were recomputed using higher precision arithmetic, forcing the residual below 10^{−50}. For each solution, the Hamiltonian *H* of the system
was computed, which represents the interaction energy of the point vortices on the sphere. For every example below except the one pictured in the final figure in this article, all circulations are 1. The normalized Hamiltonian *H*/*n*^{2} is more useful for comparison of configurations with differing *n*, but even then there will still be a dependence on *n*. For example, the even distribution of point vortices along a great circle, all with circulation 1, is an equilibrium configuration. The values of 8*πH*/*n*^{2} for these configurations of *n*= 100, 200 and 400 unit-circulation vortices are −0.0461, −0.0265 and −0.0150, respectively, tending to zero in the limit . At the other extreme, the energy associated with a uniform distribution of vorticity across the entire sphere is −0.19315 [3], so we could expect that point vortex configurations can take on slightly lower energies.

## 4. Stationary configurations approximating vortex sheets

The first few examples of point vortex equilibria displayed in figures 1–3 are ones that have all the vortices lying along curves, forming a discrete approximation to a vortex sheet equilibrium. The vortex positions are displayed as projected on the complex plane, *z*=(*X*+i*Y*)/(1+*Z*). Each small open circle is centred at a vortex position. The projection distorts inter-vortex distances, as two points on the sphere separated by a small distance *δ* will be separated on the complex plane by distance *δ*(1+*r*^{2})/2, where *r* is the distance to the origin. Another aspect of the projection is that every circular arc on the sphere projects to a circular arc or straight line segment on the plane.

Figure 1 displays a stationary configuration of 240 point vortices on the sphere. Slightly more than half the vortices are on a closed curve that is a slightly distorted circle on the sphere, the remaining 100 vortices are on two vortex sheets in the centre of the figure. These sheets are nearly circular arcs on the sphere, and the vortex density decreases near the ends. It is possible that this configuration is a discrete approximation to a stationary vortex sheet configuration. To test this conjecture, similar point vortex solutions to equations (2.3) were computed with *n*=120, 240 and 360. These configurations are superimposed in figure 2, smaller circles corresponding to larger *n*, showing close alignment that is consistent with the existence of a vortex sheet solution. The normalized energy increases with increasing *n*, as was observed for the great circle configuration earlier.

There are likely a great variety of vortex sheet equilibria, with corresponding point vortex equilibria. Not only are the sheets observed to take on various shapes on the sphere (similar to the situation described in O'Neil [9]), but also the distribution of vorticity between sheets can vary continuously. This is illustrated in figure 3, where four stationary configurations of 200 identical point vortices are displayed, with different distributions of vortices between the various sheets in the solutions. The configuration pictured on the upper left consists of four vortex sheets, each a great semicircle on the sphere, that meet symmetrically at two opposite points. Additional solutions were computed for larger *n*, and a pattern of convergence similar to that in figure 2 was observed.

The next set of solutions to equations (2.3) were computed by starting with a random vortex configuration, and the subsequent stationary configurations are more complex. To give a better view of the vortex distributions on the sphere, these configurations are displayed by projecting the upper and lower hemispheres onto the complex plane using complementary projections; the centres of the left and right parts of each figure correspond to the points (0,0,1) and (0,0,−1), respectively. The unit circle in each projection is the intersection of the plane and the sphere so that vortices on the equator *Z*=0 are displayed twice. Vortices plotted on the outside of the unit circle are located in the other hemisphere, so appear inside the circle in the complementary part of the figure. A vortex sheet that crosses from one hemisphere to the other will be seen to cross the unit circle in the left and right halves of the figure at the same relative position. As before, the stereographic projection results in some radial distortion of the vortex separations: inter-vortex distances are accurately portrayed at the unit circle, but are halved at the origin.

Figures 4–6 show stationary arrangements of identical point vortices that were computed using (3.1) starting from random initial configurations. Specifically, each vortex location was chosen independently using a random variable that is uniformly distributed on the sphere. These random configurations tend to contain both large voids (regions of no vorticity) and vortex clusters in which the vortex separations are very small. The initial configuration was then projected onto the plane and the numerical method applied. In the resulting stationary configuration, there are no vortex clusters, and the distance from each vortex to its nearest neighbour along the vortex sheets seems to be relatively regular. The intersecting vortex sheets form a net-like configuration of cells that are generally empty of vorticity, although some contain one or two vortices, and others appear to be filled with a diffuse distribution of vorticity. In figure 5, most cells are empty or contain only a few vortices. Figure 4 has a similar appearance in the Northern Hemisphere (left circle), but the other hemisphere is filled with vortices that display some linear organization and one significant void, but no clear division into cells. Figure 6 shows a configuration with no clear organization into cells at all in the right circle, and a few cells in the upper hemisphere that contain significant diffuse vorticity. The energy of this configuration is slightly lower than the two previous examples.

## 5. Stationary vortex foam

Each of the stationary configurations in figures 4 and 5 could be termed a ‘vortex foam’ since the vorticity is organized predominately into a net-like arrangement of many connected curves, with most (but not all) of the resulting cells containing no vorticity. If stationary states are given this vague classification, it is necessary to explore their relationship to other stationary states, and to the continuum limit. In this section, it is proposed that this relationship be understood in terms of the interaction energy of the vortices, which corresponds to the kinetic energy of the fluid in the continuum limit. The effect of the organization of the vorticity on the fluid motion is discussed using equation (5.1) below. It was noticed in the course of computation of many equilibrium states that the numerical method frequently converged to a vortex foam when the initial configuration was random; a possible explanation for this is suggested. Finally, vortex systems with mixed circulations are discussed.

How do these stationary vortex foams fit into the vast landscape of all stationary states of the vortex system on the sphere? Interaction energy is one possible way to place them along a ranking or spectrum. States with just a few vortex sheets and higher energy (as in figure 1) could be reached along this spectrum by reducing the number of cells in a vortex foam such as that of figure 5, rearranging the vortex positions as necessary to keep all vortex velocities at zero. On the other hand, by increasing the number of cells and by allowing diffuse areas of vorticity within some cells, one approaches the other end of the spectrum in which the vorticity is spread most uniformly across the sphere, the vortices achieve the maximum possible separation, and interaction is the lowest [3–5]. The configuration in figure 6 could be taken as an example of such a stationary vortex state: it has higher energy than the minimum, but is not organized enough to be termed a vortex foam.

The interaction energy of the point vortex system corresponds to that part of the kinetic energy of the fluid flow that depends on the point vortex positions. High-energy configurations will have higher-velocity fluid flow over a greater portion of the sphere than low-energy configurations. The fluid flow itself is strongly determined by the number and the distribution of the vortex sheets. In a vortex foam, the vortex sheets have the general effect of dividing the sphere into distinct cells. In each cell, the fluid velocity field has a single zero in the middle and increases in magnitude to a maximum at the bounding vortex sheets. The vorticity is positive, so the circulation in each cell is counterclockwise when viewed from outside the sphere. The fluid flow can be computed from (2.1), or the flow at complex coordinate *ζ* may be found directly in terms of the complex positions *z*_{j} and the circulations *Γ*_{j} of the point vortices as follows. Suppose that and let *U*(*ζ*) be given by the formula
It is not difficult to show that the fluid velocity at the point on the sphere corresponding to *ζ* has longitudinal (north–south) component *v*_{s}=*Im*(*U*) and latitudinal component , the positive directions being south and east, respectively. The streamlines of the flow of a stationary configuration, when projected onto the complex plane, coincide with the level curves of the function
5.1A plot of streamlines of a stationary vortex foam with both small and large cells is shown in figure 7. It is evident simply from the number of streamlines that large cells have circumferential flows that are larger both in area and in fluid speed. These features would persist in a continuum limit of stationary vortex sheets (if this limit did in fact exist).

One may ask if the point vortex foam configurations have counterparts in the continuum limit. That is, are there stationary vortex foams consisting of many intersecting vortex sheets that divide the sphere into many cells? A decisive answer is not likely to come from numerically generated point vortex configurations. Nevertheless, one could take a configuration similar to that in figure 5 and interpolate additional vortices into the linear structures in an effort to better resolve a vortex sheet. Figure 8 shows one such configuration containing 446 identical vortices. Two interpolation steps were used, with the numerical method applied after each to bring all vortex velocities back to zero. It is not clear at some points whether vortex sheets intersect, or whether one sheet terminates before meeting the other. It has been shown that stationary point vortex configurations in the plane may have two or more stationary configurations that differ only at the point of intersection of the vortex sheets [17]. The figure does suggest the plausibility of stationary vortex foams made up of both intersecting vortex sheets and terminating vortex sheets (with density tending to zero at the ends).

It is a natural next step to consider the stability and probability of vortex foams. Are the configurations stable under small perturbations? The answer is almost certainly no, given the instability of vortex sheets on the sphere [8]. How common (or uncommon) are these stationary foams? The answer might determine whether they would be observed as quasi-steady states of the dynamical system. In this system, there is no momentum so state space and configuration space are the same. The only energy is interaction energy, so that (for identical vortices) state-space volume in a given energy range decreases at both low *and* high energies. For this reason, the most probable states, in the context of configurations of randomly positioned vortices on the sphere (i.e. the infinite temperature state in the canonical ensemble) are those of intermediate energy. Whether the same is true of *stationary* configurations is not known, but if it were, then we could expect the stationary vortex foams to be in some sense the most probable outcome of the numerical method (3.1). This has been observed (at least qualitatively) in the course of computing several hundred equilibria at moderate values of *n*, and is a subject of continuing investigation.

All the vortex systems so far have contained identical vortices. The final example of a vortex equilibrium configuration on the sphere, shown in figure 9, uses a system with 95 positive vortices and 65 negative ones, circulations chosen so that the total circulation is zero. There is a technical reason not to choose vortices with identical magnitudes: the +1, −1 circulation vortex system contains many stationary states in the plane [14,18], and the numerical method tries to collapse mixed sign vortices together, yielding unstable performance. With the vorticities adjusted to curtail this possibility, the resulting stationary configuration exhibits some one-dimensional features. However, the vortex signs alternate along these curves so that there is no fluid velocity jump across them and they do not have the character of vortex sheets. There are few if any distinct cells, although there are large areas with no vorticity. There is also some pronounced clustering of vortices, despite the fact that the circulations do not allow for any exact equilibrium clusters. (The clustering becomes more extreme as the circulations become closer in absolute value.) This configuration is typical of the results for systems with mixed sign vorticity. No vortex foam states were found for these systems. Moreover, the reasoning in the preceding paragraph does not apply here, since there is no lower bound for energy in a mixed-circulation system.

- Received October 20, 2012.
- Accepted November 19, 2012.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.