A Appendix

# A Tale of Two Distributions: From Few To Many Vortices In Quasi-Two-Dimensional Bose-Einstein Condensates

## 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 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 of Havelock [Phil. Mag. 11, 617 (1931)] illustrating that vortex polygons in the form of a ring are unstable for . Additionally, we reconcile this modification with the recent identification of symmetry breaking bifurcations for the cases of . We also briefly discuss the case of a ring of vortices surrounding a central vortex (so-called configuration). We finally examine the opposite limit of large and illustrate how a coarse-graining, continuum approach enables the accurate identification of the radial distribution of vortices in that limit.

Vortex dynamics, Bose-Einstein condensates, vortex lattices

## I 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); (2); (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 He (6), electron columns confined in Malmberg-Penning traps (7) and even magnetized, millimeter sized disks rotating at a liquid-air interface (8). 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 (9), the consideration of higher numbers of vortices, e.g., in Ref. (10) and even of asymmetric equilibria thereof, e.g., in Ref. (11). Much of this activity has been summarized, e.g., in the review (12) and the book (13). A more recent exposition tying the classical fluid point vortex problem of Refs. (12); (13) with the BEC-oriented motivation of the present work is given in Ref. (14).

In the context of BECs, vortices have been one of the focal themes of numerous investigations which have now been summarized in many reviews (15); (16); (17); (18); (19). A considerable volume of associated experimental efforts was concerned with the creation schemes of such structures (20); (21); (22); (23); (24), especially in the form of unit charge vortices. However, there were also efforts to produce the (dynamically unstable) vortices of higher topological charge (25) and observe their decay. Furthermore, numerous studies focused on vortex lattices with a large number of vortices (26); (27); (28). 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 Ref. (29) (and their theoretical examination in Ref. (30)), it has recently ignited considerable interest due to the systematic formation and observation of counter-circulating vortex dipoles (31); (23); (32); (33), tripoles (34) and also sets of 2-, 3-, 4- (and more generally controllably few-) vortices (35).

One of the features that are especially interesting in connection to this BEC context concerns the relevance and usefulness of the modeling 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 (32); (33), but also as a tool to unveil subtle bifurcation and symmetry breaking phenomena as in the case of few co-rotating vortices (35) (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 center), or the destabilization of vortex ring formations for small vortex number .

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 (23); (32); (35); (33) 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; see, e.g., Ref. (36) 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 (37). 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 Ref. (37) (and references therein).

The present 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 setups create the vortices through the imparting of angular momentum to the system (17). 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 . 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 Ref. (9) (see also Ref. (38)). 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 in the ring formation, here due to the precessional term, even the 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 to , whereas the work of Ref. (9) 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 Ref. (35).

• 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 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 configuration, where vortices form a ring, while is located at the trap center. 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 Refs. (31); (23); (32); (33) or Ref. (35) 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 Sec. II, we offer the general theoretical setup in line with the earlier works (such as Refs. (32); (35)) that corroborated its correspondence with experimental results. In Sec. III, we present the analysis of the polygonal vortex ring and how it connects to the stability conclusions of recent few vortex studies such as Refs. (35); (39). In Sec. IV, we explore the opposite limit of large and the corresponding coarse-grained, continuum model. Again a comparison is offered, this time with numerical computations within that limit. Finally, in Sec. V, we summarize our findings and present a number of conclusions and possible future directions. The appendices present a number of technical details associated, e.g., with the stability of the vortex configuration.

## Ii Theoretical Setup

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 , where denotes the planar position of the -th vortex. This equation reads:

 ˙zj=if(∣∣zj∣∣)zj+ic∑k≠jzj−zk∣∣zj−zk∣∣2,   j=1…N. (1)

In the right hand side of Eq. (1) the first term represents the precession of the vortex around the trap center. Here we consider only vortices of a unit charge (given their dynamical stability) and assume that all have the same charge (without loss of generality we assume this to be ), as in the context of the recent experiments of Ref. (35). 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 Refs. (32); (35); see also the detailed analysis/comparison with single vortex precession experiments in Ref. (23). In particular, in line with these works, we will choose:

 f(r)=a1−r2. (2)

Here, the radius has been normalized to the Thomas-Fermi radius (roughly tantamount to the spatial extent of the BEC), while the factor representing the precession frequency of the vortex very near the center of the trap can be absorbed into a rescaling of time (rendering time dimensionless).

On the other hand, the second term in Eq. (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 Ref. (35), where it was explained that a “typical” experimentally relevant value for is . 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, e.g., in Ref. (40). 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 Eq. (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 . These satisfy the equation

 ωxj=f(∣∣xj∣∣)xj+c∑k≠jxj−xk∣∣xj−xk∣∣2, (3)

which will be the main focus of our considerations in what follows.

As a (partially numerical) aside, we will also consider the following “aggregation” equation:

 ˙xj=(f(∣∣xj∣∣)−ω)xj+c∑k≠jxj−xk∣∣xj−xk∣∣2. (4)

By construction, the relative equilibrium of Eq. (1) corresponds to the relaxational equilibrium of Eq. (4) and vice-versa. Moreover there is an intimate connection between the stability of the two models. In Appendix A.1 we prove the following result.

###### Theorem II.1

Suppose that an equilibrium of Eq. (4) is stable. Then the relative equilibrium of Eq. (1) is (neutrally) stable. The converse is also true in the following two cases: (i) or (ii)  for all

The case (i) of this theorem was shown in Ref. (41); this is reproduced in Appendix A.1. However the proof in Ref. (41) does not work for a general case of non-homogeneous , and a more general approach is taken here. On the flip side, we can only show the stability of Eq. (4) implies stability of Eq. (1); we do not know how to prove the converse (nor do we have counter-examples).

## Iii Ring solutions

The prototypical configuration that one expects to identify as a stable equilibrium in the case of small vortex number (motivated by the corresponding result in the absence of precession (9); (10)) is the “ring” configuration with the vortices sitting at the vertices of a canonical polygon. Assuming such a relative equilibrium to Eq. (1) with radius , we can write it in our complex notation as:

 zj(t)=Rexp(iωt)exp(2πiNj). (5)

In this case, we can compute:

 ∑k≠jzj−zk∣∣zj−zk∣∣2=1Rexp(iωt+2πiNj)N−12, (6)

where we have used the identity

 N−1∑k=111−exp(−2πiNk)=N−12.

 ω=f(R)+c(N−1)2R2. (7)

As the case (ii) of Theorem II.1 [of Appendix A.1] shows, the ring for the vortex model Eq. (1) is stable if and only if it is stable for the aggregation model Eq. (4). The full characterization of linear stability is given by the following theorem (following the procedure used, e.g., in Refs. (42); (43)).

###### Theorem III.1

Consider the ring solution for Eq. (1), of radius as given by Eq. (5), where the frequency is given by Eq. (7). Suppose that is odd. Then the ring is stable provided that

 f′(R)R+c8R2(N−1)(N−7)<0,

and is unstable if the inequality is reversed. Suppose that is even. Then the ring is stable provided that

 f′(R)R+c8R2(N2−8N+8)<0,

and is unstable if the inequality is reversed.

The ring is generically unstable if and is an increasing function.

Proof. Because of Theorem II.1 case (ii), it is sufficient to consider the stability of the steady state of the aggregation equation (4). Consider small perturbations of this state, of the form

 xk(t)=Rexp(2πikN)(1+hk(t)),   |hk|≪1,

where is a small, complex-valued, perturbation. After some algebra we obtain for the evolution of the small perturbations

 dhjdt=(f′(R)R2+f(R)−ω)hj+f′(R)R2¯hj+c∑k≠jexp(2π(k−j)/N)¯hj−¯hk4R2sin2(π(k−j)N), (8)

where the overbar denotes complex conjugation. Using the following Fourier mode decomposition for the perturbation:

 hj(t)=ξ+(t)eim2πj/N+ξ−(t)e−im2πj/N,   m∈N, (9)

and collecting like terms in and the system (8) decouples into a sequence of subproblems

 ξ′+ =(f′(R)R2+f(R)−ω)ξ++f′(R)R2¯ξ−+σ+¯ξ−, (10) ξ′− =(f′(R)R2+f(R)−ω)ξ−+f′(R)R2¯ξ++σ−¯ξ+, (11)

where

 σ+ ≡c∑k,k≠jei2π(k−j)/N−eim2π(k−j)/N4R2sin2(π(k−j)/N), σ− ≡c∑k,k≠jei2π(k−j)/N−e−im2π(k−j)/N4R2sin2(π(k−j)/N).

Using the following identity:

 Missing or unrecognized delimiter for \left (12)

it is possible to write

 σ≡σ+=σ−=c2R2(m−1)(N−m−1). (13)

Taking the conjugate of Eq. (11), the system can be written as

 Unknown environment 'array% (14)

whose eigenvalues are given by

 λ±(m)=f′(R)R2+f(R)−ω±(f′(R)R2+c2R2(m−1)(N−m−1)),  m=0…N−1.

Using Eq. (7) this simplifies to

 λ−(m)=c2R2[−(N−1)−(m−1)(N−m−1)],λ+(m)=f′(R)R+c2R2[(m−1)(N−m−1)−(N−1)]. (15)

Note that ; this mode corresponds to rotation invariance. Moreover is positive for all integers and attains its max at Therefore correspond to stable eigendirections and the instability threshold is obtained by setting , where denotes the floor function.

Suppose that is odd. Then we substitute into Eq. (15) to obtain

 λ+((N−1)/2)=f′(R)R+c8R2(N−1)(N−7),   N odd.

If is even we substitute which yields

 λ+(N/2)=f′(R)R+c8R2(N2−8N+8),   N even.

The ring is stable for the aggregation equation (and by extension for the actual vortex model) if the right hand side of the above two equations (respectively, for odd and even) is negative and unstable otherwise. If this recovers the threshold.

Some observations are important to make here. In particular, we wish to examine the special cases of 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 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 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 , we can express the corresponding eigenvalue in the general form , where for odd and for even. A key observation now (which directly connects the present work with the experimental observations of Ref. (35) and the computational/theoretical analysis of Ref. (39)) is that although for these cases the configuration is not generically unstable, nevertheless it may be unstable for sufficiently large . 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 and in the case of the ring configuration acquires the especially simple form .1 In that light and taking as in Eq. (2) with , the eigenvalue whose zero crossing will determine the potential instability of a configuration with reads:

 λ=2LN(1−LN)2+cN8L(N2−8N+s).

It is then straightforward to infer that the critical angular momentum given by setting satisfies

 16L2+c(N−L)2(N2−8N+s)=0

which yields a critical of the form:

 Lcr=cN(N2−8N+s)+4N√c(8N−N2−s)16+c(N2−8N+s).

Remarkably, this expression yields the relevant critical angular momenta for all cases of to in direct agreement with the expressions given in Refs. (35); (39) for the symmetry breaking bifurcations due to the destabilization of the ring configuration. Using the notation

 r21≡√c√c+2,r22≡√c√c+√2, (16)

we have that:

 Lcr,N=2=2r21, Lcr,N=3=3r22, Lcr,N=4=4r22, (17) Lcr,N=5=5r22, Lcr,N=6=6r21.

It is rather intriguing that the critical angular momentum for the different cases exhibits an apparent pattern although not a clearly definite one. Interestingly, even for , while the dominant eigenvalue is always positive (as indicated above), even the next one can be seen to cross at , extending this interesting pattern. We also note in passing that for the case of , the relevant critical point had not been previously identified in Ref. (35) or Ref. (39).

However, a final observation is also in order. As discussed in Ref. (39), the interval (of ) of dynamical stability of the small configurations does not coincide with the interval of for which these constitute the ground state of the system. In the case of , the stability threshold and ground state asymmetry threshold coincide (this is a supercritical pitchfork point), but in other cases such as and , asymmetric configurations (such as isosceles triangles for or rhombic configurations for ) acquire lower energy than the polygonal ring distinctly before its loss of stability threshold (39). 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 Ref. (11)), 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 examined in Ref. (39) (and obviously also for larger ) is the so-called vortex configuration, consisting of vortices on the polygonal ring and one more at the center. For the classical vortex problem [ in Eq. (1)], the stability of this configuration was analyzed in Ref. (44); see also Ref. (45). In Appendix A.2 we show the following generalization to the full BEC problem:

###### Theorem III.2

Consider the configuration of Eq. (1) consisting of vortices uniformly distributed on a ring of radius and angular velocity given by Eq. (5) plus a vortex at the origin. Then satisfies

 ω=f(R)+c2R2(N+1). (18)

Let

 λ∗+(m)=λ+(m)−2c/R2 (19)

where as given by Eq. (15) and let

 Missing or unrecognized delimiter for \left (20)

The configuration is stable if and only if   and all eigenvalues of are negative.

In particular, this theorem shows that the configuration is stable if the -ring is stable and if in addition the matrix is negative definite.

When both and configurations are marginally stable [the former due to when the latter due to the eigenvalue of crossing zero when ]; the configuration is stable for and is unstable for or This is in agreement with the results in Ref. (44). When is increasing as in Eq. (2), both and configurations lose their marginal stability and become unstable, so that configuration becomes unstable for any or .

## Iv Continuum limit

Having considered the case of small , we now explore the opposite limit of large . Notice, however, that our results for the ring and the vortex configuration are entirely general and the dynamical stability thereof is obtained for any . Yet, these configurations can only be stable when is sufficiently small, as discussed above, e.g., for the ring state (and even then for sufficiently small ). Hence, in the opposite limit of large , we expect that a substantially different vortex distribution will arise. This expectation is confirmed not only by the well-known vortex lattice observations of, e.g., Refs. (26); (27); (28), but also by the direct numerical evolution of the aggregation equation (4); see the top left panel of Fig. 1. Recall that the aggregation equation has the benefit of relaxing to equilibrium attractors corresponding to the marginally stable equilibria of our original Eq. (1). It is then particularly relevant to attempt to identify the distribution of such a large number of vortices . Here, we will use techniques similar to Ref. (46) to derive the limiting density profile.

Following, e.g., the discussion of Ref. (47), we coarse-grain by defining the particle density to be

 ρ(x)=∑k=1…Nδ(x−xk), (21)

where is the Dirac-delta function. It is then straightforward to rewrite the aggregation Eq. (4) as where we define the continuum limit of the velocity as

 v(x)≡(f(r)−ω)x+c∫R2x−y|x−y|2ρ(y)dy.

The term is relevant, in particular, to the precessional dynamics of interest in quasi-two-dimensional trapped BECs. Here represents the radial variable and the density normalization condition reads

 ∫R2ρ(x)dx=N.

In the limit of large , conservation of mass then yields the following continuity equation:

 ρt+∇⋅(vρ)=0. (22)

Assuming in this large limit a radially symmetric density, we note that for any smooth radial function we have the following identity:

 ∫R2x−y|x−y|2g(|y|)dy=x2πr2∫r0g(s)sds,

so that [with slight abuse of notation ],

 v=(f(r)−ω+2πcr2∫r0ρ(s)sds)x. (23)

Let be the bracketed expression in Eq. (23) so that and note that

 ∇⋅(vρ)=∇⋅(V(r)ρx)=1r(Vρr2)r

so that becomes

 (rρ)t+(Vrρr)r=0. (24)

Now define

 u(r)=∫r0ρ(s)sds.

Integrating Eq. (24) we obtain

 ut+Vrur=0. (25)

Recall that we have

 Vr =r(f(r)−ω)+2πcr∫r0ρ(s)sds =r(f(r)−ω)+2πcru.

Thus we obtain the following characteristics for Eq. (25)

 drdt =r(f(r)−ω)+2πcru, dudt =0. (26)

Now suppose that the initial density is radially symmetric and has finite support of radius Then we have:

 u(R)=∫R0rρdr=N2π; (27)

the corresponding characteristic then evolves according to

 dRdt=R(f(R)−ω)+cNR. (28)

In particular at the steady state and with we obtain the equation for the support radius,

 ω=a1−R2+cNR2. (29)

From Eq. (26), at the steady state we have . From Eq. (27) we then obtain or

 ρ=1πc(ω−a(1−r2)2). (30)

Note that for fixed , Eq. (29) has either zero or (one at the critical point or) two solutions for If it has two solutions then is stable and is unstable, as is easily deduced from Eq. (28). The threshold occurs by setting to obtain

 ωc=(√a+√cN)2,R2c=√cN√a+√cN.

Thus two solutions exist when the one with smaller support is stable and the one with bigger support is unstable. Interestingly, the density vanishes precisely at . Hence, among the two solutions only the stable one with is physically relevant, while the unstable one with cannot be computationally obtained (and presumably physically observed since it would involve negative vortex densities for ). In order to compute the steady-state distribution of the vortices, we first evolved Eq. (4) until it settled to an equilibrium state [Fig. 1(a)]. We then computed the Voronoi tessellation of the plane using the Matlab function voronoi [Fig. 1(b)]. This tessellation assigns to each vortex a region (Voronoi cell) which consists of all points in the plane that are closer to than to any other vortex. We then approximated the density distribution at by where is the area of the Voronoi cell associated to The resulting distribution (extrapolated linearly between the points) is plotted in Fig. 1(c). Finally, in Fig. 1(d) we plot the radial density which we computed by taking the average of along . We compare this to the distribution from the analytical expression of Eq. (30). 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 (e.g., ) as well, although naturally it becomes progressively worse as is decreasing.

In a future work, we plan to study the stability of the steady state (30). Numerical computations of Eq. (4) show that solution (30) is indeed a stable equilibrium for the aggregation model. The corresponding relative equilibrium of the BEC model (1) is then neutrally stable and has vibrational or the so-called Tkachenko modes (27); (48); (49). We plan to extend the techniques in this section to compute the vibrational modes in the continuum limit of large .

## V Conclusions and Future Challenges

In summary, in the present work we have revisited two opposite limits of the quasi-two-dimensional co-rotating vortex dynamics in Bose-Einstein condensates. 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 and the large limit of such -vortex configurations. In the former case, we obtained vortex configurations in the form of polygonal rings. We generalized the classical result of Ref. (9) unveiling that the ring is generically unstable for in the case of monotonic precessional frequency dependence on the distance from the trap center. Moreover, we showed that the critical contribution of the precessional term creates the potential for stable asymmetric, as well as other configurations even for , for sufficiently high angular momentum. In that light, we also mentioned in passing the vortex configuration, whose stability is analyzed in Appendix A.2. The opposite limit of large is quite interesting in its own right. Since polygonal configurations are already highly unstable for sufficiently large , a fundamentally different distribution is expected for large . 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 (see Refs. (23); (35)) suggests the feasibility of looking at controllably small numbers of (up to ) as well as at the regime of large regime experimentally.

The results in Sec. IV generalize a similar computation for classical vortex dynamics (41). The methods used here and in Ref. (41) are borrowed from the literature on biological swarming, see, e.g., Refs. (46); (47). Similar techniques were also recently used to study predator-swarm dynamics (50). It is hoped that further developments in the swarming literature will help to shed light on the behavior 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 (since the experimental possibilities reported in Ref. (35) allow the “dialing in” of different numbers of vortices, e.g., between 1 and 11) to explore further the case of intermediate-size clusters i.e., between and . There, identifying the potential -vortex ring polygons, or that of rings or the examination of different ground state configurations would be relevant to examine. On the other hand, in the case of large , 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 limits could, in principle be accessible as well. Furthermore, examinations of multi-component (e.g., pseudo-spinor) settings in Refs. (27); (28), of potentials of different symmetry (such as square optical lattices, which can again induce structural phase transitions (51)) 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; see, e.g., Ref. (52) where different lattices and where their dynamics and internal (Tkachenko) modes were observed (53). The point vortex BEC model (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”, see, e.g., Refs. (53); (54). It would be interesting to see if 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 (55); (56); (57); (58) and in that of many such rings would be a broad direction of considerable importance for future studies.

## Vi Acknowledgements

We thank the anonymous referees for their valuable comments, which significantly improved the paper. T.K. was supported by NSERC grant no. 47050. P.G.K. and R.C.G. gratefully acknowledge the support of NSF-DMS-0806762 and NSF-DMS-1312856 and NSF-CMMI-1000337.

## Appendix A Appendix

### a.1 Relation between stability of aggregation and BEC equations

We prove Theorem II.1 here. Linearize Eq. (4) around the steady state by using with We then obtain the system

 ˙η=Dη+L¯η. (31)

Here is the perturbation vector; overbar denotes the complex conjugate; is a symmetric complex matrix whose entries are

 Ljk=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩c(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ξj−ξk)2,     j≠k,f′(∣∣ξj∣∣)ξ2j2∣∣ξj∣∣−∑k≠jc(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ξj−ξk)2,       j=k,

and is a diagonal real matrix whose entries are

 Djj=f(∣∣ξj∣∣)+f′(∣∣ξj∣∣)∣∣ξj∣∣2−ω. (32)

By taking the complex conjugate of Eq. (31) we obtain a closed system of ODEs given by

 {∂tη=Dη+L¯η,∂t¯η=D¯η+¯Lη. (33)

Linearizing around the steady state equilibrium, we find that the eigenvalues of the zero equilibrium of Eq. (33) are given by the matrix

 A=[DL¯LD].

Performing a similar analysis the relative equilibrium of Eq. (1), we find that its eigenvalues are given by the matrix

 B=iJA=[iDiL−i¯L−iD],

where

 J=[I00−I].

Next, we show that if all eigenvalues of are strictly negative, then all eigenvalues of are purely imaginary. Since is Hermitian, we may write where is a diagonal matrix whose diagonal entries are the eigenvalues of and is unitary. Assume that all eigenvalues of are negative. Then we can write where is a real diagonal matrix, so that . Note that in general, the spectrum of matrices and is the same, so that has the same spectrum as the matrix whose eigenvalues are purely imaginary since is Hermitian.

To show the converse, note that under either conditions (i) or (ii) of Theorem II.1, the matrix given by Eq. (32) is a multiple of identity so we may write where is a constant. In this case, the eigenvalues of are given by where is an eigenvalue of whereas the eigenvalues of are given by It follows that is purely imaginary if and only if which is if and only if

### a.2 N+1 state: ring solution with a vortex at the center

Here we prove Theorem III.2. Similar to the ring steady state, we consider the relative equilibrium of the aggregation model with vortices;  on the ring and one at the center. As in Ref. (44), we will actually consider a slightly more general problem where the central vortex has weight whereas other vortices have weight ; Theorem III.2 will follow by taking The starting point is

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩˙zj=(f(∣∣zj∣∣)−ω)zj+c∑k≠jzj−zk∣∣zj−zk∣∣2+bzj−zN+1∣∣zj−zN+1∣∣2,   j=1…N,˙zN+1=(f(|zN+1|)−ω)zN+1+cN∑k≠jzN+1−zk|zN+1−zk|2.

As before, we make the ansatz

 ⎧⎨⎩zj(t)=Rexp(2πiNj),  j=1…N,zN+1(t)=0.

Then satisfies

 ω=f(R)+c(N−1)2R2+bR2; (34)

setting recovers formula (18).

Next we consider perturbations to the vortex configurations. As before we perturb the steady state as

 xk(t) =Rξk(1+hk(t)),   |hk|≪1,  k=1…N xN+1(t) =RhN+1(t),

where we defined

 ξ≡exp(2πi/N),

to obtain

 dhjdt=(f′(R)R2+f(R)−ω)hj+f′(R)R2¯hj+cN∑k≠jξk−j¯hj−¯hk4R2sin2(π(k−j)N)−b¯hj−ξj¯hN+1R2, dhN+1dt=(f(0)−ω)hN+1+cN∑k=1ξk¯hkR2.

The solution decomposes into a product of two subspaces:

Subspace 1: Use the ansatz

 hj(t)=ξ+(t)ξmj+ξ−(t)ξ−mj, m∈N,  j=1…N;      hN+1=0,

and collecting like terms in and the system (8) decouples into a sequence of subproblems

 ξ′+= (f′(R)R2+f(R)−ω)ξ++(f′(R)R2−bR2)¯ξ−+¯ξ−c∑k,k≠jξk−j−ξm(k−j)4R2sin2(π(k−j)/N), ξ′−= (f′(R)R2+f(R)−ω)ξ−+(f′(R)R2−bR2)¯ξ++¯ξ+c∑k,k≠jξk−j−ξ−m(k−j)4R2sin2(π(k−j)/N),

and, as previously, we obtain

 λ±(m)=f′(R)R2+f(R)−ω±(f′(R)R2−bR2+c2R2(m−1)(N−m−1)),  m=0…N−1.

Using Eq. (34) yields

 λ+(m)=f′(R)R+c2R2{(m−1)(N−m−1)−(N−1)}−2bR2,

with for all As in Theorem III.1, this expression is maximized when . Setting recovers Eq. (19).

Subspace 2: we use the ansatz

 hj(t)=ξ+(t)ξj+ξ−(t)ξ−j,   j=1…N;      hN+1=η(t),

which yields

 ξ′+ =(f′(R)R2+f(R)−ω)ξ++(f′(R)R2−bR2)¯ξ−+b¯ηR2, ξ′− =(f′(R)R2+f(R)−ω)ξ−+(f′(R)R2−bR2)¯ξ+, dηdt =(f(0)−ω)η+cN∑k=1¯ξ+R2,

or

 Unknown environment '%

Substituting into the matrix above yields Eq. (20).

### 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) in p. 514 of the review (14), this is not identical to the angular momentum of the Bose-Einstein condensate, which is well-established, e.g., in the context of the corresponding PDE model of the Gross-Pitaevskii equation (1); (2); (3). In our considerations here, we will use the term for the former and not the one for the latter.

### References

1. C.J. Pethick and H. Smith, Bose-Einstein condensation in dilute gases, Cambridge University Press, Cambridge, 2002.
2. L.P. Pitaevskii and S. Stringari, Bose-Einstein Condensation, Oxford University Press, Oxford, 2003.
3. P.G. Kevrekidis, D.J. Frantzeskakis, and R. Carretero-González, Emergent Nonlinear Phenomena in Bose-Einstein Condensates, Springer-Verlag, Berlin, 2008.
4. W.T. Kelvin, Mathematical and Physical Papers, Cambridge University Press, 1878, Vol. IV, p. 135.
5. L. Onsager, Statistical Hydrodynamics, Il Nuovo Cimento 6, 279 (1949).
6. E. Yarmchuk, M. Gordon, and R. Packard, Observation of stationary vortex arrays in rotating superfluid helium, Phys. Rev. Lett. 43, 214 (1979).
7. D. Durkin and J. Fajans, Experiments on two-dimensional vortex patterns, Phys. Fluids 12, 289 (2000).
8. B. Grzybowski, H. Stone, and G. Whitesides, Dynamic self-assembly of magnetized, millimetre-sized objects rotating at a liquid-air interface, Nature 405, 1033 (2000); B. Grzybowski, H. Stone, and G. Whitesides, Dynamics of self assembly of magnetized disks rotating at the liquid-air interface, Proc. Natl. Acad. Sci. 99, 4147 (2002).
9. T. Havelock, The stability of motion of rectilinear vortices in ring formation, Philos. Mag. 11, 617 (1931).
10. L.J. Campbell and R. Ziff, A catalog of two-dimensional vortex patterns, Los Alamos Scientific Laboratory Report No. LA-7384-MS (1978); L.J. Campbell and R.M. Ziff Vortex patterns and energies in a rotating superfluid, Phys. Rev. B 20, 1886 (1979).
11. H. Aref and D. Vainchtein, Point vortices exhibit asymmetric equilibria, Nature 392, 769 (1998).
12. H. Aref, P. Newton, M. Stremler, T. Tokieda, and D. Vainchtein, Vortex crystals, Adv. Appl. Mech. 39, 1 (2003).
13. P. Newton, The -vortex problem: analytical techniques, Springer-Verlag, New York, 2001.
14. P.K. Newton, G. Chamoun, Vortex lattice theory: a particle interaction perspective, SIAM Rev. 51, 501 (2009).
15. R.J. Donnelly, Quantized Vortices in Helium II, Cambridge University Press, New York, 1991.
16. A.L. Fetter and A.A. Svidzinsky, Vortices in a trapped dilute Bose-Einstein condensate, J. Phys.: Condens. Matter 13, R135 (2001).
17. A.L. Fetter, Rotating trapped Bose-Einstein condensates, Rev. Mod. Phys. 81, 647 (2009).
18. P.G. Kevrekidis, R. Carretero-González, D.J. Frantzeskakis, and I.G. Kevrekidis, Vortices in Bose-Einstein Condensates: Some Recent Developments, Mod. Phys. Lett. B 18, 1481 (2004).
19. M. Tsubota, K. Kasamatsu, and M. Kobayashi, Quantized vortices in superfluid helium and atomic Bose-Einstein condensates, arXiv:1004.5458.
20. J.E. Williams and M.J. Holland, Preparing topological states of a Bose-Einstein condensate, Nature 401, 568 (1999).
21. R. Onofrio, C. Raman, J.M. Vogels, J.R. Abo-Shaeer, A.P. Chikkatur, and W. Ketterle, Observation of Superfluid Flow in a Bose-Einstein Condensed Gas, Phys. Rev. Lett. 85, 2228 (2000).
22. C. Weiler, T. Neely, D. Scherer, A. Bradley, M. Davis, and B.P. Anderson, Spontaneous vortices in the formation of Bose-Einstein condensates, Nature 455, 948 (2008).
23. D.V. Freilich, D.M. Bianchi, A.M. Kaufman, T.K. Langin, and D.S. Hall, Real-Time Dynamics of Single Vortex Lines and Vortex Dipoles in a Bose-Einstein Condensate, Science 329, 1182 (2010).
24. D.R. Scherer, C.N. Weiler, T.W. Neely, and B.P. Anderson, Vortex Formation by Merging of Multiple Trapped Bose-Einstein Condensate, Phys. Rev. Lett. 98, 110402 (2010).
25. A.E. Leanhardt, A. Görlitz, A.P. Chikkatur, D. Kielpinski, Y. Shin, D.E. Pritchard, and W. Ketterle, Imprinting Vortices in a Bose-Einstein Condensate using Topological Phases, Phys. Rev. Lett. 89, 190403 (2002); Y. Shin, M. Saba, M. Vengalattore, T.A. Pasquini, C. Sanner, A.E. Leanhardt, M. Prentiss, D.E. Pritchard, and W. Ketterle, Dynamical Instability of a Doubly Quantized Vortex in a Bose-Einstein Condensate, Phys. Rev. Lett. 93, 160406 (2004).
26. C. Raman, J.R. Abo-Shaeer, J.M. Vogels, K. Xu, and W. Ketterle, Vortex Nucleation in a Stirred Bose-Einstein Condensate, Phys. Rev. Lett. 87, 210402 (2001).
27. I. Coddington, P. Engels, V. Schweikhard, and E.A. Cornell, Observation of Tkachenko Oscillations in Rapidly Rotating Bose-Einstein Condensates, Phys. Rev. Lett. 91, 100402 (2003);
28. V. Schweikhard, I. Coddington, P. Engels, S. Tung, and E.A. Cornell, Vortex-Lattice Dynamics in Rotating Spinor Bose-Einstein Condensates, Phys. Rev. Lett. 93, 210403 (2004).
29. K.W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard, Vortex Formation in a Stirred Bose-Einstein Condensate, Phys. Rev. Lett. 84 806 (2000).
30. Y. Castin and R. Dum, Bose-Einstein condensates with vortices in rotating traps, Eur. Phys. J. D 7, 399 (1999).
31. T.W. Neely, E.C. Samson, A.S. Bradley, M.J. Davis, and B.P. Anderson, Observation of Vortex Dipoles in an Oblate Bose-Einstein Condensate, Phys. Rev. Lett. 104, 160401 (2010).
32. S. Middelkamp, P.J. Torres, P.G. Kevrekidis, D.J. Frantzeskakis, R. Carretero-González, P. Schmelcher, D.V. Freilich, and D.S. Hall, Guiding-center dynamics of vortex dipoles in Bose-Einstein condensates, Phys. Rev. A 84, 011605 (2011).
33. P. Kuopanportti, J.A.M. Huhtamäki, and M. Möttönen Size and dynamics of vortex dipoles in dilute Bose-Einstein condensates, Phys. Rev. A 83, 011603(R) (2011).
34. J.A. Seman, E.A.L. Henn, M. Haque, R.F. Shiozaki, E.R.F. Ramos, M. Caracanhas, P. Castilho, C. Castelo Branco, P.E.S. Tavares, F.J. Poveda-Cuevas, G. Roati, K.M.F. Magalhaes, and V.S. Bagnato, Three-vortex configurations in trapped Bose-Einstein condensates, Phys. Rev. A 82, 033616 (2010).
35. R. Navarro, R. Carretero-González, P.J. Torres, P.G. Kevrekidis, D.J. Frantzeskakis, M.W. Ray, E. Altuntaş, and D.S. Hall, Dynamics of Few Co-rotating Vortices in Bose-Einstein Condensates, Phys. Rev. Lett. 110, 225301 (2013).
36. V. Koukouloyannis, G. Voyatzis, and P. G. Kevrekidis, Dynamics of three noncorotating vortices in Bose-Einstein condensates, Phys. Rev. E 89, 042905 (2014).
37. D. Yan, R. Carretero-González, D. J. Frantzeskakis, P. G. Kevrekidis, N. P. Proukakis, and D. Spirn, Exploring vortex dynamics in the presence of dissipation: Analytical and numerical results, Phys. Rev. A 89, 043613 (2014).
38. H. Cabral and D. Schmidt, Stability of relative equilibria in the problem of n+1 vortices. SIAM J. Math. Anal. 31, 231 (2000).
39. A.V. Zampetaki, R. Carretero-González, P.G. Kevrekidis, F.K. Diakonos, and D.J. Frantzeskakis, Exploring Rigidly Rotating Vortex Configurations and their Bifurcations in Atomic Bose-Einstein Condensate, Phys. Rev. E 88, 042914 (2013).
40. S. McEndoo and Th. Busch, Small numbers of vortices in anisotropic traps, Phys. Rev. A 79, 053616 (2009).
41. Yu. Chen, T. Kolokolnikov, and D. Zhirov, Collective behaviour of large number of vortices in the plane, Proc. Roy. Soc. A 469, 20130085 (2013).
42. T. Kolokolnikov, H. Sun, D. Uminsky, and A.L. Bertozzi, Stability of ring patterns arising from two-dimensional particle interactions, Phys. Rev. E 84, 015203(R) (2011).
43. A.L. Bertozzi, J. von Brecht, H. Sun, T. Kolokolnikov, and D. Uminsky, Ring Patterns and Their Bifurcations in a Nonlocal Model of Biological Swarms, Comm. Math. Sci. (to appear).
44. H. Cabral and D. Schmidt, Stability of Relative Equilibria in the Problem of Vortices, SIAM J. Math. Anal. 31, 231 (2000).
45. A. Barry, G. Hall, and C. Wayne, Relative Equilibria of the -Vortex Problem, J. Nonlin. Sci. 22, 63 (2012).
46. R.C. Fetecau, Y. Huang, and T. Kolokolnikov, Swarm dynamics and equilibria for a nonlocal aggregation model, Nonlinearity 24(10), 2681 (2011).
47. A.J. Bernoff and C.M. Topaz, A Primer of Swarm Equilibria, SIAM J. Appl. Dyn. Syst. 10, 212 (2011).
48. L.O. Baksmaty, S.J. Woo, S. Choi, and N.P. Bigelow, Tkachenko Waves in Rapidly Rotating Bose-Einstein Condensates, Phys. Rev. Lett. 92, 160405 (2004).
49. G. Baym, Tkachenko Modes of Vortex Lattices in Rapidly Rotating Bose-Einstein Condensates, Phys. Rev. Lett. 91, 110402 (2003).
50. Y. Chen and T. Kolokolnikov, A minimal model of predator-swarm interactions, J. Royal Soc. Interface 11:20131208 (2014)
51. H. Pu, L.O. Baksmaty, S. Yi, and N.P. Bigelow, Structural Phase Transitions of Vortex Matter in an Optical Lattice, Phys. Rev. Lett. 94, 190401 (2005).
52. M. Cipriani, M. Nitta, Crossover between Integer and Fractional Vortex Lattices in Coherently Coupled Two-Component Bose-Einstein Condensates, Phys. Rev. Lett., 111, 170401 (2013)
53. T. Simula, Collective dynamics of vortices in trapped Bose-Einstein condensates, Phys. Rev. A, 87, 023630 (2013)
54. T.P. Simula, and K. Machida, Kelvin-Tkachenko waves of few-vortex arrays in trapped Bose-Einstein condensates, Phys. Rev. A, 82, 063627 (2010)
55. N.S. Ginsberg, J. Brand, and L.V. Hau Observation of Hybrid Soliton Vortex-Ring Structures in Bose-Einstein Condensates, Phys. Rev. Lett. 94, 040403 (2005).
56. S. Komineas and J. Brand, Collisions of Solitons and Vortex Rings in Cylindrical Bose-Einstein Condensates, Phys. Rev. Lett. 95, 110401 (2005).
57. S. Komineas, Vortex rings and solitary waves in trapped Bose-Einstein condensates, Eur. Phys. J. Spec. Top. 147, 133 (2007).
58. M. Konstantinov, Chaotic phenomena in the interaction of vortex rings, Phys. Fluids 6, 1752 (1994).
101717