## Abstract

We investigate the dynamics of *N* point vortices in the plane, in the limit of large *N*. We consider *relative equilibria*, which are rigidly rotating lattice-like configurations of vortices. These configurations were observed in several recent experiments. We show that these solutions and their stability are fully characterized via a related *aggregation model* which was recently investigated in the context of biological swarms. By using this connection, we give explicit analytical formulae for many of the configurations that have been observed experimentally. These include configurations of vortices of equal strength; the *N*+1 configurations of *N* vortices of equal strength and one vortex of much higher strength; and more generally, *N*+*K* configurations. We also give examples of configurations that have not been studied experimentally, including *N*+2 configurations, where *N* vortices aggregate inside an ellipse. Finally, we introduce an artificial ‘damping’ to the vortex dynamics, in an attempt to explain the phenomenon of crystallization that is often observed in real experiments. The diffusion breaks the conservative structure of vortex dynamics, so that any initial conditions converge to the lattice-like relative equilibrium.

## 1. Introduction

The dynamics of *N* point vortices in a plane is a classical problem that goes back to the works of Helmholtz [1], who first described this model as a fluid analogue to the *N*-body problem in celestial mechanics. Fundamentally, point vortices correspond to singularities in an ideal irrotational flow. These singularities in turn characterize the flow itself. Vortex dynamics have been reproduced in many physical experiments, starting with those of Meyer [2] of floating needle magnets in an applied magnetic field. Yarmchk *et al.* [3] obtained stable vortex lattices in a rotating ^{4}He superfluid. Vortex dynamics have also been observed in Bose–Einstein condensates [4–6] and in electron columns confined in a Malmberg–Penning trap [7], which allow for a high degree of control over initial vortex configurations [8]. Vortex dynamics were also reproduced experimentally using a system of small rotating discs [9–11]. See reference [12] for a recent review of these experiments. Some examples of vortices in nature include geophysical flows [13] such as Jupiter's red spot [14], intensifying hurricanes [15] and plasma flows [16].

As first described by Helmholtz, each vortex generates a rotational velocity field that advects all other vortices. This yields a system of ordinary differential equations (ODEs) for vortex centres *z*_{j} which we will refer to as the *vortex model*,
1.1Here, and 2*πγ*_{k} is the circulation of the *k*th vortex; we will assume throughout the paper that *γ*_{k}>0. There is an extensive literature that describes vortex dynamics of a small number of vortices, or special vortex configurations, such as vortices arranged on a ring (see [17–19] for the overview of the field). The vortex model is a Hamiltonian system of many variables. As such, general initial conditions with four or more vortices typically result in chaotic dynamics. However, there are many special arrangements, which lead to the rigidly rotating configurations of vortices, called *relative equilibria* [12,17,20,21]. These equilibria can be either stable [12,22] or unstable [23]. For example, *N* vortices of equal strength arranged uniformly along a ring is a basic relative equilibrium, which is stable for *N*≤6 and unstable for *N*≥8, with *N*=7 being the threshold case [24–26]. The stable relative equilibria are of particular physical importance, and are often observed in experiments, even when starting with arbitrary initial conditions. In real experiments, there is usually some form of dissipation present, which disturbs the underlying Hamiltonian structure of vortex dynamics and, in practice, causes a decay towards a stable lattice-type relative equilibrium over a long time. For example, this process was observed in experiments with magnetized electron columns [7,8,27] and has been dubbed ‘crystallization’ in [8].

Motivated by rotating ^{4}He superfluid experiments [3], Campbell & Ziff [22] classified many stable configurations for a *small* number of vortices of equal strength. Our goal here is to describe the stable configuration in the continuum limit of a *large* number of vortices . The key observation is that the model (1.1) is intimately connected to the following *aggregation model*,
1.2Here, *ω* is the angular velocity of the relative equilibrium. Model (1.2) was recently investigated in studies [28,29] in connection with swarm formations in biology. There is a one-to-one correspondence between the steady states *x*_{j}(*t*)=*ξ*_{j} of (1.2) and the relative equilibrium *z*_{j}(*t*)=e^{ωit}*ξ*_{j} of (1.1). Moreover, as we show below, this correspondence carries over to stability: the equilibrium *x*_{j}(*t*)=*ξ*_{j} is asymptotically stable for the aggregation model (1.2) if and only if the relative equilibrium *z*_{j}(*t*)=e^{ωit}*ξ*_{j} is stable (neutrally, in the Hamiltonian sense) for the vortex model (1.1). While the steady states and their local stability are the same for the vortex and aggregation model, the latter has much simpler dynamics, and many results can be explicitly obtained in the large-*N* limit as we show below.

Throughout the paper, we make use of the fact that *ω*>0, which follows from our original assumption that *γ*_{j}>0 as we now show. Indeed, we have
1.3where *x*_{j}(*t*)=*ξ*_{j} is the equilibrium state of (1.2) corresponding to the relative equilibrium of (1.1). This can be seen by multiplying the right-hand side of (1.2) by *γ*_{j}*ξ*_{j}, summing over *j*, and setting the sum to zero. The quantity is the total angular vortex momentum, and is the angular impulse.^{1}

## 2. Equilibrium and stability of vortices of equal strength

Consider the special case of *N* vortices of equal strength *γ*_{k}=*γ* so that the aggregation model (1.2) becomes
2.1where *a*=*Nγ*_{k}/*ω* and where we rescaled the time The distinguished 1/*N* scaling makes it possible to take the mean-field limit which yields a non-local partial differential equation (PDE) [28–30]
2.2Here, *ρ*(*x*,*t*) approximates the particle density normalized so that represents the fraction of particles inside a domain *D* with .

The system (2.2) was analysed in detail in recent studies [28,29]. It was shown that in the limit the density *ρ*(*x*,*t*) approaches a steady state which is constant inside a disc of radius and is zero outside such a disc: if and as otherwise. For the vortex model, this result implies that for large *N*, the stable relative equilibrium for (1.1) with *γ*_{k}=*γ* consists of particles uniformly distributed inside a disc of radius
2.3The radius *R*_{0} was previously derived in the physics literature, see for example [31] and in [32, eqn. 25]. The stability can also be ascertained: as was shown in [29], the uniform disc is the global attractor for the continuum limit (2.2) of the aggregation model. An immediate consequence is that the uniform disc is the only stable relative equilibrium of the vortex model (1.1) with *γ*_{k}=*γ* in the large *N* limit. While in the derivation of (2.3), we assumed that *R*_{0} is *O*(1), the formula actually holds even when *γ*/*ω*=*O*(1). This is because one can always rescale the space to make *R*_{0} of *O*(1). Indeed, the rescaling , eliminates *a*/*N* from equation (2.3); hence, these results are independent of scaling.

Figure 1 shows that the continuum limit radius (2.3) provides a very good estimate of the equilibrium radius even for relatively small *N*=25; as expected, larger *N* results in an even better agreement. Below, we will re-derive this result as a special case of the *N*+1 configuration.

## 3. *N*+1 problem

A well-studied case in the literature going back to Havelock [25] consists of a ring of *N* vortices of equal strength surrounding a single vortex of a much higher strength. This is the fluids analogue of the ‘rings of saturn’ problem in celestial mechanics that was first studied by Maxwell [33]. Cabral & Schmidt [26] established the stability of an *N*+1 uniform ring, provided that the strength of the central vortex is between a certain lower and an upper bound. Below the lower bound, the ring has high-mode instabilities that cause the ring to deform into an annulus. This is illustrated in the left-top panel of figure 2. Above the upper bound, there is a low-mode instability that causes the ring to deform into a stable non-uniform ring-like state [20], as illustrated in the second row of figure 2. The non-uniform state of this type was also observed experimentally in [9] and in [8] (figure 2). For the experiments with Bose–Einstein condensates, it has been shown that higher charge vortices are both experimentally and theoretically potentially unstable, owing to break-up to lower topological charges [34–36], which makes the *N*+1 configurations more difficult to achieve.

In this section, we will construct the *N*+1 states in the large *N* limit. Label *η*=*x*_{N+1} and rewrite (1.2) as
3.1and
3.2where *a*=*Nγ*_{k}/*ω*, *k*=1…*N*; *b*=*γ*_{N+1}/*ω* and where we rescaled the time . Taking the mean-field limit as before yields a PDE for the density distribution of small vortices *ρ*(*y*) which satisfies
3.3We make the following ansatz for the steady state:
where *B*_{R}(*z*) denotes a disc of radius *R* centred at *z*; and *η*∈*B*_{R0}(0) are chosen such that *B*_{R1}(*η*)⊂*B*_{R0}(0). The normalization yields . Using the identity
we find that for *x*∈*B*_{R0}(0)\*B*_{R1}(*η*),
At the steady state of (3.3), *v*(*x*)=0=*η*_{t} so that *aρ*_{0}*π*=1; . Solving for *R*_{1},*R*_{0}, we obtain the following result.

### Result 3.1

Define , and suppose that *η* is any point such that *B*_{R1}(*η*)⊂*B*_{R0}(0). Then, the equilibrium solution for (3.3) is constant inside *B*_{R0}(0)\*B*_{R1}(*η*) and is zero outside.

This result is illustrated in figure 2, first row, where the boundaries of discs *B*_{R1}(*η*),*B*_{R0}(0) are shown by dashed lines. Excellent agreement between the continuum limit result is observed. Unlike the *N*+0 problem, the relative equilibrium for the *N*+1 problem is non-unique: any choice of *η* yields a steady state as long as |*η*|<*R*_{0}−*R*_{1}.

In the limit , note that so that the steady-state approaches a (possibly non-uniform) ring solution. Assuming that *η*∼0 is at its equilibrium, it was shown in [20] that the evolution of the small vortices is given by *x*_{j}(*t*)∼*R*_{0}e^{iθj(t)} where *θ*_{j} satisfies, after a rescaling ,
3.4In the mean-field limit , the density distribution *ρ*(*θ*) for the angles *θ*_{j} satisfies
3.5where PV denotes the principal value integral, and . Using a formal expansion
3.6note that up to rotations, the steady-state density *ρ*(*θ*) for which *v*=0 must be of the form
3.7This formula was first derived in [37] using a similar technique. The free constant *α* corresponds to the non-uniqueness of the steady state of the continuum limit (3.5) of the *N*+1 problem. This steady state is in fact stable (see appendix A), and the density is strictly positive whenever |*α*|<1. However, for any *finite* *N*, the density *ρ* is quantized; this has the effect of choosing a specific constant *α*. The discrete locations are well approximated by
3.8(figure 2, row 4). The relative equilibrium exhibits a gap along the unit circle. In appendix A, we show the following result.

### Result 3.2

In the limit *a*≪1 and *N*≫1, the equilibrium *θ*_{j} for the problem (3.4) is approximated by solving (3.8), where *α*∼1+*A* *N*^{−1/3} with *θ*_{1}∼−*π*+*B* *N*^{−1/3}, where the constants *A*≈2.0699802 and *B*≈4.122044 satisfy equations (A4). The gap in the steady state, given by *Θ*_{gap}=2(*θ*_{1}+*π*), is of the size
3.9This result is illustrated in figure 2, second row. Despite the low-power scaling *O*(*N*^{−1/3}) for the gap, formula (3.9) is very effective even for relatively few vortices *N*=5.

## 4. *N*+*K* problem

It is straightforward to generalize result 1 to the situation where there are *K* large vortices and *N* small ones. In analogy to the *N*+1 case, we let *a*=*Nγ*_{k}/*ω*, *k*=1…*N*; *b*_{k}=*γ*_{N+k}/*ω*, *k*=1…*K* and rescale the time . The mean-field limit for the velocity *v* of the density distribution *ρ* of vortices then becomes a coupled system
4.1The generalization of result 1 then yields

### Result 4.1

Let , *k*=1…*K* and . Suppose *η*_{1}…*η*_{K} are such *B*_{R1}(*η*_{1})…*B*_{RK}(*η*_{K}) are all disjoint and are contained inside *B*_{R0}(0). The equilibrium density for (4.1) is constant inside and is zero outside.

In analogy to the *N*+1 problem, such an equilibrium has *K* free parameters *η*_{1}…*η*_{K}. Some examples are given in figure 3. Unlike the *N*+1 problem, not all values of parameters are allowed in result 3. For example, for the *N*+2 problem, such solutions exist if and only if *R*_{1}+*R*_{2}<*R*_{0} or . In the opposite case, the two smaller discs no longer fit inside the big disc; this is illustrated in figure 3, row 2.

When *K*>2 and in the limit , the integral terms in (4.1) disappear at the leading order, so that *K* large vortices ‘decouple’ from the *N* small vortices, and their behaviour is simply given by a lower-dimensional *K*-vortex problem obtained by setting *a*=0. The *N* small equilibria aggregate at several points as illustrated in figure 3 (middle left and bottom). These special points correspond to the stagnation points of the velocity field (4.1) with *a*=0. As we now show, the shape of the swarm around these points can be computed explicitly. For simplicity, we concentrate on a special case of two big vortices of equal strength *b*_{1}=*b*_{2}=*b* and one aggregation of small vortices. This configuration is shown in figure 3, middle right panel. The relative equilibrium position of the two big vortices *η*_{1},*η*_{2} satisfies
By rotation, assume *η*_{1},*η*_{2} are on the *x*-axis so that . For the small vortices and in the limit , we set *x*_{j}=*X*+*a*^{1/2}*ξ*_{j}+*O*(*a*). Collecting the *O*(1) terms, we then obtain
This yields . At the next order, collecting the *O*(*a*^{1/2}) terms, we then obtain system
4.2where the bar denotes the complex conjugate. The system (4.2) was recently analysed in [38] in the context of singularly perturbed aggregation kernels. Using a complex variables method, it was shown that *ξ*_{j} concentrate uniformly inside an *ellipse* whose radii are and and whose major axis is parallel to the *x*-axis. This ellipse, along with the full solution, is plotted in figure 3, right-middle panel; good agreement is observed.

## 5. Crystallization

Owing to the conservative nature of the idealized vortex dynamics (1.1), random initial conditions typically result in chaotic motion, and stable relative equilibria are only neutrally stable: small perturbations neither grow nor decay. By contrast, in experiments [3,7–9], arbitrary initial conditions often converge to an asymptotically stable lattice-like state. Similar stable lattices have also been observed in Bose–Einstein condensates [4]. The emergence of such lattices can be explained by adding a small ‘damping’ term that destroys the Hamiltonian structure. Here, we propose the following phenomenological model that incorporates this friction:
5.1The constant *μ*≥0 models damping effects. The term *γ*_{k}(*i*+*μ*)(*z*−*z*_{k}/|*z*−*z*_{k}|^{2}) corresponds to a velocity field generated by an outwards spiralling source that adds local repulsion of strength *μ*. The term −*μωz*_{j} keeps the vortices confined near the origin; it arises naturally for vortices confined to a circular domain via boundary effects. The original vortex model (1.1) is a special case of (5.1) obtained by setting *μ*=0. Moreover, the aggregation model (1.2) is also a special case of (5.1) by taking after rescaling the time .

The specific form of (5.1) is motivated by the fact that any relative equilibrium of (5.1) is also a relative equilibrium of (1.1) for any *μ*; and vice versa. Furthermore, we now show that the relative equilibrium of the vortex model (1.1) is stable if and only if it is a stable for the damped system (5.1), and if and only if the corresponding equilibrium *x*_{j}(*t*)=*ξ*_{j} of the aggregation model (1.2) is stable. Note that systems (1.1), (1.2) and (5.1) are all invariant under rotations; hence, there is always a zero eigenvalue of the relative equilibrium that corresponds to the rotation invariance. So, we define stability to mean that all *other* eigenvalues are non-positive. Because the vortex model (1.1) is Hamiltonian, its eigenvalues come in pairs ±*λ*. The (neutral) stability in this case means that all eigenvalues (except for the zero rotational mode) are purely imaginary.

We linearize around the relative equilibrium by setting where *η*≪1. We then obtain the 2*N*×2*N* linear system of ODEs
5.2where the overbar denotes complex conjugation; *η*=(*η*_{1},…*η*_{N}) and *L* is an *N*×*N* matrix with for *j*≠*k* and with .

Eliminating from (5.2), we obtain
5.3Setting in (5.3), we find that *λ* satisfies
5.4where *σ* is an eigenvalue of with corresponding eigenvector *v*. Note that where is a positive diagonal matrix, and is a symmetric matrix. It follows that is similar to the hermitian matrix , so that . Hence, recalling the assumption *μ*,*ω*>0, we obtain that *Re*(*λ*)<0 if and only if *σ*<*ω*^{2}.

The original vortex model (1.1) corresponds to taking *μ*=0 so that . If *σ*<*ω*^{2}, then *λ* is purely imaginary which corresponds to the neutral stability. Otherwise, is a positive eigenvalue, and the relative equilibrium is unstable.

An analogous computation for the aggregation model (1.2) shows that the eigenvalues of the steady state *x*_{j}(*t*)=*ξ*_{j} are similarly given by
5.5so that with both *λ*<0 if and only if *σ*<*ω*^{2} (note that (5.5) can also obtained from (5.4) by taking the limit after rescaling *λ*→*μλ*). This establishes the equivalence of linear stability of the steady state *x*_{j}(*t*)=*ξ*_{j} of the aggregation model (1.2), the neutral stability of the corresponding relative equilibrium *z*_{j}(*t*)=e^{iωt}*ξ*_{j} of the vortex model (1.1), and the linear stability of the model with damping (5.1). We remark that this method is similar to the argument that was used in [39].

## 6. Discussion

Vortex dynamics is a very old subject, dating back to the 1850s. In this paper, we explored the connection between vortex dynamics and a model of biological swarming. This connection allows us to obtain many new as well as existing results in the large *N* limit. For example, previous works on the *N*+1 configurations [20,26] consider only the case where the *N* vortices form a ring around the big vortex. Here, we treat the more general case where the equilibrium is bounded between two circles. The previously known cases can be recovered as limiting case when the radius of the inner circle approaches the radius of the outer circle. Moreover, we also computed asymptotically the gap in the *N*+1 non-uniform ring using a similar approach. This result is also new.

Many open questions remain. The stability of *N*+*K* configurations has not been studied, except when *K*=0 or for a very restricted case of the *N*+1 problem (see appendix A). For the *N*+2 problem, the steady state as described in result 2 requires that the two discs inside the swarm are disjoint. But there are other configurations, such as those illustrated in figure 3, top-middle panel, where the two discs overlap by an *O*(1) amount. Numerical simulations indicate that the amount of the overlap is independent of initial conditions, provided that the initial conditions are such that the small vortices are all on one side of the two big vortices.

Numerical simulations of the crystallization model (5.1) indicate that the steady state attained in figure 4 is a global attractor. This remains to be shown. Model (5.1) is the simplest possible model that incorporates damping while still preserving the relative equilibrium of the original vortex model (1.1). It would be interesting to incorporate the damping effects in a more systematic way starting from first principles. See [40] for work in this direction.

## Acknowledgements

We thank the anonymous referees for a careful reading, very useful suggestions including formula (1.3), and for finding several mistakes in the first draft of the paper. Their comments helped to improve the paper significantly. We thank James von Brecht for fruitful discussions. T.K. was supported by a grant from AARMS CRG in Dynamical Systems and NSERC grant no. 47050. Some of the research for this project was carried out while Y.C. and T.K. were supported by the California Research Training Programme in Computational and Applied Mathematics (NSF grant no. DMS-1045536).

## Appendix A

**(a) Stability of the N+1 solution (3.7)**

We now show that (3.7) is a stable steady state of the continuum equation (3.5). For simplicity, we restrict the discussion to symmetric solutions that have the Fourier expansion
Using (3.6), the velocity is then given by
Substituting into *ρ*_{t}+(*ρv*_{θ})_{θ}=0, we obtain the following infinite system of ODEs for the Fourier coefficients
A1This system admits a steady state *C*_{m}=0, *m*≥2 with *C*_{1} arbitrary. Linearizing about this state yields the infinite Jacobian matrix

Applying the Gershgorin Theorem to the columns of *J*, its eigenvalues are non-positive provided that |*C*_{1}|<*C*_{0}/2, which is precisely the condition that steady-state density is positive.

**(b) Derivation of result 2**

The continuum density is approximated by a discrete density , where *θ*_{k} are given by (3.8). The idea is to choose *α* such that the velocity at the boundary is zero: d*θ*_{1}/d*t*=0, where d*θ*_{1}/d*t* is given by (3.4). We then estimate the discrete sum in (3.4) by (d*θ*_{1}/d*t*)∼
where *θ*_{1+1/2} is given by (3.8) with *j*=1+1/2; refer to the bottom of figure 2. This integral can be evaluated exactly using the identities
We then expand using the following scaling:
where *A*,*B*,*C* are *O*(1). Setting d*θ*_{1}/d*t*=0 and expanding in Taylor series yields, at the leading order,
A2The expressions (3.8) for *x*_{1} and *x*_{1+1/2} yield to leading order,
A3Together, (A2) and (A3) give a system of three algebraic equations for *A*,*B*,*C* independent of *N* to leading order. Two of the variables can be eliminated by defining *u*=*C*/*B*; *v*=*A*/*B*^{2}. Eliminating *v*, we then obtain
A4

- Received February 8, 2013.
- Accepted May 22, 2013.

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