Single and Multiple Vortex Rings in Three-Dimensional Bose-Einstein Condensates: Existence, Stability and Dynamics

TL;DR

This paper investigates the existence, stability, and dynamics of single and multiple vortex rings in three-dimensional Bose-Einstein condensates, extending previous bifurcation studies to the Thomas-Fermi regime and analyzing their stability and interactions.

Contribution

It extends the analysis of vortex ring states to the large chemical potential regime and explores their stability, dynamics, and interactions using a particle picture approach.

Findings

Vortex rings exist and are stable in the Thomas-Fermi limit.

Multi-ring states exhibit specific equilibrium radii and inter-ring distances.

Instability scenarios include breakup, reconnections, and vortex line formation.

Abstract

In the present work, we explore the existence, stability and dynamics of single and multiple vortex ring states that can arise in Bose-Einstein condensates. Earlier works have illustrated the bifurcation of such states, in the vicinity of the linear limit, for isotropic or anisotropic three-dimensional harmonic traps. Here, we extend these states to the regime of large chemical potentials, the so-called Thomas-Fermi limit, and explore their properties such as equilibrium radii and inter-ring distance, for multi-ring states, as well as their vibrational spectra and possible instabilities. In this limit, both the existence and stability characteristics can be partially traced to a particle picture that considers the rings as individual particles oscillating within the trap and interacting pairwise with one another. Finally, we examine some representative instability scenarios of the multi-ring dynamics including breakup and reconnections, as well as the transient formation of vortex lines.

Equations

$$i\hslash {\psi }_t=-\frac{{\hslash }^2}{2m}{\nabla }^2\psi +V\psi +g|\psi {|}^2\psi -\mu \psi ,$$

$$V(\vec{r})=\frac{1}{2}m({\omega }_r^2{r}^2+{\omega }_z^2{z}^2),$$

$$\psi (\vec{r},t)={\psi }_0(\vec{r})+ϵ\left(a(\vec{r})e^{\lambda t}+b^{\ast }(\vec{r})e^{{\lambda }^{\ast }t}\right).$$

$$E$$

$$E_{\text{VR-T}}$$

$${\tilde{E}}_{\text{VR-VR}}=4\pi \sum _{\begin{array}{c}i,j=1\\ i\ne j\end{array}}^N{m}_i{m}_j\sqrt{r_i{r}_j}C(k_{ij}),$$

$$C(k_{ij})$$

$$k_{ij}^2$$

$${\dot{p}}_i=-\frac{\partial E}{\partial q_i}\text{and}{\dot{q}}_i=\frac{\partial E}{\partial p_i}.$$

$$(p_i,q_i)=(2\pi m_i{r}_i^2,z_i).$$

$$r_c=\sqrt{\frac{2\mu }{3{\omega }_r^2}},$$

$${\dot{r}}_i$$

$${\dot{z}}_i$$

$${\dot{r}}_i^{(\text{VR-T})}$$

$${\dot{z}}_i^{(\text{VR-T})}$$

$$L$$

$$\alpha$$

$${\dot{r}}_i^{(\text{VR-VR})}$$

$${\dot{z}}_i^{(\text{VR-VR})}$$

$$W_{ij}$$

Full text

I Introduction

Bose-Einstein condensates (BECs) of ultracold atomic gases have for around two decades now book1 ; book2 ; dumitr1 ; dumitr2 captured the interest not only of the atomic, molecular, and optical physics communities, but also considerably so that of the nonlinear waves emergent ; book_new . This, to a significant degree, is due to the effective nonlinearity introduced at the lowest-order mean-field theory book1 ; book2 ; book_new , which leads to a nonlinear Schrödinger type equation abl1 ; abl2 ; abl3 ; abl4 , referred to as the Gross-Pitaevskii equation (GPE). Depending on the sign of the $s$-wave scattering length, this results in a self-defocusing or self-focusing nonlinearity in the equation which, in turn, admits an array of relevant nonlinear wave structures. The latter comprise, but are not limited to, bright expb1 ; expb2 ; expb3 , gap gap and dark djf matter-wave solitons. In higher dimensions the main structures are vortices fetter1 ; fetter2 , in two dimensions, as well as vortex lines and rings, in three dimensions komineas_rev .

Our focus in the present work will be based on vortex rings (VRs) arising in three-dimensional (3D) repulsive BECs. Such structures have been studied rather extensively in theoretical, computational, and experimental works, with a number of reviews having emerged from this activity emergent ; komineas_rev ; book_new . Despite the wide body of knowledge existing from the physics of fluids and superfluids saffman ; Pismen1999 and the much earlier experimental observation of these structures in helium donnelly ; Rayfield64 ; Gamota73 , the amount of experimental attention that VRs have received in BECs is less than might be expected given their importance. Nevertheless, they have been experimentally observed in a range of different studies including, their realization through the decay of planar dark solitons in two-component BECs Anderson01 , creation via density engineering Shomroni09 , spontaneous emergence via the collision of symmetric defects Ginsberg05 and their detection through unusual collision outcomes of dark solitonic structures sengstock . Recently, intriguing six-vortex-ring and six-vortex-line cages have been predicted as transient states in the decay of spherical-shell dark solitons dsnew .

In recent years, a program of exploring coherent structures in a two-pronged way has naturally emerged and has been summarized, e.g., in Ref. book_new . On the one hand, it is possible to study nonlinear waves in the vicinity of the linear (noninteracting) limit of the GPE, namely that of the quantum harmonic oscillator. This limit, while physically less relevant given its association with small atom numbers, is very insightful towards the waveforms that are possible via different combinations of the linear eigenfunctions; incidentally, this regime is also more prone to important quantum-fluctuation phenomena prouk . In the context of VRs and related states, this approach has been utilized, e.g., in Refs. danaila ; ourpra1 to construct single and multiple VR structures. On the other hand, a complementary approach that is certainly relevant experimentally is the exploration of the highly nonlinear, large atom number regime. This is known as the Thomas-Fermi (TF) limit, in which the structures become narrower as the healing length, which constitutes their characteristic length scale, decreases. Here, the coherent waves can be thought of as individual “particles” that have kinematics and inter-particle dynamics that can be approximated by suitable particle models. Although attempts have been made to explore a single VR as such a particle, to describe its vibrational modes jackson99 ; fetterpra ; horng ; ourpra2 , as well as to explore multiple VRs in a homogeneous (untrapped) medium konstantinov ; Shashikanth03 ; caplan , to the best of our knowledge no such attempt has been made for the case of multiple trapped VRs. It is the purpose of this work to contribute to this direction by developing a systematic approach to study this problem. We thus show that multiple trapped VRs are especially important: this is due to the fact that the interplay of the intrinsic dynamics of each VR, with the trapping and inter-VR interactions, produce the possibility of stationary multi-VR states that we will argue —based on their stability properties— are experimentally accessible. It is also shown that multiple VRs also produce intriguing vibrational dynamics, within both the stable and unstable regimes.

The layout of the paper is as follows. In Sec. II.1 we outline the numerical formulation of the three-dimensional GPE and its Bogoliubov-de Gennes (BdG) spectral stability analysis. In the theoretical formulation of Sec. II.2, we present the “particle picture” (PP) that we will use to provide insights into the single and multiple VR steady states. Then, in Sec. III, we discuss our numerical results; firstly, for the cases of the single and double VRs, and finally for that of a triple VR state that we can also systematically construct and probe. We make direct comparisons between the treatments to illustrate the qualitative ability of the theory to capture multi-VR states. Our existence and stability computations, together with the PP analysis, are complemented with a selection of dynamical manifestations of both the vibrational modes —such as those that describe the relative motion of multiple VRs— and the intriguing instabilities in certain regimes. Finally, we summarize our findings, discuss some open problems and present a number of possibilities for future studies in Sec. IV.

II Theoretical and Computational Model Setup

II.1 The Gross-Pitaevskii equation

In the framework of the lowest-order mean-field theory, for sufficiently low temperatures, the dynamics of a 3D repulsive BEC is well-described by the GPE book1 ; book2 ; emergent ; book_new ,

[TABLE]

where $\psi(x,y,z,t)$ is the macroscopic wavefunction, $\mu$ is the chemical potential and $g=4\pi\hbar^{2}a_{s}/m$, with $a_{s}$ being the $s$-wave scattering length and $m$ the atomic mass. We use a harmonic trap of the form

[TABLE]

where $r^{2}=x^{2}+y^{2}$, while $\omega_{r}$ and $\omega_{z}$ are the trapping frequencies along the radial and axial directions, respectively. Note that the potential has rotational symmetry with respect to the $z$-axis. In our simulations we study the two-VR state with $\omega_{z}=\omega_{r}$, while we explore the three-VR state with $\omega_{z}=\frac{2}{3}\omega_{r}$. In these settings, the two states have chemical potentials at the linear (non-interacting) limit equal to $\mu_{c}=\frac{7}{2}\hbar\omega_{r}$ and $\frac{10}{3}\hbar\omega_{r}$, respectively. This occurs when the ring-dark-soliton state becomes energetically degenerate with either the two or three-planar-dark-soliton state, enabling combinations of two states with a relative phase difference of $\pi/2$ that give rise to the two-VR or three-VR state, respectively danaila ; ourpra1 . These states can then be numerically followed all the way from the above mentioned linear limit to the large chemical potential, highly nonlinear limit. We have also performed the corresponding BdG spectral analysis of the states as a function of the chemical potential using the linearization ansatz,

[TABLE]

Here, $\psi_{0}$ denotes the single- or multi-VR stationary state whose stability is sought, $\epsilon$ is a formal perturbation parameter and $\lambda$ denotes the eigenvalue corresponding to the eigenvector $(a,b)^{T}$, with $(\cdot)^{T}$ denoting the transpose and $(\cdot)^{*}$ complex conjugation.

We temporally evolve the GPE, Eq. (1), using two independent methods, providing a strong test of our numerics. In one method, we utilize a real-space product scheme with a finite-element discrete-variable representation. This is based on a split-operator approach, and a Gauss-Legendre quadrature is implemented within each element Schneider2006 ; dvr-exp . In the other method, a split-step operator is performed on a Fast-Fourier-Transform (FFT) grid. Our BdG calculations are made feasible by implementing a Fourier-Hankel scheme that utilizes azimuthal symmetry, much as was done in Refs. Ronen2006 ; Blakie2012 . In addition to allowing one to treat 3D functions as 2D, numerically, this allows the diagonalization of each angular momentum subspace individually.

In what follows, we are interested in both the properties of the GPE solutions themselves, such as the equilibrium radial or axial positions, and those of their excitations. We will also attempt to connect the GPE solutions with our particle picture results.

II.2 The particle picture of vortex rings in a trap

We now focus on the limit of large chemical potentials $\mu$, where we can provide a theoretical analysis of the in-trap dynamics and interactions of multiple VRs. For simplicity, in all that follows we have used dimensionless units (see, e.g., Refs. emergent ; book_new ) where time is measured in units of inverse trap frequency ($1/\omega_{r}$), length scales in units of harmonic oscillator length ($a_{r}=\sqrt{\hbar/m\omega_{r}}$) and energy units of $\hbar\omega_{r}$. In this TF limit, there exists a well known approximation to the ground state of the GPE given by $\psi_{\text{TF}}=\sqrt{\max[\mu-V(\vec{r}),0]}$. Here, our aim is to explore vortical excitations (in particular, VRs) on top of this ground state. As indicated in the introduction, studies have independently considered each of the following:

• (i)

the motion of a single VR in a homogeneous medium roberts ,

• (ii)

the effect of a trap on a single VR JacksonMcCannAdams ; jackson99 ; fetterpra ; horng , and the

• (iii)

interaction of multiple VRs in the absence of a trap konstantinov ; Shashikanth03 ; caplan .

Another aim of this work, in addition to exploring the existence, stability, and dynamics of these states, is to explore the effective PP model arising from combining these different ingredients together, and its usefulness in capturing the essential static properties corresponding to equilibrium configurations of multiple VRs in a trap as well as the periodic oscillations ensuing from initial conditions away from these equilibria.

For a set of co-axial VRs along the $z$-axis, a naïve approach to combine the above-mentioned VR-VR and VR-trap contributions would consist of simply adding the corresponding reduced dynamics at the level of the effective ordinary differential equations (ODEs) on the VR radii $r_{i}$ and positions $z_{i}$. However, perhaps somewhat surprisingly, this approach turns out to produce non-Hamiltonian ODEs because the two main contributions, namely the VR-VR interaction konstantinov and VR-trap interaction JacksonMcCannAdams , originate from energy terms with different canonical variables (see below). Therefore, this approach —although capable of reasonably predicting the positions for stationary multi-VR configuration (results not shown here)— fails to describe the actual dynamics of trapped multi-VR configurations. In fact, the ensuing VR dynamics for this naïve approach give rise to damped/anti-damped orbits (depending on the VR charges) that tend to spuriously spiral in/out of the condensate.

In order to derive a self-consistent, Hamiltonian set of equations for trapped multi-VR configurations it is necessary to start from the Hamiltonian formulation of the different interaction energies involved and then obtain the equations of motion through a common set of canonical variables. Let us then consider the following energies: (i) the VR-trap energy, denoted by $E_{\text{VR-T}}$, described in Ref. JacksonMcCannAdams and (ii) the VR-VR energy, denoted by $E_{\text{VR-VR}}$, described in Ref. konstantinov . Importantly, both $E_{\text{VR-T}}$ and $E_{\text{VR-VR}}$ contain the VR self-induced velocity that is responsible for a single VR to always have an intrinsic velocity. Therefore, we construct the total energy of the system with the following combination that only includes the self-induced interaction once:

[TABLE]

where $E_{\mbox{\scriptsize VR-VR}}^{\text{self}}$ corresponds to the contribution to the energy originating from the self-induced velocity in the untrapped VR-VR description konstantinov . Specifically, using Refs. JacksonMcCannAdams and konstantinov , to describe these contributions yields the following energies. The VR-T contribution for a single VR at position $(r_{i},z_{i})$ inside an isotropic ($\omega=\omega_{r}=\omega_{z}$) trapping potential with TF radius $R_{\perp}=\sqrt{2\mu}/\omega$ yields JacksonMcCannAdams

[TABLE]

with $r_{0}^{2}=r_{i}^{2}+z_{i}^{2}$ and where $\xi=1/\sqrt{a\mu}$ is obtained from the $r\approx 0$ asymptotics of the vortex core density $\rho(r)\approx a\mu^{2}r^{2}$ and $a=0.82226$ was computed numerically by fitting this asymptotic expression. On the other hand, the VR-VR contributions, without self-induced velocity terms, for a set of $N$ VRs of charge $m_{i}$ and position $(r_{i},z_{i})$ yield konstantinov

[TABLE]

where

[TABLE]

and ${\cal K}$ and ${\cal E}$ are, respectively, the complete elliptic integrals of the first and second kind and $k_{ij}$ is their respective elliptic modulus.

Having the total energy (4), i.e., the Hamiltonian for the system of interacting VRs, we can obtain the equations of motion as a set of coupled ODEs for the vortex positions using the corresponding Hamilton’s equations:

[TABLE]

We follow the choice of canonical variables $(p_{i},q_{i})$ from Ref. konstantinov :

[TABLE]

It is precisely due to this choice of canonical variables that it is not possible to apply the naïve approach of adding the resulting equations on motion from Refs. konstantinov and horng as the former uses the canonical variables (10) while the latter uses $(p_{i},q_{i})=(r_{i},z_{i})$. For completeness, the resulting ODEs are included in the Appendix.

III Results

We begin the discussion of our numerical findings by presenting existence results for the case of single and multiple VRs. The top row of panels in Fig. 1 depicts illustrative examples of the steady-state configurations for one and two VRs. However, if the VRs are not placed at their stationary position, they start to oscillate as depicted in the bottom row of panels in Fig. 1. These panels depict the trajectories for one VR (left panel) and two oppositely charged VRs (right panel) in an isotropic parabolic trap (2) with $\omega_{r}=\omega_{z}=1$ and chemical potential $\mu=50$. The thick and thin trajectories correspond, respectively, to our PP and the GPE dynamics. The figure suggests that the dynamics for a single VR is qualitatively and quantitatively very well described by the PP. The case of two VRs (right panel) indicates a good qualitative match between PP and original GPE dynamics, yet the position of the steady state configurations seems to be slightly shifted.

In order to have a broader sense of the validity of the PP, let us follow the steady state VR positions for one and two VRs as the chemical potential is varied. Figure 2 depicts the comparison between the numerical solutions of the GPE (solid lines) and the PP (dashed lines). For a single VR (top panel), the PP does an excellent job at predicting the stationary position of the VR and its functional dependence on the chemical potential $\mu$. Note that the $\mu\rightarrow\infty$ asymptotic prediction in the TF limit from Refs. kaper ; dsnew for the equilibrium radius $r_{c}$:

[TABLE]

is slightly higher than the one predicted by the PP.

The bottom panel of Fig. 2 depicts the comparison between the GPE and PP results for the two-VR steady state configuration.

In this case, one can see that the trends predicted by the PP can qualitatively follow that of the full GPE, although some quantitative disparity remains. This can be attributed to the following causes:

• (i)

the PP is an amalgamation of different contributions stemming from different approaches; it would be useful, although technically seemingly especially tedious to consider an approach incorporating all three effects concurrently.

• (ii)

The accumulation of errors for the three different contributions, since each one of them involves corresponding approximations.

• (iii)

the interaction between the VRs is modulated by density variations, a feature that is not captured in the effective PP in the form considered herein.

This last point is also an issue that affects similar particle approaches for vortices in 2D, and has been discussed in Refs. busch ; middelkamp . Nevertheless, we conclude that the PP approach, while less quantitatively dependable, can be used to provide a qualitative handle on the trends of stationary multi-VR characteristics.

We now turn to the exploration of the spectrum of the multi-VR states. Although a brief discussion of this can be found in Ref. ourpra1 , in the vicinity of the linear limit, here we consider the relevant spectrum more systematically for a wide parameter range. The relevant spectrum for the two-VR solution is shown in Fig. 3, and is seen to bear numerous similarities to that of the single VR (discussed in Ref. ourpra2 ). In particular, the modes associated with dynamical instability are fairly similar to those of the single VR. The most significant one among them, associated with the widest parameter range of instability, is related to a quadrupolar mode having an $n=2$ azimuthal dependence $e^{in\theta}$, as we will also see below in the dynamical-evolution results. On the other hand, the mode immediately above this one, for large $\mu$ in the spectrum, is azimuthally symmetric ($n=0$) and is associated with a stable relative motion of the two vortex rings. The trajectory of this mode is illustrated in Fig. 4. This excitation was initiated by mixing the two-VR stationary state with the corresponding BdG excitation; the resulting wavefunction was then renormalized to preserve the total atom number and subsequently evolved in time according to the GPE. Since this excitation has $n=0$, these rings have azimuthal symmetry, and the cores remain circular and coaxial with the $z$-axis throughout the motion. In the ($r,z$) cross section of Fig. 4(c), the upper ring orbits clockwise while the lower ring orbits counterclockwise in such a way that their $z$ motion remains in-phase. The fitted period for this oscillation is $T=27.1$ (in units of $1/\omega_{r}$), which compares favorably with the BdG prediction $T=27.9$ (see black circle in Fig. 3). Note that the chemical potential is not well-defined in the absence of a stationary state, of which was (in dimensionless units) $\mu\sim 16$ before the addition of the excitation.

Contrary to the single VR studied in Ref. ourpra2 , the existence of an additional ring means that there is a second $n=0$ vortex excitation. This is shown in Fig. 5, and exhibits a similar motion to the first excitation but now the measured period is $T=11.9$, compared to the BdG period of $T=12.0$ (see black triangle symbol in Fig. 3). Furthermore, the relative phase of the ring motion differs such that their $r$ positions are now in-phase, while the $z$ positions now perform an out-of-phase oscillation. Finally, the branches that level off at large $\mu$ are the non-vortex excitations of the underlying ground state.

We also provide a similar spectral perspective in the case of the three-VR solution in Fig. 6. The top quartet and the middle isocontour panel illustrate a characteristic example of this state for (dimensionless) chemical potential $\mu=12$. In the cut along the central ($x,y$) plane, a ring-like structure is clearly discernible, while a vertical ($y,z$)-cut reveals three pairs of opposite-charge vortices, alternating along the $y$-axis, as well as along the $z$-axis. This implies that the three rings, given the alternating nature of the constituent triple soliton pattern, result in an equilibrium triplet of vortex rings of charge $+1$, $-1$ and $+1$, up to parity reversals of all three. The middle panel clearly showcases, through its density isocontours, the nature of the pattern. The bottom panel in the figure depicts the BdG spectrum for the three-VR solution as the chemical potential is varied. While some of the oscillatory instabilities associated with complex eigenvalues appear over narrow parameter intervals, for small values of $\mu<10$, there is an instability arising from the collision of two eigenmodes near $\mu=5$ that seems to persist for $\mu>10$ and indeed for the entire parametric interval of chemical potentials that we have explored. We will explore this dynamical mode in the direct numerical simulations that follow, as it seems the most pertinent one to the potential stability of multi-VR states in the TF regime of large $\mu$. We note that this instability might be cured by adjusting the trap aspect ratio. Again, we have checked that the modes that level off are bulk excitations of the underlying ground state.

We now explore some of the key dynamical features of multi-VRs through direct numerical simulations. In Fig. 7, we examine a rather intriguing example of the dynamical instability of the two-VR state for $\mu=10$. In this figure we show the 3D positions of the vortex core as red curves foster and their projection onto the $(x,y)$, $(x,z)$ and $(y,z)$ planes are shown as black curves. 2D density contours are also projected on these three planes and the contours correspond to 0.25, 0.5, and 0.75 of the maximum density at each time. We observe that the rings initially deform in a quadrupolar fashion, and once this deformation becomes sufficiently severe, they split and reconnect as (nearly) co-axial rings perpendicular to their original orientation. The rotated double rings remain robust for a considerable time interval, as can be seen in frames Fig. 7(c)–(e) which span $\Delta t\sim 25$. We have also observed that, depending on the initial conditions, two VRs can undergo this effective rotation several times, sometimes transiently breaking into four vortex lines, as shown in Fig. 7(f), before reforming as a rotated double VR. Eventually, though, the VRs break up into a vortex tangle.

Lastly, we examine a prototypical example of the instability associated with the three-VR state. We probe, in particular, the unstable mode that was observed in Fig. 6 to be the persistent cause of instability for large values of the chemical potential. This instability in the case of $\mu=9$ is illustrated in Fig. 8. There, it can be observed that the triple VR becomes subject to a symmetry-breaking deformation of the rings. As time evolves, the resulting undulations are amplified and lead to intense Kelvin mode excitations along each of the rings. Subsequently, two of these rings may combine, as can be seen in Fig. 8(d). These two rings then separate again, as shown in Fig. 8(e), and eventually all three rings again evolve into a vortex tangle as shown in Fig. 8(f).

IV Discussion and Future Challenges

In the present work, we have extended the analytical and numerical analysis concerning the dynamics of single and multiple VRs in harmonically confined Bose-Einstein condensates. We showed that in such systems it is controllably possible to form and establish the existence of states involving one-, two-, three- or more VRs, essentially at will over a wide range of chemical potentials —or, equivalently, atom numbers. Furthermore, we provided a theoretical formulation based on the energetics of the different processes involving the rings at the particle level. We demonstrated that for the 1-, 2- and 3-VR states considered herein this approach, encompassing the self-induced ring translation, the trap induced ring oscillation and the inter-ring interaction, could capture qualitatively and in some cases quantitatively the ring steady states. The stability for these multi-VR steady states was also probed at the level of a full Bogoliubov-de Gennes analysis which revealed both the instabilities of, e.g., 2- and 3-VR states, but also the collective relative motions of the rings. Finally, some selected examples of the dynamics of the VRs were displayed, demonstrating both their potential for coherent multi-VR motions, and also their rich instability scenarios containing examples of VR breakups and recombinations. Such events include ring mergers and subsequent splits or redistribution of the vorticity in the form of vortex lines.

There are numerous questions that remain open and constitute interesting directions for future work.

First, in the context of a single vortex ring, identifying an approach, perhaps involving an adiabatic invariant as in Refs. ba ; kono ; kamch ; ai ), that would give a more accurate description of its equilibrium and overall dynamics would be helpful both in the realm of single and in that of multiple VRs.

Secondly, regarding the case of more than one VR, it is important to mention that, although the PP gave a qualitative description of the size and location of the multiple VR steady-state configuration, it fails to accurately describe the precise location of the fixed points in a quantitative fashion. Nonetheless, it is interesting to note that, despite this shortcoming, the reduced PP equations of motion are able to capture (results not shown here): (a) all the different modes of vibration of the VRs (for two VRs the in-phase and our-of-phase oscillations and for three VRs the three normal modes of vibration) and (b) the trend corresponding to the actual frequency of vibration for these oscillatory vibrational modes as the chemical potential in varied.

In the same vein, it would be interesting, albeit technically demanding, to attempt to incorporate the density modulations in the inter-ring interactions, in a similar way to what was proposed for ordinary vortices in Ref. busch . A perhaps more straightforward, although computationally demanding possibility, would be to explore the dynamics of VRs in a way similar to that of Ref. baggaley . In this work, the authors use the full Biot-Savart law in order to explore the interactions between numerous initially co-axial vortex rings in a uniform medium with the aim of studying their leapfrogging behavior. The Biot-Savart law encompasses two out of the three features we are considering here, namely the ring self-action and the inter-ring interaction. It does not include the effect of the trap, yet this can be accounted for through the work of Ref. fetterpra . The combination of these approaches would then enable at an effective, rather than full GPE, level the examination of not solely co-axial vortex rings as is done here but rather allowing their deformation (beyond co-axiality) while considering their motion. Comparing such a sophisticated approach with the full Gross-Pitaevskii dynamics would be especially informative.

Naturally, also, extensions of the present considerations to a higher number of components, including the potential formation of more exotic structural configurations such as Skyrmions ruost1 ; ruost2 would be of interest as well for future work.

Progress along some of these directions, as relevant, will be reported in future publications.

Appendix A Effective ODEs for trapped VRs

Here we include the explicit effective ODEs describing the dynamics of $N$ VRs with radii $r_{i}$ and vertical positions $z_{i}$ inside an isotropic trap with strength $\omega$ and TF radius $R_{\perp}=\sqrt{2\mu}/\omega$ for chemical potential $\mu$. By using Hamilton’s equations (9), the dynamics for each VR is given by:

[TABLE]

where the contributions due to the trap are encapsulated in the VR-T terms and the vortex-vortex interactions in the VR-VR terms. The VR-T contributions yield

[TABLE]

where

[TABLE]

On the other hand, the VR-VR contributions yield

[TABLE]

where

[TABLE]

and $C(k_{ij})$ is defined in Eq. (7).

Acknowledgements.

W.W. acknowledges support from NSF-DMR-1151387 and from the Office of the Director of National Intelligence (ODNI), Intelligence Advance Research Projects Activity (IARPA), via MIT Lincoln Laboratory Air Force Contract No. FA8721-05-C-0002. R.N.B. is supported by the QUIC grant of the Horizon2020 FET program and by Provincia Autonoma di Trento. R.N.B., C.T., L.A.C. and P.G.K. acknowledge support by Los Alamos National Laboratory, which is operated by LANS, LLC for the NNSA of the U.S. DOE and, specifically, Contract No. DEAC52-06NA25396. R.C.G. gratefully acknowledges the support of NSF-DMS-1309035 and PHY-1603058. P.G.K. gratefully acknowledges the support of NSF-DMS-1312856 and PHY-1602994, from the ERC under FP7, Marie Curie Actions, People, International Research Staff Exchange Scheme (IRSES-605096), and the Stavros Niarchos Foundation via the Greek Diaspora Fellowship Program. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purpose notwithstanding any copyright annotation thereon. We thank the Texas A&M University for access to their Ada cluster.