## Abstract

Motivated by the recent successes of particle models in capturing the precession and interactions of vortex structures in quasi-two-dimensional Bose–Einstein condensates, we revisit the relevant systems of ordinary differential equations. We consider the number of vortices *N* as a parameter and explore the prototypical configurations (‘ground states’) that arise in the case of few or many vortices. In the case of few vortices, we modify the classical result illustrating that vortex polygons in the form of a ring are unstable for *N*≥7. Additionally, we reconcile this modification with the recent identification of symmetry-breaking bifurcations for the cases of *N*=2,…,5. We also briefly discuss the case of a ring of vortices surrounding a central vortex (so-called *N*+1 configuration). We finally examine the opposite limit of large *N* and illustrate how a coarse-graining, continuum approach enables the accurate identification of the radial distribution of vortices in that limit.

## 1. Introduction

The advent of Bose–Einstein condensates (BECs) has offered intriguing twists to a number of explorations regarding nonlinear waves, their dynamics and their mutual interactions [1–3]. The realm of point vortices constitutes a canonical example of this type. Their exploration has been a fascinating topic, garnering considerable attention starting from the fundamental contribution of Lord Kelvin [4], extending to their critical role in turbulent dynamics proposed by Onsager [5] and reaching up to more recent explorations in a diverse range of fields. The latter include (but are not limited to) patterns forming in rotating superfluid ^{4}He [6], electron columns confined in Malmberg-Penning traps [7] and even magnetized, millimetre-sized discs rotating at a liquid–air interface [8,9]. Numerous theoretical advances have also been made by considering the ordinary differential equations (ODEs) describing the vortex particles. Among them, we briefly note the classical examination of the stability of vortices in ring formation [10], the consideration of higher numbers of vortices in [11,12] and even of asymmetric equilibria thereof in [13]. Much of this activity has been summarized, for example, in the review [14] and in the book [15]. A more recent exposition tying the classical fluid point-vortex problem of Aref *et al*. [14] and Newton [15] with the BEC-oriented motivation of this work is given in [16].

In the context of BECs, vortices have been one of the focal themes of numerous investigations, which have now been summarized in many reviews [17–21]. A considerable volume of associated experimental efforts was concerned with the creation schemes of such structures [22–26], especially in the form of unit charge vortices. However, there were also efforts to produce the (dynamically unstable) vortices of higher topological charge [27,28] and observe their decay. Furthermore, numerous studies focused on vortex lattices with a large number of vortices [29–31]. Recently, the significant advancements of experimental vortex creation and visualization techniques have led to a new thrust towards the study of small clusters of vortices. While this direction was initiated early on with the creation of few same-circulation vortices in [32] (and their theoretical examination in [33]), it has recently ignited considerable interest due to the systematic formation and observation of counter-circulating vortex dipoles [25,34–36], tripoles [37] and also sets of 2-, 3-, 4- (and more generally controllably few) vortices [38].

One of the features that are especially interesting in connection to this BEC context concerns the relevance and usefulness of the modelling of the vortices as point particles whose positions are described by ODEs. Such an approach has been used in order to not only offer a detailed quantitative description of the vortex dynamics (in comparison to the experiment) as in the case of the vortex dipole [35,36], but also as a tool to unveil subtle bifurcation and symmetry-breaking phenomena as in the case of few co-rotating vortices [38] (again corroborating experimental observations). It is in that light that we consider a deeper understanding of the features of such ODEs of particular relevance and importance within this system. On the other hand, in connection to the classical and intensely studied point-vortex problem, the BEC setting offers an intriguing twist. Namely, not only should one consider the pairwise interaction between the vortices, but the phenomenology is significantly affected by the precessional motion of *each* vortex within the parabolic external trap confining the BEC atoms. It is the combination of these two key features that gives rise to numerous unprecedented phenomena, such as the existence of an equilibrium vortex dipole (with the vortices located at a suitable distance from the trap centre), or the destabilization of vortex polygons for small vortex number *N*.

A natural question is the one concerning the applicability of this point-vortex method to BECs and the potential advantages that this approach may hold in comparison to the well-established approach of the mean-field so-called Gross–Pitaevskii partial differential equation (PDE) model of this setting. While works combining theoretical observations based on this method and actual laboratory experiments and comparisons between the two [25,35,36,38] strongly suggest the usefulness of the method, let us add a few more comments in that regard. In particular, comparisons of the vortex-particle ODEs with the Gross–Pitaevskii PDE have been given not only for simpler, oscillatory dynamics, but also for more complex chaotic dynamics (e.g. [39] as a recent example). Finally, not only cases where only condensate atoms are present have been considered, but more recently cases with thermal atoms have also been studied [40]. From these works, there is an emerging understanding of the settings where this point-vortex approximation may be most suitable. In particular, large chemical potentials make the vortex increasingly more localized and hence its internal structure progressively less relevant. Furthermore, the quantum pressure term is effectively accounted for in the form of these models and should not *a priori* pose a problem. Perhaps the most intricate and less controllable aspect is that of the vortex–sound interactions, which partially depends on the precise form and the symmetry of the initial conditions; for a relevant, very recent discussion, see Yan *et al*. [40] and references therein.

This work aims to explore some of these features in the case of co-rotating vortices (i.e. vortices of only one charge sign). This is the most typical experimental situation, given that most set-ups create the vortices through the imparting of angular momentum to the system [19]. More specifically, we intend to examine the two opposite limits of experimental tractability:

— On the one hand, we consider the case of small vortex numbers

*N*. In this case, the canonical expectation is that the vortices will lead to the formation of a polygonal ring, given the dynamical stability of such a ring. Here, we will give a systematic stability analysis of the ring formation generalizing the classical work of Havelock [10] (see also [41]). However, this will have two important side conclusions. On the one hand, it will be shown that while the classical result is that vortices become unstable for*N*>7 in the ring formation, here due to the precessional term, even the*N*=7 case is*always*unstable, non-trivially modifying the classical result. However, there are more surprises; we will see that the eigenvalues of the ring formation critically depend on the ring radius and may even lead to instability for the cases*N*=2–6, whereas the work of Havelock [10] predicts generic stability for the classical point-vortex problem. In fact, it will be argued that these destabilization events are exactly the ones recently identified by varying the angular momentum of the vortex system in [38].— On the other hand, the opposite limit of large vortex number is equally interesting (and experimentally accessible, per the vortex lattice experiments discussed above). In that case, we present a coarse-graining description developing a continuum model for the vortex distribution and its stationary form. We identify the radial form of this distribution. By finding the stable equilibrium of the ODEs for the case of large

*N*and developing a vortex counting algorithm that enables the identification (from the particle results) of this distribution, we obtain excellent agreement with the prediction of the coarse-grained model.

We believe that these findings will shed light on the theoretical analysis of vortices in quasi-two-dimensional BECs, identifying some of the complications and subtleties arising due to the presence of the (critical for the present dynamics) precessional term. As an aside, we also show how some of the relevant techniques can be generalized in other settings, e.g. by computing the stability of the so-called *N*+1 configuration, where *N* vortices form a ring, while 1 is located at the trap centre. It should be noted here that a principal motivation of this work concerns settings other than the ones where the rotation of the entire condensate has rendered a multi-vortex state the ground state of the system (minimizing its free energy due to the presence of the vortices). More specifically, instead, we envision a situation such as that of [25,34–36] or [38] where multiple vortices have been created through a suitable external driving or quenching of the condensate, yet they represent an excited, potentially dynamically stable state of the system.

Our presentation is structured as follows. In §2, we offer the general theoretical set-up in line with the earlier works [35,38] that corroborated its correspondence with experimental results. In §3, we present the analysis of the polygonal vortex ring and how it connects to the stability conclusions of recent few-vortex studies [38,42]. In §4, we explore the opposite limit of large *N* and the corresponding coarse-grained, continuum model. Again a comparison is offered, this time with numerical computations within that limit. Finally, in §5, we summarize our findings and present a number of conclusions and possible future directions. Appendix A presents a number of technical details associated, for example, with the stability of the *N*+1 vortex configuration.

## 2. Theoretical set-up

The starting point for our considerations will consist of the vortex equation of motion which can be written in the form of a single complex-valued ODE for *z*_{j}=*X*_{j}+*iY* _{j}, where (*X*_{j},*Y* _{j}) denotes the planar position of the *j*th vortex. This equation reads
*S*=+1), as in the context of the recent experiments of Navarro *et al.* [38]. While for many of our considerations, we will keep this precession term as general as possible, for a number of concrete calculations we will assume the same form as used in the recent experimental considerations of Middelkamp *et al.* [35] and Navarro *et al.* [38]; see also the detailed analysis/comparison with single vortex precession experiments in [25]. In particular, in line with these works, we will choose
*a* representing the precession frequency of the vortex very near the centre of the trap can be absorbed into a rescaling of time (rendering time dimensionless).

On the other hand, the second term in equation (2.1) represents the vortex–vortex interaction, i.e. the velocity field at the location of a vortex due to the presence of other surrounding vortices. Our adimensionalization of the model follows that of Navarro *et al.* [38], where it was explained that a ‘typical’ experimentally relevant value for *c* is 0.1. It should be noted here that this constant is effectively taking into account the non-uniformity of the background density (due to the presence of the trap) through which the vortices are interacting. However, one can consider more elaborate (integral) functional forms, explicitly taking into account the density modulation, as described in [43]. Given the successes of the simpler set-up in capturing the recent experimental observations and the amenability of the corresponding functional forms to our analytical considerations, we will indeed proceed to consider the precession and interaction terms as presented in equation (2.1).

It is natural to expect in the considered case of co-rotating vortices that both the precession and interaction effects lead to rotation of the vortices in the same direction. In that light, no genuine equilibria may exist but instead we may seek relative (in a rotating frame) equilibria of the form *z*_{j}(*t*)=*e*^{iωt}*x*_{j}. These satisfy the equation

As a (partially numerical) aside, we will also consider the following ‘aggregation’ equation:
*z*_{j}(*t*)=*e*^{iωt}*x*_{j} of equation (2.1) corresponds to the relaxational equilibrium *x*_{j}(*t*)=*x*_{j} of equation (2.4) and vice versa. Moreover, there is an intimate connection between the stability of the two models. In appendix A(a), we prove the following result.

### Theorem 2.1

*Suppose that an equilibrium x*_{j}*(t)=ξ*_{j} *of equation* (*2.4*) *is stable. Then the relative equilibrium z*_{j}(*t*)=*e*^{iωt}*ξ*_{j} *of equation* (*2.1*) *is* (*neutrally*) *stable. The converse is also true in the following two cases*: (*i*) *f*(*r*)=*const. or* (*ii*) |*ξ*_{j}|=*const. for all j and f*(*R*) *is decreasing.*

Case (i) of this theorem was shown in [44]; this is reproduced in appendix A(a). However, the proof in Ref. [44] does not work for a general case of non-homogeneous *f*(*r*), and a more general approach is taken here. On the flip side, we can only show the stability of equation (2.4) implies stability of equation (2.1); we do not know how to prove the converse (nor do we have counter-examples).

## 3. Ring solutions

The prototypical configuration that one expects to identify as a stable equilibrium in the case of small vortex number *N* (motivated by the corresponding result in the absence of precession [10–12]) is the ‘ring’ configuration with the vortices sitting at the vertices of a canonical polygon. Assuming such a relative equilibrium to equation (2.1) with radius *R*, we can write it in our complex notation as
*R* satisfies

As case (ii) of theorem 2.1 (of appendix A(a)) shows that the ring for the vortex model equation (2.1) is stable if and only if it is stable for the aggregation model equation (2.4). The full characterization of linear stability is given by the following theorem (following the procedure used in [45,46]).

### Theorem 3.1

*Consider the ring solution for equation* (*2.1*), *of radius R as given by equation* (*3.1*), *where the frequency ω is given by equation* (*3.3*). *Suppose that N is odd. Then the ring is stable provided that*
*and is unstable if the inequality is reversed. Suppose that N is even. Then the ring is stable provided that*
*and is unstable if the inequality is reversed.*

*The ring is generically unstable if N≥7 and f*(*r*) *is an increasing function.*

### Proof

Because of theorem 2.1 case (ii), it is sufficient to consider the stability of the steady state *h*_{k} is a small, complex-valued, perturbation. After some algebra, we obtain for the evolution of the small perturbations
*e*^{im2πj/N} and *e*^{−im2πj/N}, system (3.4) decouples into a sequence of 2×2 subproblems
_{−}(0)=0; this mode corresponds to rotation invariance. Moreover, (*m*−1)(*N*−*m*−1) is positive for all integers 1≤*m*≤*N*−1 and attains its max at *m*=*N*/2. Therefore, λ_{−}(*m*) correspond to stable eigendirections and the instability threshold is obtained by setting λ_{+}(⌊*N*/2⌋)=0, where ⌊⋅⌋ denotes the floor function.

Suppose that *N* is odd. Then we substitute *m*=(*N*−1)/2 into equation (3.11) to obtain
*N* is even, we substitute *m*=*N*/2 which yields
*N* odd and even) is negative and unstable otherwise. If *f*^{′}=0, this recovers the *N*=7 threshold. ▪

Some observations are important to make here. In particular, we wish to examine the special cases of *N*=2,…,7 which are well known to be canonical examples where the polygonal ring is of interest as a configuration even in the absence of precession.

We start with the *N*=7 case. We can see that contrary to the case where the precession is absent (wherein this case is *critical* with the corresponding eigenvalue being neutral), here the presence of precession has a crucial impact. In particular, it adds a positive part (for *f*(*R*) increasing, as is the case in our BECs) to the pertinent eigenvalue rendering the corresponding configuration *generically* unstable. This is a particular trait of our BEC vortex ring configuration.

For the cases with *N*<7, we can express the corresponding eigenvalue in the general form λ=*f*^{′}(*R*)*R*+(*c*/8*R*^{2})(*N*^{2}−8*N*+*s*), where *s*=7 for *N* odd and *s*=8 for *N* even. A key observation now (which directly connects this work with the experimental observations of Navarro *et al*. [38] and the computational/theoretical analysis of Zampetaki *et al*. [42]) is that although for these *N*<7 cases, the configuration is *not* generically unstable, nevertheless it *may* be unstable for sufficiently large *R*. To phrase this differently, we can parametrize the vortex system by the angular momentum of the vortex particles, which is a conserved quantity of the dynamics. The angular momentum is defined as *L*=*NR*^{2}.^{1} In that light and taking *f*(*r*) as in equation (2.2) with *a*=1, the eigenvalue whose zero crossing will determine the potential instability of a configuration with *N*=2,…,6 reads
*L* of the form

Remarkably, this expression yields the relevant critical angular momenta for *all* cases of *N*=2–5 in direct agreement with the expressions given in [38,42] for the symmetry-breaking bifurcations due to the destabilization of the ring configuration. Using the notation
*N*=7, while the dominant eigenvalue is always positive (as indicated above), even the next one λ_{(N−3)/2} can be seen to cross 0 at *N*=6, the relevant critical point had not been previously identified in [38] or [42].

However, a final observation is also in order. As discussed in [42], the interval (of *L*) of dynamical stability of the small *N* configurations does *not* coincide with the interval of *L* for which these constitute the ground state of the system. In the case of *N*=2, the stability threshold and ground-state asymmetry threshold coincide (this is a supercritical pitchfork point), but in other cases such as *N*=3 and *N*=4, asymmetric configurations (such as isosceles triangles for *N*=3 or rhombic configurations for *N*=4) acquire lower energy than the polygonal ring distinctly before its loss of stability threshold [42]. Namely, the ring configuration becomes a local (rather than global) minimum of the energy clearly before its destabilization. Unfortunately, these asymmetric configurations which are stabilized by the presence of our precession term (in its absence such asymmetric configurations are unstable, as discussed, e.g. in [13]), do not have a general closed form that would enable their stability analysis. Nevertheless, another symmetric configuration that emerges as relevant for the ground state of the system, at least in the case of *N*=5 examined in [42] (and obviously also for larger *N*) is the so-called *N*+1 vortex configuration, consisting of *N* vortices on the polygonal ring and one more at the centre. For the classical vortex problem ( *f*=0 in equation (2.1)), the stability of this configuration was analysed in [47]; see also [48]. In appendix A(b), we show the following generalization to the full BEC problem.

### Theorem 3.2

*Consider the N+1 configuration of equation* (*2.1*) *consisting of N vortices uniformly distributed on a ring of radius R and angular velocity ω given by equation* (*3.1*) *plus a vortex at the origin. Then R satisfies*
*Let
*
*where* λ_{+}*(m) as given by equation* (3.11), *and let*
*The N+1 configuration is stable if and only if* *and all eigenvalues of M*_{0} *are negative.*

In particular, this theorem shows that the *N*+1 configuration is stable if the *N*-ring is stable and if in addition the matrix *M*_{0} is negative definite.

When *f*=0, both 9+1 and 3+1 configurations are marginally stable (the former due to *N*=9, the latter due to the eigenvalue of *M*_{0} crossing zero when *N*=3); the *N*+1 configuration is stable for 3<*N*<7 and is unstable for *N*>9 or *N*<3. This is in agreement with the results in [47]. When *f* is increasing as in equation (2.2), both 9+1 and 3+1 configurations lose their marginal stability and become unstable, so that *N*+1 configuration becomes unstable for any *N*≥9 or *N*≤3.

## 4. Continuum limit

Having considered the case of small *N*, we now explore the opposite limit of large *N*. Note, however, that our results for the ring and the *N*+1 vortex configuration are entirely general and the dynamical stability thereof is obtained for any *N*. Yet, these configurations can only be stable when *N* is sufficiently small, as discussed above, e.g. *N*<7 for the ring state (and even then for sufficiently small *L*). Hence, in the opposite limit of large *N*, we expect that a substantially different vortex distribution will arise. This expectation is confirmed not only by the well-known vortex lattice observations of, for example, [29–31], but also by the direct numerical evolution of the aggregation equation (2.4) (figure 1*a*,*b*). Recall that the aggregation equation has the benefit of relaxing to equilibrium attractors corresponding to the marginally stable equilibria of our original equation (2.1). It is then particularly relevant to attempt to identify the distribution of such a large number of vortices *N*. Here, we will use techniques similar to Ref. [49] to derive the limiting density profile.

Following, the discussion of Bernoff & Topaz [50], we coarse-grain by defining the particle density to be
*δ*(*x*) is the Dirac-delta function. It is then straightforward to rewrite the aggregation equation (2.4) as *f*(*r*)=*a*/(1−*r*^{2}) is relevant, in particular, to the precessional dynamics of interest in quasi-two-dimensional trapped BECs. Here, *r*=|*x*| represents the radial variable, and the density normalization condition reads
*N*, conservation of mass then yields the following continuity equation:

Assuming in this large *N* limit a radially symmetric density, we note that for any smooth radial function *g*(*r*), we have the following identity:
*ρ*(*y*)=*ρ*(|*y*|)),
*V* (*r*) be the bracketed expression in equation (4.3) so that *v*=*V* (*r*)*x* and note that
*ρ*_{t}+∇⋅(*vρ*)=0 becomes
*R*. Then we have
*r*=*R* then evolves according to
*f*(*R*)=*a*/(1−*R*^{2}), we obtain the equation for the support radius,

From equation (4.6), at the steady state we have *u*=−*r*^{2}( *f*(*r*)−*ω*)/(2*πc*). From equation (4.7), we then obtain *ρ*=−(1/2*πc*)(1/*r*)(∂/∂*r*)(*r*^{2}[*f*(*r*)−*ω*]) or
*ω*, equation (4.9) has either zero or (one at the critical point or) two solutions for *R*. If it has two solutions *R*_{−}<*R*_{+}, then *R*_{−} is stable and *R*_{+} is unstable, as is easily deduced from equation (4.8). The threshold occurs by setting ∂*ω*/∂*R*=0 to obtain
*ω*>*ω*_{c}: the one with smaller support is stable and the one with bigger support is unstable. Interestingly, the density *ρ*(*R*) vanishes precisely at *R*=*R*_{c}. Hence, among the two solutions only the stable one with *R*_{−}<*R*_{c} is physically relevant, while the unstable one with *R*_{+}>*R*_{c} cannot be computationally obtained (and presumably physically observed as it would involve negative vortex densities for *R*>*R*_{c}). In order to compute the steady-state distribution of the vortices, we first evolved equation (2.4) until it settled to an equilibrium state (figure 1*a*). We then computed the Voronoi tessellation of the plane using the Matlab function `voronoi` (figure 1*b*). This tessellation assigns to each vortex *x*_{j} a region (Voronoi cell) which consists of all points in the plane that are closer to *x*_{j} than to any other vortex. We then approximated the density distribution at *x*_{j} by *ρ*(*x*_{j})=1/*area*_{j}, where *area*_{j} is the area of the Voronoi cell associated to *x*_{j}. The resulting distribution (extrapolated linearly between the points) is plotted in figure 1*c*. Finally, in figure 1*d* we plot the radial density *ρ*(*r*), which we computed by taking the average of *ρ*(*x*) along |*x*|=*r*. We compare this to the distribution from the analytical expression of equation (4.10). There is an excellent agreement between the two corroborating the value of our theoretical prediction. Figure 2 shows that this agreement persists for smaller values of *N* (e.g. *N*=25) as well, although naturally it becomes progressively worse as *N* is decreasing.

In a future work, we plan to study the stability of the steady state (4.10). Numerical computations of equation (2.4) show that solution (4.10) is indeed a stable equilibrium for the aggregation model. The corresponding relative equilibrium of the BEC model (2.1) is then neutrally stable and has vibrational or the so-called Tkachenko modes [30,51,52]. We plan to extend the techniques in this section to compute the vibrational modes in the continuum limit of large *N*.

## 5. Conclusion and future challenges

In summary, in this work we have revisited two opposite limits of the quasi-two-dimensional co-rotating vortex dynamics in BECs. Motivated by the recent success of particle models in capturing experimental features of both the counter- and co-rotating vortex case, we have attempted to examine in detail (fully analytically, wherever possible) both the small *N* and the large *N* limit of such *N*-vortex configurations. In the former case, we obtained vortex configurations in the form of polygonal rings. We generalized the classical result of Havelock [10] unveiling that the ring is generically unstable for *N*≥7 in the case of monotonic precessional frequency dependence on the distance from the trap centre. Moreover, we showed that the critical contribution of the precessional term creates the potential for stable asymmetric, as well as other configurations even for *N*=2,…,6, for sufficiently high angular momentum. In that light, we also mentioned in passing the *N*+1 vortex configuration, whose stability is analysed in appendix A(b). The opposite limit of large *N* is quite interesting in its own right. As polygonal configurations are already highly unstable for sufficiently large *N*, a fundamentally different distribution is expected for large *N*. This distribution was identified in a radial form, by looking at the corresponding continuum equation and was corroborated numerically. An ongoing collaboration with the group that has made critical earlier experimental contributions in this theme [25,38] suggests the feasibility of looking at controllably small numbers of *N* (up to 11) as well as at the regime of large *N* regime experimentally.

The results in §4 generalize a similar computation for classical vortex dynamics [44]. The methods used here and in [44] are borrowed from the literature on biological swarming (e.g. [49,50]). Similar techniques were also recently used to study predator–swarm dynamics [53]. It is hoped that further developments in the swarming literature will help to shed light on the behaviour of BEC (and in particular, the stability of vortex lattices) and vice versa.

There are numerous directions in which we foresee that this activity can be extended. On the one hand, it would be particularly interesting (as the experimental possibilities reported in [38] allow the ‘dialling in’ of different numbers of vortices, e.g. between 1 and 11) to explore further the case of intermediate-size clusters i.e. between *N*=5 and *N*=11. There, identifying the potential *N*-vortex ring polygons, or that of *N*+1 rings or the examination of different ground-state configurations would be relevant to examine. On the other hand, in the case of large *N*, our preliminary computations (via fixed point iterations of a Newton scheme) reveal a large number of excited-state configurations. It will be interesting to explore in future studies whether these are generically unstable or whether additional dynamically stable large *N* limits could, in principle be accessible as well. Furthermore, examinations of multi-component (e.g. pseudo-spinor) settings in [30,31], of potentials of different symmetry (such as square optical lattices, which can again induce structural phase transitions [54]) motivate analogous considerations/extensions at the level of our particle model.

While the mean-field theory is successful at predicting the large-scale vortex density distribution, it does not capture the fine structure of the BEC lattice itself (e.g. [55]) where different lattices and where their dynamics and internal (Tkachenko) modes were observed [56]. The point-vortex BEC model (2.1) is an approximation to the full system more accurately described by a three-dimensional Gross–Pitaevskii model, while neglecting the vortex core structure. The three-dimensionality can lead to more complex configurations such as ‘Olympic rings’ (e.g. [56,57]). It would be interesting to see whether the techniques of this paper could also be applied to such configurations. Finally, the examination of trapped, interacting vortex rings in three dimensions both in the context of few [58–61] and in that of many such rings would be a broad direction of considerable importance for future studies.

## Funding statement

T.K. was supported by NSERC grant no. 47050. P.G.K. and R.C.G. gratefully acknowledge the support of NSF-DMS-0806762, NSF-DMS-1312856 and NSF-CMMI-1000337.

## Acknowledgements

We thank the anonymous referees for their valuable comments, which significantly improved the paper.

## Appendix A

### (a) Relationship between stability of aggregation and Bose–Einstein condensate equations

We prove theorem 2.1 here. Linearize equation (2.4) around the steady-state *x*_{j}(*t*)=*ξ*_{j} by using *x*_{j}(*t*)=*ξ*_{j}+*η*_{j}(*t*) with |*η*_{j}|≪1. We then obtain the system
*η*=(*η*_{1}…*η*_{N})^{T} is the perturbation vector; overbar denotes the complex conjugate; *L* is a symmetric complex matrix whose entries are
*D* is a diagonal real matrix whose entries are
*N* ODEs given by
*z*_{j}(*t*)=*e*^{iωt}*ξ*_{j} of equation (2.1), we find that its eigenvalues are given by the matrix
*A* are strictly negative, then all eigenvalues of *B* are purely imaginary. As *A* is Hermitian, we may write *E* is a diagonal matrix whose diagonal entries are the eigenvalues of *A*, and *U* is unitary. Assume that all eigenvalues of *A* are negative. Then we can write *E*=−*Q*^{2}, where *Q* is a real diagonal matrix, so that *M*_{1}*M*_{2} and *M*_{2}*M*_{1} is the same, so that *B* has the same spectrum as the matrix

To show the converse, note that under either condition (i) or condition (ii) of theorem 2.1, the matrix *D* given by equation (A 2) is a multiple of identity so we may write *D*=*dI*, where *d*=*f*′(*R*)*R*/2−*c*(*N*−1)/*R*^{2} is a negative constant. In this case, the eigenvalues of *A* are given by *B* are given by _{B} is purely imaginary if and only if *ϵ*<*d*^{2}, which is if and only if

### (b) *N*+1 state: ring solution with a vortex at the centre

Here, we prove theorem 3.2. Similar to the ring steady state, we consider the relative equilibrium of the aggregation model with *N*+1 vortices; *N* on the ring and one at the centre. As in [47], we will actually consider a slightly more general problem where the central vortex has weight *b*, whereas other vortices have weight *c*; theorem 3.2 will follow by taking *b*=*c*. The starting point is
*R* satisfies
*b*=*c* recovers formula (3.14).

Next, we consider perturbations to the *N*+1 vortex configurations. As before, we perturb the steady state as

*Subspace 1:* use the ansatz
*e*^{im2πj/N} and *e*^{−im2πj/N}, system (3.4) decouples into a sequence of 2×2 subproblems
_{−}(*m*)≤0 for all *m*. As in theorem 3.1, this expression is maximized when *m*=⌊*N*/2⌋. Setting *b*=*c* recovers equation (3.15).

*Subspace 2:* we use the ansatz
*b*=*c* into the matrix above yields equation (3.16).

## Footnotes

↵1 It is worth noting here that while this is the standard form for the angular momentum in the context of point vortices (cf. eq. (2.13) [16]) this is not identical to the angular momentum of the BEC, which is well-established, e.g. in the context of the corresponding PDE model of the Gross–Pitaevskii equation [1–3]. In our considerations here, we will use the term for the former and not the one for the latter.

- Received January 16, 2014.
- Accepted May 7, 2014.

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