The dynamics of vortex ring pairs in the homogeneous nonlinear Schrödinger equation is studied. The generation of numerically exact solutions of traveling vortex rings is described and their translational velocity compared to revised analytic approximations. The scattering behavior of co-axial vortex rings with opposite charge undergoing collision is numerically investigated for different scattering angles yielding a surprisingly simple result for its dependence as a function of the initial vortex ring parameters. We also study the leapfrogging behavior of co-axial rings with equal charge and compare it with the dynamics stemming from a modified version of the reduced equations of motion from a classical fluid model derived using the Biot-Savart law.

## I. INTRODUCTION

One of the most widespread and interesting models for studying the emergence, dynamics and interactions of coherent structures is the nonlinear Schrödinger equation (NLSE).^{1} The NLSE is a paradigm for the evolution of coherent structures since it is a universal model describing the evolution of complex field envelopes in nonlinear dispersive media.^{2} This universality stems from the fact that the NLSE is the prototypical, lowest-order, nonlinear partial differential equation (PDE) that describes the dynamics of modulated envelope waves in a nonlinear medium.^{3} As such, it appears in a wide variety of physical contexts, ranging from optics^{4–7} to fluid dynamics and plasma physics^{8} to matter waves,^{9,10} while it has also attracted much interest mathematically.^{1,11–15}

The general form of the NLSE (with the lowest order nonlinearity) can be written, in non-dimensional units, as

where *a* is an adimensional parameter (usually *a* = 1 or *a* = 1/2 depending on the choice of scaling) and *s* represents the strength of the nonlinear interaction and determines whether the NLSE is attractive (*s* > 0) or repulsive (*s* < 0). The NLSE is a model for the evolution of the mean-field wavefunction in Bose-Einstein condensates (BECs) at low temperatures,^{16,17} namely, a superfluid. In the BEC setting, the NLSE also contains an external trapping potential term that is often spatially parabolic or periodic in nature and stems from the magnetic or optical confinement of the atoms.

The NLSE allows for the prediction and description of a wide range of (nonlinear) excitations depending on the sign of the nonlinearity and the dimensionality of the system, the latter of which can be tuned by the ratio of trapping strengths along the different spatial directions.^{9} In this paper, we are interested in the fully three-dimensional (3D) regime with repulsive nonlinearity (*s* < 0) whose basic nonlinear excitations correspond to vortex lines and vortex rings. A vortex line is the 3D extension of a two-dimensional (2D) vortex by (infinitely and homogeneously) extending the solution into the axis perpendicular to the vortex plane. Vortex lines might be rendered finite in length if the background where they live is bounded by the externally confining potential. In that case, vortex lines become vorticity “tubes” that are straight across the background or bent in U and S shapes, depending on the aspect ratio of the background.^{18,19} If a vortex line is bent enough to close on itself, or if two vortex lines are close enough to each other, they can produce a vortex ring.^{20} Vortex rings are 3D structures whose core is a closed loop with vorticity around it^{21} (i.e., a vortex line that loops back into itself). Vortex rings have been observed experimentally in superfluid helium^{22,23} as well as in the context of BECs in the decay of dark solitons in two-component BECs,^{24} direct density engineering^{25} (see also the complementary theoretical proposals of Refs. 26 and 27), and in the evolution of colliding symmetric defects.^{28} They have also been argued to be responsible for unusual experimental collisional outcomes of structures that may appear as dark solitons in cigar-shaped traps.^{29} Numerical studies of experimentally feasible vortex ring generation in BECs have also been explored in the context of flow past an obstacle,^{30,31} Bloch oscillations in an optical trap,^{32} collapse of bubbles,^{33} instability of 2D rarefaction pulses,^{34} flow past a positive ion^{35,36} or an electron bubble,^{37} crow instability of two vortex pairs,^{38} and collisions of multiple BECs.^{39,40} Vortex rings have an intrinsic velocity,^{41} which can be overcome by counteracting the velocity with a trapping potential,^{42} adding Kelvin mode perturbations,^{43,44} or by placing the vortex ring in a co-traveling frame (see Sec. II).

In the present paper, motivated by the above abundant interest in this theme both from a theoretical perspective and from that of ongoing experimental efforts in BECs, we are interested in the dynamical evolution of vortex rings and in particular in their interactions. In an effort to exclusively capture the vortex ring dynamics and not the influence of the external potential (the latter is a natural subject for future work), we assume the absence of any trapping potential and thus consider the background supporting the vortex rings as homogeneous. Namely, we consider vortex rings embedded in a constant density background. On this homogeneous background we revisit the intrinsic velocity of a single vortex ring that serves as the self-interacting term that will be supplemented by the interaction terms between different vortex rings. The main goal of our work is to describe vortex ring interactions and scattering collisions in a reduced (i.e., “effective”) manner. The motivation for this approach stems from the fact that full numerical simulations in the 3D NLSE can be very time consuming especially when examining scenarios involving multiple vortex rings in large domains. The first dynamical scenario that we study corresponds to the scattering of colliding vortex rings of opposite charge. Using a large set of full 3D numerical simulations (based on efficient GPU codes, see below), we are able to distill a very simple, phenomenological, rule describing the scattering angle of vortex rings as a function of the vortex ring radii and the initial collisional offset. Then, based on the interaction of vortex filaments from classical fluids, we revisit a reduced, ordinary differential equation, model for the interaction of co-axial vortex rings based on the Biot-Savart law, adapting it to the superfluid (BEC) case. The ensuing reduced dynamics for the leapfrogging of vortex ring pairs is then favorably compared to full 3D simulations.

Our paper is structured as follows. In Sec. II, we describe a procedure to generate numerically exact traveling vortex ring solutions in a homogeneous background, and test the resulting ring's translational velocity against suitably revised known analytic approximations. Section III is devoted to describing scattering scenarios from the collisions of vortex rings, where we derive a very simple phenomenological relationship for the scattering angle as a function of the initial offset distance and radii of the colliding rings. In Sec. IV, we use effective equations of motion for the interaction of co-axial vortex rings to describe the leapfrogging behavior of two vortex rings and compare the results with direct simulations of the NLSE. Finally, Sec. V summarizes our results and prompts a few avenues for future exploration.

Before we embark on the relevant analysis, we provide some details on our computational methods which pertain to all the numerical simulations given below. For all 3D NLSE simulations, we use high-order explicit finite-difference schemes on a uniform grid with cell-spacing *h* = Δ*x* = Δ*y* = Δ*z* and constant time-step *k* = Δ*t*. Time-stepping is accomplished with the standard 4th-order Runge-Kutta method (RK4), while spatial differencing is performed with the 4th-order 2-step high-order compact scheme described in Ref. 45. The time-step *k* is set based on the stability bounds of the overall scheme as discussed in Ref. 46. Since the stability forces the time-step to be proportional to *h*^{2}, the overall error of the scheme is *O*(*h*^{4}). Due to the constant density background of the problem, we utilize a recently developed modulus-squared Dirichlet boundary condition.^{47} The simulations are computed using the NLSEmagic code package^{48,49} which contains algorithms written in C and CUDA, and are primarily run on NVIDIA GeForce GTX 580 and GeForce GT 650M GPU cards.

## II. GENERATION OF A NUMERICALLY EXACT NON-STATIONARY VORTEX RING

There exist various techniques to generate approximations to vortex ring configurations.^{50} However, in order to achieve accurate results for simulating the interactions of multiple vortex rings in a homogeneous background, it is necessary to first be able to generate numerically “exact” solutions of the individual rings. Since vortex rings have an intrinsic transverse velocity due to their topological structure, standard steady-state methods (such as imaginary time integration^{42} and nonlinear equation solvers) for obtaining solutions of coherent structures cannot be directly applied, unless the rings are in a configuration where they are at steady-state (such as in a magnetic trap^{42}). In order to be able to generate rings in a homogeneous background, we apply a co-moving background velocity to render the vortex ring stationary.

In order to apply the back-flow, we need a very accurate value of the vortex ring's transverse velocity. The transverse velocity of a single vortex ring in the NLSE has been studied analytically^{41,51,52} and the results were shown to be consistent with numerical simulations of the NLSE.^{53} For the NLSE in the form of Eq. (1), with *s* < 0 and *V* = 0, an asymptotically approximate velocity for the ring is given by^{41}

where *d* is the ring's radius, *m* is its the charge, and *r*_{c} is the vortex ring's core radius defined as

where ξ is the healing length given by

Ω plays the role of the frequency and is tantamount to the system's chemical potential. The value *L*_{0}(*m*) is referred to as the vortex core parameter and is defined as the convergent part of the energy-per-unit-length of a vortex line in the NLSE.^{41} While *L*_{0}(*m*) depends on the vortex ring charge *m*, it is independent of the NLSE parameters. In order to determine the velocity of the vortex ring as precisely as possible, *L*_{0}(*m*) must be computed accurately. Since there are various values of *L*_{0}(*m*) reported in the literature (cf. values given in Refs. 41,54,55), we briefly show in the Appendix how *L*_{0}(*m*) is numerically computed and give its values for *m* ∈ [1, 10]. Since we are only focusing on vortex rings of charge |*m*| = 1 (indeed, stable higher-charge vortex rings are not known to exist^{56}), we use the value *L*_{0}(1) ≈ 0.380868.

The initial condition for a single vortex ring solution in cylindrical coordinates (see Fig. 1) with radius *d* centered at (*r*, *z*) = (0, 0) can be described by

where

and the function *g*(*r*, *z*) is an axisymmetric 2D profile which may contain additional phase information of the solution. This 2D profile can be approximated by the modulus of the radial solution of a 2D NLSE vortex in the (*r*, ϕ^{+}) plane and is described by

where

*f*(ρ) is the numerical solution to the radial steady-state equation

which is derived from inserting the form of a vortex solution Ψ = *f*(ρ) exp [*i*(*m*ϕ^{+} + Ω*t*)] into the 2D NLSE in polar coordinates.

In order to obtain a numerically exact vortex ring initial condition, we would like to solve a steady-state equation for Ψ_{0} using the above approximation as an initial seed. This can be accomplished by noting that a steady-state solution to the NLSE of Eq. (1) in cylindrical coordinates in a co-moving frame in the *z*-direction with velocity *c* is given by^{57}

where *Z* = *z* − *ct* and *U*(*r*, *Z*, θ) solves the time-independent NLSE

Therefore, by imposing a counter-flow with a velocity equal and opposite to that of the intrinsic velocity of the vortex ring by the following Galilean boost:

where *c* is the analytical approximation of the vortex ring velocity given by Eq. (2), we can then solve Eq. (5) for *U*. We use the nonlinear Newton-Krylov solver nsoli^{58} to find the solution, using a central-difference discretization of the spatial derivatives, along with the modulus-squared Dirichlet boundary conditions mentioned in Sec. I. The initial condition of the vortex ring solution is then found by removing the added counter-flow velocity

The numerically exact vortex ring solutions can be used to compare the analytical approximations of the vortex ring's transverse velocity Eq. (2) to direct simulations of the NLSE. In order to track the vortex rings during the simulations, a 2D (*y*, *z*) cut of the computational grid is extracted, and a center-of-mass calculation of the modulus of the vorticity in the region of the ring's intersection points in the cut are tracked. In Fig. 2, we show the results of integrating vortex rings until *t* = 50 for radii *d* ∈ [2, 10] and comparing their tracked velocity to Eq. (2). It is seen that the direct velocity results match the analytical approximation very accurately; the disagreement is never more than 1.5% in our studies. This illustrates that values obtained from the asymptotic velocity equation, Eq. (2), match direct integration even for rings with relatively small radii. We note that although the velocity *c* is used to generate the vortex rings initially, the fact that the rings maintain the velocity over a long integration time with no noticeable noise or counter-flow in the solution confirms that the simulations demonstrate the validity of Eq. (2).

We note here that the transverse velocity discussed in this section is that of an *unperturbed* vortex ring. It is possible to perturb a vortex ring (or line) to produce oscillations along the ring called Kelvin modes (or Kelvons).^{59–62} A typical event that can produce Kelvin modes occurs when two vortex rings scatter off of each other (see Sec. III). These Kelvin modes not only have their own dynamics, interactions, and decay properties,^{63,64} but they can also self-interact within a vortex ring, resulting in a reduction or even reversal of its transverse velocity.^{43,44}

## III. SCATTERING VORTEX RINGS

We now present a quantitative analysis of the scattering of two unit-charge (|*m*| = 1) vortex rings in the NLSE. An initial qualitative study of the scattering of such rings was performed in Ref. 53, where the authors focused on head-on collisions of the rings at different angled orientations. In contrast, we focus here exclusively on offset co-planar collisions, not studied in Ref. 53.

As shown in Ref. 53, two colliding vortex rings with zero offset (axisymmetric) expand and annihilate each other's topological charge resulting in an axisymmetric decaying rarefaction pulse. An example of such a case is shown in Fig. 3 where two vortex rings of radius *d* = 6 (one with charge *m* = 1 and one with *m* = −1) are positioned a distance of 6*d* apart from each other in the *z*-direction and allowed to collide. We study here the scattering that results when the colliding vortex rings are offset from each other in the planar direction. Our numerical experiments show that, typically, the two rings merge and then separate into two new vortex rings moving away from each other at an angled trajectory, each exhibiting large quadrupole Kelvin oscillations.

To quantify the relationship between offset distance and scattering angle, the vortex rings are positioned at a distance of 6*d* apart in the *z*-direction, and in the *y*-direction they are placed with a planar offset distance (impact parameter) defined as *q* (see Fig. 4(a)). The scattering angle θ_{s} is defined as the angle away from vertical so that if the vortex rings travel in a straight trajectory, the angle is 0, while if they scatter at right-angles to each other, the angle is π/2 (see Fig. 4(a)). Multiple simulations of colliding rings were performed with planar offsets ranging from *q*/*d* = 0 to *q*/*d* = 3 in increments of 0.05 for vortex rings of radius *d* = 6, 8, 10. We found three different qualitative scenarios as depicted and described in Fig. 5. The centers of the vortex rings in the (*y*, *z*) plane are tracked and the scattering angle recorded. The results are shown in Fig. 6. It is clear from the results that there are three distinct sections in the relationship between scattering angle and offset, which correspond to qualitatively distinct topological outcomes as follows:

For offsets with

*q*/*d*∈ [0, 1/2) the two rings annihilate each other's topological charge, sending out a rarefaction pulse, see Fig. 5(a), similar to the cases of head-on collision shown in Fig. 3. However in this case, as the rings approach each other, they rotate causing the collisions to occur at the scattering angle. Another distinction of these collisions is that the resulting rarefaction pulse is not axisymmetric as in the head-on collision case, but more concentrated in the scattering direction.When the offset has

*q*/*d*∈ (1/2, 2), the two rings undergo merging immediately followed by splitting into internally excited rings with quadrupole oscillations (see Fig. 5(b)). Interestingly, when the offset is set to*q*/*d*= 3/2 or larger, the scattering angle becomes negative, but the qualitative behavior remains the same.For larger offset such that

*q*/*d*> 2, the two rings no longer merge and separate into new rings, but are merely perturbed by each other in a fly-by scenario causing their trajectories to travel at a small negative angle from incidence (i.e., get weakly deflected “inward”), see Fig. 5(c). Since, in this case, the vortex rings primarily interact through their decaying tails (i.e., the decay of the density dip away from the vortex core), it is expected that the scattering in this region also decays to zero as the offset increases. Thus, as observed in the figure, when the offset increases the interaction decreases until approximately*q*/*d*= 3 where the rings essentially ignore each other's presence and travel unperturbed.

Noting that the scattering angle as a function of offset in Fig. 6 appears linear in the first two qualitatively distinct intervals, simple phenomenological approximations of the relationships can be formulated. For the annihilation section [*q*/*d* ∈ (0, 1/2)], the angle of the rarefaction pulse's scattering angle can be approximated by

while for the scattering rings [*q*/*d* ∈ (1/2, 2)], the angle of the resulting rings can be approximated as

These approximations are shown along with the simulation data in Fig. 6. It is interesting to note that, despite the complex dynamics involved in the scattering of vortex rings, we are able to formulate a very simple, phenomenological scattering rule for co-planar collisions.

It is interesting to contrast the 3D scenario of co-planar vortex ring scattering to its 2D counterpart. If one focuses on the dynamics restricted on the plane defined by the vortex ring axes, each vortex ring induces a vortex pair in 2D. The case of scattering between 2D (point-) vortex pairs has been studied in some detail.^{65–67} The most notable differences between the ensuing dynamics for 2D vortex pairs and that of 3D vortex rings scattering is that the dynamics in 2D is apparently richer due to the fact that truly 2D vortices are more “free” to move (respecting the vortex-vortex interactions) while the vortex pairs produced by the cut of the vortex ring in the plane are more tightly linked to each other through the ring. This extra freedom for 2D vortices allows for a rich variety of dynamical evolution scenarios. For instance, depending on the initial configuration (initial positions, orientations, and vortex pair internal distances) the collision dynamics of two vortex pairs can display: (a) transient bound vortex pairs where the two vortex pairs interact by orbiting around each other and exchange partners for a few periods until scattering away from each other as two distinct vortex pairs, (b) transient interactions between the four vortices that are chaotic and lead to the eventual expulsion of two distinct vortex pairs resulting in chaotic scattering, or (c) transient bound states where three of the four vortices lock into a translational and precessional motion, while the fourth one moves separately, before the two pairs again scatter away, separating indefinitely.

## IV. LEAPFROGGING VORTEX RINGS

Let us now consider the interaction between two same sign unit-charge vortex rings (*m*_{1} = *m*_{2} = 1) aligned along the *z* axis. In classical fluids, the analogous setup exhibits leapfrogging behavior.^{68,69} Leapfrogging has also been recently reported in models of ferromagnetic fluids.^{70} In Fig. 7, we show an example of two equal-charge vortex rings of radius *d* = 10 and initially separated by a distance of *Q* = *Z*_{2}(0) − *Z*_{1}(0) = 5 in the *z*-direction undergoing leapfrogging dynamics. In the course of our numerical investigation, we have found that when *Q* is small (around two or three times the vortex ring's core radius), the resulting dynamics can be quite different. For example, instead of leapfrogging, the rings can form a bound pair and travel together maintaining their radii and common speed. While such dynamical features are intriguing in their own right, for the scope of the current work, we only use results from vortex rings separated with large enough values of *Q* to result in leapfrogging dynamics.

In order to formulate a reduced model of the leapfrogging interaction of the rings, we start with reviewing the situation in the classical fluid case. There, a vorticity field ** ω** will induce a velocity field according to the Biot-Savart law

^{71,72}

where the integral is computed over all of space. When considering a vorticity distribution of infinite strength but having a constant circulation κ along a closed curve **R** parametrized by ℓ, Eq. (8) becomes

where **s** is the unit tangent vector. Equation (9) can be evaluated to find the effect of a vortex ring on another vortex ring. Incorporating the single vortex ring kinematics of Eq. (2) and applying the Biot-Savart law to *N* co-axial rings results in the following reduced ordinary differential equations (ODEs) for their effective motion:^{68,69,73,74}

where

and

*R*_{i} and *Z*_{i} are the radius and *z*-coordinate of the *i*th ring [see Fig. 4(b)], κ_{i} is its circulation or “vortex strength” which, in our case, is a constant for each ring, *r*_{c, i} is its core radius which is time-dependent and follows the relation

*C*is a constant, which for the case of a classical fluid is given as

*C*= 1/4.

We note that the terms involving *U* relate to the interaction of the rings, while the first term of Eq. (10) represents the innate translational velocity of each ring. If one compares this velocity with that of the quantum vortex ring velocity of Eq. (2), we see that they are quite similar. We suggest that it is possible to use the results of Eqs. (10) and (11) for the vortex rings in the NLSE by simply substituting the quantum fluid values such that *C* = *L*_{0}(*m*) − 1, κ = 4π*am*, and *r*_{c} from Eq. (3) into Eqs. (10) and (11). We point out, that in contrast to the classical fluid leapfrogging picture, in a quantum fluid, the healing length being constant forces *r*_{c} to be, approximately, time-independent. The conservation of the core width during vortex ring evolution will be briefly addressed below.

In Figure 8, we show a series of (*R*, *Z*) orbits—in the center-of-mass frame—as predicted by Eqs. (10) and (11) for different initial separation distances along with those obtained through tracking the vortex rings in the integration of the full 3D NLSE system (1). As it is clear from the figure, the (*R*, *Z*) orbits from the reduced ODE system are in very good agreement with the ones obtained from the full PDE. These results stress the usefulness of the reduced system towards capturing the interaction dynamics of co-axial vortex rings. We note that while the ODE trajectories are periodic to numerical accuracy, the PDE results are only approximately periodic (only one “period” is shown here). The deviation from periodicity, in the co-moving reference frame, in the full PDE numerics might be attributed to several factors including: internal mode (Kelvin) excitations and stability features of the actual leapfrogging orbits.^{75}

Figure 8 suggests the existence of a stable, co-moving, equilibrium at the center of the periodic orbits for *Z* = 0 and *R* being constant. This point corresponds precisely to two overlapping (*Z*_{1} = *Z*_{2}) vortex rings of the same radius, namely, a vortex ring of charge two. In fact, the linearization about this co-moving equilibrium would yield an approximation for the period of the vortex ring leapfrogging in the limit of small initial separations. However, as vortex rings of higher charge are unstable and break up into two unit charge vortex rings, we will not study them further here—with the exception of using the theoretically predicted velocity of the doubly charged vortex ring as an approximation for the average leapfrogging velocity.

It is interesting to note that the vortex core size in the NLSE (1) should be approximately constant provided the vortex core is not significantly perturbed. In fact, if the vortex core size is invariant during evolution, then the total length of the vorticity tube should be conserved due to conservation of angular momentum. Therefore, by measuring the total vortex tube length, it is possible to monitor the degree to which the vortex core width deviates from its unperturbed value. In the case of leapfrogging vortex rings, where one vortex is forced inside the other ring, there is a strong mutual perturbation that can affect the vortex core width. In Fig. 9, we depict the sum of the vortex radii for different case examples of two vortex rings undergoing leapfrogging for both the effective ODE system and the original PDE. The sum *R*_{1} + *R*_{2} is proportional to the total vortex tube length *L* = 2π(*R*_{1} + *R*_{2}). As it is clear from the figure, the total vortex length is not conserved but it oscillates around a mean value with a relatively small amplitude. For the effective ODE model, the total vortex tube length varies between 4% and 12%, while for the full 3D model it varies between 3% and 8% for the cases depicted in Fig. 9. This result is a clear indication that the vortex core is indeed perturbed by the strong mutual interaction between vortex rings during the leapfrogging dynamics.

It is also worth mentioning that the results in Figs. 8 and 9 suggest that although the reduced ODEs are able to capture very well the shape of the leapfrogging orbits in the (*R*, *Z*) plane, there appears to be a systematic offset (between the ODEs of Eqs. (10) and (11) and the full 3D NLSE) as regards the period for these orbits. Figure 9 shows that for small initial vortex ring separations [see panel (a)] the offset in the period is larger than for large initial vortex ring separations [see panel (c)]. This can be qualitatively understood from the fact that an initial small separation between vortices induces stronger interactions between vortex rings bringing them closer to each other and thus strongly perturbing the vortex core which was assumed to have zero width in the Biot-Savart law yielding the ODE approximation.

Two relevant quantities of the leapfrogging rings that we can use to further validate the reduced model are the average ring-pair *z*-velocity and the rings' period of oscillation. We measure the average ring-pair *z*-velocity by tracking the *z* value of the center position of a (*y*, *z*) cut of the rings over time and then perform a least-squares linear fit. This is done for both the full NLSE simulations, and for the ODE integration of the reduced system. In Fig. 10, the results are presented as a function of *d* for a few values of the initial *z* separation distance, *Q*. We observe that the average velocity of the leapfrogging rings obtained in the PDE model is always lower than that obtained in the effective ODE model, the percent difference ranges between 1.6% and 7.2%. It is interesting to compare the average velocity of two leapfrogging rings with the velocity predicted for a vortex ring of charge *m* = 2 from Eq. (2). As depicted by the black thin line in the top panel of Fig. 10, the velocity of a vortex ring of charge *m* = 2 (evaluated at the initial radii for the leapfrogging vortex ring pair) closely resembles the average velocity of a vortex ring pair under leapfrogging evolution. This stems from the fact that, in the far field, a vortex ring pair can be approximated by a vortex ring of charge *m* = 2. In fact, after closer inspection of the results in the top panels of Fig. 10, the velocity of the doubly charged vortex ring is closer to the velocity for the vortex pairs when their separation is small, i.e., when they are closer to each other and thus closer to a doubly charged vortex ring. While, this relation could be further explored by taking the average radius along one leapfrogging period instead of the initial radii, the main phenomenology does not change significantly.

In the bottom panel of Fig. 10, we show the results of measuring the period of oscillation of the leapfrogging rings from both the full NLSE simulations and the ODE reduced model for the same ring separations and radii. We see that the period of oscillation for the PDE model is always longer than that of the ODE model, and the percent difference in period ranges from 4.2% to 28.5%, which, while clearly worse than the velocity comparison, is not unreasonable. It is interesting to note that while the period of the leapfrogging rings changes depending on the *z*-distance between the rings, it remains fairly constant for vortex ring pairs of *different* radii at the same separation distance.

From the velocity and period results, we see that there are obviously some differences between the quantum and classical fluid regimes including the modified values that result in lower average pair velocities and longer periods of oscillation than those predicted by our current reduced equations of motion. These are worthy of further investigation, in order to improve the accuracy of the reduced equations of motion. Nevertheless, the trends of the two sets of results are very similar, and their magnitudes are comparable, which suggests that our effective description provides a valuable means of examining the numerical (and, presumably, also the experimental) data and, as such, is useful in its own right as well as worthy of further investigation.

Finally, it is worth pointing out that there have been recent efforts in describing the leapfrogging dynamics of vortex ring bundles in classical fluids using the vortex filament method.^{75} In particular, it is found that the leapfrogging dynamics is weakly unstable for a small number of vortex rings and progressively more unstable for larger numbers of vortex rings. In our numerics, we have found similar instabilities, that depending on the separation and radii of the vortex rings, may not be visible until tens of leapfrogging periods. Therefore, although weakly unstable, the time required for the instability to manifest itself might be long enough to be able to observe leapfrogging during the typical lifetime of current BEC experiments.

It is also important to mention that the 2D equivalent of leapfrogging has been studied in the case of point vortices.^{76,77} This 2D leapfrogging corresponds to the dynamics of two vortex dipoles (vortex and anti-vortex pairs) and it is found to display three distinct regions depending on the distance between pairs (relative to the internal distance between vortices in each pair) as follows. (a) The small separation region does not yield leapfrogging. (b) Intermediate separations yield *unstable* leapfrogging. (c) Large separation distances yield *stable* leapfrogging. The 2D vortex dipole leapfrogging is strongly connected to 3D vortex ring leapfrogging since any 2D cut passing through the axes of the vortex rings will produce a 2D pair of vortex dipoles. Therefore, 2D vortex dipole dynamics already encompasses the fundamental interactions responsible for vortex ring leapfrogging. Therefore, it would be interesting to test if the vortex ring leapfrogging also displays the three regions described above and whether the weak instability that we have observed can in fact be eliminated by using larger vortex ring radii. This avenue of research is currently being examined and a systematic investigation will be reported in a future publication.

## V. CONCLUSIONS

We studied the dynamics and interactions of vortex rings in the nonlinear Schrödinger equation as a model for a homogeneous Bose-Einstein condensate. Upon corroborating the well-known single vortex ring evolution in our 3D NLSE framework, we focused on two dynamical scenarios: scattering and leapfrogging. For the former, we considered the scattering collision of two co-planar, opposite charge, vortex rings. Upon their collision, these scatter at an angle depending on the initial conditions. We found a surprisingly simple phenomenological rule for the scattering angle as a function of the vortex ring radii and the initial distance (impact parameter) between the vortex rings. For the second scenario, we followed the dynamics of a pair of co-axial vortex rings of the same charge. We proposed a modification of the effective equations of motion for the vortex ring interactions in a classical superfluid. The modification was based on correcting the self-induced velocity using previous results for NLSE settings in a quantum superfluid. We compared the resulting effective equations of motion against full 3D simulations of the original NLSE model for several cases of leapfrogging evolution and found good agreement between the two. Finally, we found a monotonic decrease for the average velocity of the leapfrogging as the radii of the ring increase and/or as their initial separation increases. We also found that the leapfrogging period remains approximately constant for vortex ring pairs of different radii starting at the same separation distance, while the period increases with initial separation distance for fixed radii.

Future work will be directed towards refining the reduced equations of motion to better account for the differences between classical and quantum fluids (especially in the context of the term encoding the interaction between the different rings). By directly using the Biot-Savart law, it should be possible to obtain a more general reduced model to describe more complex vortex ring configurations including non-co-axial vortex ring interactions and potentially also incorporating features such as the presence of Kelvin waves. Such an approach would lie in between the full NLSE simulation and the ODEs representing the radius and center position of the rings utilized herein and would thus facilitate the comparison between the two, as well as an understanding of the different approximations involved in obtaining the latter from the former. It would also be interesting to understand and quantify the instability (and its range of relevance over initial conditions and system parameters) observed in the leapfrogging dynamics, not only for a pair of vortex rings but also, for vortex ring bundles.^{75} Another avenue for future exploration consists of using classical fluid-inspired effective equations of motion (Biot-Savart law) to obtain reduced ODEs describing the effects of slow-down created by internal excitation (Kelvin) modes^{43,44} and, in more general terms, the effects of these modes on the interactions and scattering between vortex rings. We are currently exploring these directions and will report on them in future works.

## ACKNOWLEDGMENTS

We gratefully acknowledge the support of NSF-DMS-0806762 and NSF-DMS-1312856, NSF-CMMI-1000337, as well as from the (U.S.) Air Force Office of Scientific Research (USAFOSR) under Grant No. FA950-12-1-0332, the Binational Science Foundation under Grant No. 2010239, from the Alexander von Humboldt Foundation and the FP7, Marie Curie Actions, People, International Research Staff Exchange Scheme (IRSES-606096). R.M.C. and J.D.T. acknowledge the support of the Computational Science Research Center (CSRC) at SDSU.

### APPENDIX: COMPUTATION OF THE VORTEX CORE PARAMETER

The “vortex core parameter” *L*_{0}(*m*) is defined as the convergent part of the energy per unit length of a vortex line. In this appendix, we simply state the equations used to compute *L*_{0}(*m*). For full details of the derivation of the equations, we refer the reader to Ref. 57.

For the NLSE in the form of Eq. (1), the vortex core parameter is found by numerically integrating

where *f*(*r*) is found numerically by solving the ODE

and where we note that for a given *m*, *L*_{0}(*m*) is invariant over the parameters *a*, *s*, and Ω.

We compute the integrals in Eq. (A1) using a simple trapezoid integration on a numerically exact *f*(*r*) computed on a domain of *r* = [0, 500], with a spatial step size of Δ*r* = 0.0075 and NLSE parameters *a* = 1, *s* = −1, and Ω = −1. The values of *L*_{0} for charges *m* = 1, …, 10 are given in Table I. We note that the values of *L*_{0} for |*m*| > 3 have not been previously reported, and those for *m* = 2 and *m* = 3 are given with more accuracy than reported in Ref. 54.

m
. | L_{0}(m)
. | m
. | L_{0}(m)
. |
---|---|---|---|

1 | 0.380868 | 6 | 0.022954 |

2 | 0.133837 | 7 | 0.017785 |

3 | 0.070755 | 8 | 0.014247 |

4 | 0.044567 | 9 | 0.011713 |

5 | 0.030981 | 10 | 0.009833 |

m
. | L_{0}(m)
. | m
. | L_{0}(m)
. |
---|---|---|---|

1 | 0.380868 | 6 | 0.022954 |

2 | 0.133837 | 7 | 0.017785 |

3 | 0.070755 | 8 | 0.014247 |

4 | 0.044567 | 9 | 0.011713 |

5 | 0.030981 | 10 | 0.009833 |

It is intriguing to note that one can produce a very close power-law fit yielding an approximation of *L*_{0}(*m*) given by