Classical and quantum vortex leapfrogging in two-dimensional channels

The leapfrogging of coaxial vortex rings is a famous effect which has been noticed since the times of Helmholtz. Recent advances in ultra-cold atomic gases show that the effect can now be studied in quantum fluids. The strong confinement which characterizes these systems motivates the study of leapfrogging of vortices within narrow channels. Using the two-dimensional point vortex model, we show that in the constrained geometry of a two-dimensional channel the dynamics is richer than in an unbounded domain: alongsize the known regimes of standard leapfrogging and the absence of it, we identify new regimes of backward leapfrogging and periodic orbits. Moreover, by solving the Gross-Pitaevskii equation for a Bose-Einstein condensate, we show that all four regimes exist for quantum vortices too. Finally, we discuss the differences between classical and quantum vortex leapfrogging which appear when the quantum healing length becomes significant compared to the vortex separation or the channel size, and when, due to high velocity, compressibility effects in the condensate becomes significant.


I. INTRODUCTION
The leapfrogging of two co-axial vortex rings (in three dimensions) or of two vortex-antivortex pairs (in two dimensions) is a benchmark problem of vortex interaction [41] which dates back to [24]. The time evolution of this vortex configuration is striking: the vortex ring (or pair) which is ahead widens and slows down, while the ring behind contracts, speeds up, catches up with the first ring and goes ahead through it; this 'leapfrogging' game is then repeated over and over again, unless instabilities disrupt it. A number of papers have been written on different aspects of this problem, ranging from the stability [1,26,36,64] to the deformation of the vortex cores and to the effects of viscosity [60] using numerical [10,13,54] as well as experimental methods [37,40,51,73]. The most recent developments concern leapfrogging of vortex bundles [68] and helical waves [27,52,57].
Our work is motivated by recent experiments with atomic Bose-Einstein condensates, which constitute a dilute quantum fluid and provide an idealised platform to study fundamental vortex dynamics [70]. In these experiments, atomic gases are confined by suitable magnetic-optical traps and cooled to nano-Kelvin temperatures. If the atoms of the gas are bosons (i.e. have integer spin), a phase transition occurs upon cooling below a critical temperature T c , and the gas forms a macroscopic coherent quantum state [4] called a Bose-Einstein condensate (BEC). From the point of view of the hydrodynamics, a BEC has three key properties: it is superfluid (i.e. it suffers no viscous losses of kinetic energy when it flows), it is compressible, and its vorticity is concentrated to thin hollow vortex lines with fixed width a 0 and fixed circulation ±h/m where h is Planck's constant and m is the mass of a boson (while vortices with larger quanta of circulation, ±2h/m, ±3h/m, · · · , are possible, they are unstable to decay into multiple singly-charged vortices). Thus, in BECs, vortices are well-defined and identical objects, evolving in an inviscid compressible fluid.
There are several additional characteristics of atomic BECs that make them attractive for probing vortex dynamics. Firstly, the physical parameters of the fluid (including the width and speed of the vortices) are tunable, for example, through the density of the gas and the strength of the atom-atom interaction (which can be modified by means of Feshbach resonances [29]); this should be contrasted with superfluid liquid helium -historically the most studied quantum fluid -whose physical parameters are fixed by nature. Secondly, the potential experienced by the gas can be controlled through magnetic and optical fields. Such trapping is essential, on one hand, to contain the gas, and gives rise to the boundary effects which are central to this work. However, the potential can also be exploited to engineer the dimensionality of the gas -particularly, quasi-two-dimensional geometries in which vortex lines effectively become point-like vortices -and to stir and shake the condensate. Finally, recent techniques have enabled the observation of vortex lines [59] and vortex points [58] in real-time, including inference of their individual circulations.
Atomic BECs however are characterized by their small dimensions, typically from 10 to 100 times the vortex core size, for which the motion of vortices can be significantly affected by the presence of boundaries. This drawback is mainly due to the loss of atoms in the final evaporative stage of cooling the gas. There are even experiments in which, by design, the most interesting physics occurs in the most restricted region of the system, for example vortex rings nucleated in the weak link of the Josephson junction between two condensates [66,72]. The aim of the present work is to provide insight in the interpretation of current and future experimental studies of vortex dynamics in confined condensates (rather than idealised open domains), where leapfrogging dynamics, which can be established if the vortex nucleation frequency is sufficiently high, is affected by the presence of boundaries. The characteristics of leapfrogging motion in such confined systems is likely to show significant dissimilarities compared to the corresponding dynamics in unbounded systems stemming from the role played by image vortices arising from the presence of boundaries. Despite the expected impact of geometrical confinement, to the best of our knowledge the role of boundaries in leapfrogging dynamics has never been investigated in literature neither for classical nor for quantum fluids ( [31] and [28] indeed studied leapfrogging in homogeneous condensates, without boundaries). In order to assess the impact of the boundaries and disentangle the latter from other concurrent physical effects existing in quantum fluids (e.g. compressibility), in this research we compare the leapfrogging of vortices in plane channels in (i) ideal incompressible classical fluids and (ii) box-trapped Bose-Einstein condensates. In order to simplify the system under investigation, our theoretical and numerical analysis is performed in two-dimensions, employing the point vortex model for classical fluids and the Gross-Pitaevskii equation for BECs. We stress that the Gross-Pitaevskii equation has proved an excellent quantitative model of experiments with Bose-Einstein condensates at temperatures T T c ; at relatively high values of temperature, the condensate exchanges energy and particles with the thermal cloud, and the Gross-Pitaevskii equation requires modifications [5,7,8,50]. We also remark that on one hand the two-dimensional nature of the system that we consider is an idealisation (the aim is to get insight in the motion of three-dimensional vortex rings), but, on the other hand, where atomic Bose-Einstein condensates are tightly confined in one direction the system becomes effectively two-dimensional and our two-dimensional approach becomes realistic.
The article is organised as follows. In Section II, we illustrate the two theoretical models employed, namely the classical point vortex model and the Gross-Pitaevskii equation describing the dynamics of BECs in the zerotemperature limit. In Section III, we report the results obtained in both classical and quantum fluids, focusing on the role of boundaries and on the differences between classical and quantum systems. Finally, in the last Section IV, we summarise our findings and illustrate their importance in the future of quantum vortex experiments.

II. MODELS
A. Point vortex model The simplest model of our system is the classical point vortex model: a two-dimensional inviscid incompressible irrotational fluid in an infinite channel of width 2D containing two vortex-antivortex pairs (the two-dimensional analog of three-dimensional coaxial vortex rings), each of circulation ±Γ. In view of comparing the results obtained with this classical model to quantum vortices in confined BECs, the hypotheses behind the point vortex model must be carefully considered.
The classical model describes a fluid with constant density. In the bulk of the condensate, i.e. sufficiently far from boundaries or vortices, this assumption is realistic: indeed, although in past experiments condensates were usually confined by harmonic trapping potentials resulting in density gradients [14], current experimental techniques [20] allow box-like trapping potentials which lead to uniform density profiles in the bulk of the condensate as in the classical point vortex model. In particular, in the vicinity of a vortex, the classical model assumes constant density at any radial distance r to the vortex axis, including the vortex axis r = 0 itself. In Bose-Einstein condensates, a vortex is a topological defect of the phase of the governing complex wavefunction (or order parameter), as we shall describe with more details in Section II B 1. Therefore the vortex core is a thin tubular region around the vortex axis which is depleted of atoms: as r → 0, the velocity tends to infinity, as in the point vortex model, but the fluid density tends to zero. The radius of this tube is of the order of the quantum mechanical healing length ξ (see Section II B 1). A similar difference between the classical point vortex model and Bose-Einstein condensates occurs near a hard boundary: the classical model assumes that the fluid's density is constant up to the boundary; in a Bose-Einstein condensate a thin boundary region (again of the order of ξ) forms near the boundary where, in the case of box-like traps, the condensate's density rapidly drops from the bulk value to zero. We conclude that, from a geometrical point of view, the classical point vortex model can be used to model Bose-Einstein condensates provided that vortex-vortex and vortex-boundary distances are larger than the healing length ξ.
From a dynamical point of view, the assumption of constant density implies that the classical point vortex model neglects sound waves which are radiated away by quantum vortices when they accelerate [3]. The point vortex model, in fact, is based on the classical ideal Euler equation which conserves energy. In the low temperature limit T /T c 1 that guarantees the validity of the Gross-Pitaevskii equation, the total energy of a BEC is constant, but transformation of incompressible kinetic energy of the vortex configuration into compressible kinetic energy of the field of sound waves (or vice versa) is permitted. This dynamical difference between the classical point vortex model and the Gross-Pitaevskii model is, physically, perhaps the most significant, and will be addressed while discussing the results in Section III.
Despite these approximations, we reckon the model captures the essential ingredient of our problem: the motion of quantised irrotational vortices in presence of boundaries. Indeed, the classical point vortex model in a circular disk has been already used with success to model two-dimensional turbulence in low temperature trapped condensates, for example, by Simula et al. [61]. It must also be noticed that Mason et al. [38] have shown that the motion of a realistic vortex at distance d to a boundary can be described in terms of a classical image vortex even if ξ is comparable to d (although a small correction is needed to account for the density depletion in the boundary region). In the suitable physical limits, we hence expect the point vortex model to correctly describe the impact of boundaries on the leapfrogging of quantised vortices.

Equations of motion
Our physical domain under investigation is a two-dimensional infinite strip C R 2 defined as C R 2 = {(x, y) ∈ R 2 : (x, y) ∈ (−∞, ∞) × (0, 2D)}, which hereafter we will refer to as the channel. We assume the flow to be twodimensional, i.e. the velocity vertical component v z = 0 and the horizontal components v x and v y only depend on horizontal coordinates x and y and time t. The incompressibility assumption implies that the continuity equation can be written as follows (1) The velocity field v can hence be expressed as the curl of vector field Ψ which, given the two-dimensionality of the flow, has non-vanishing components only in the z direction, Ψ = (0, 0, ψ(x, y, t)). The velocity components have hence the following expressions in terms of the function ψ which is often denominated The irrotationality of the flow implies that the velocity field can be expressed via a potential function ϕ, i.e.
leading to the following relations for the components v i = ∂ i ϕ. Equations (1) and (2) imply that both ϕ and ψ satisfy Laplace equation, ∆ϕ = ∆ψ = 0, and the following equalities between their spatial derivatives: Equations (3) and (4) coincide with the well known Cauchy-Riemann relations for the complex function Ω(z) := ϕ+iψ, where z = x + iy. Hence, following basic complex analysis, the function Ω(z), denominated complex potential, is an analytical complex function on the simply connected open domain C = {z ∈ C : 0 < m z < 2D} C. As a consequence, Ω(z) is differentiable and its derivative is the so-called complex velocity. In the framework of complex potentials, the impermeable boundary conditions for ideal fluids correspond in our channel C to the following constraint: The description of incompressible and irrotational flows of ideal fluids via the complex potential-based formulation is particularly useful in the present work as it allows the employment of conformal mapping techniques for the derivation of the analytical expression of the complex potential Ω(z) describing the velocity field induced by a point vortex in our channel C. The essential steps for this derivation are as follows. The necessary ingredients are mainly two: (a) the knowledge of the complex potential Θ(ζ) describing the flow induced by a point vortex in a simply connected open subset D of the complex plane, with ζ ∈ D C; and (b) the construction of a conformal map ζ = f (z) transforming our channel C onto the domain D.
Conformal maps f are transformations defined on the complex plane which preserve angles. Such maps are performed by analytical complex functions with non-vanishing derivative, i.e., in the present case, f (z) = 0 for all z ∈ C. The requirement D not coinciding with the entire complex plane C, is fundamental in order to exploit the Riemann Mapping Theorem which ensures the existence of the conformal map f mapping C onto D. Once Θ(ζ) and f (z) are determined, the complex potential Ω(z) for a vortex flow in C is obtained by transforming the potential Θ(ζ) via the conformal map f −1 (ζ), i.e.
The reasons why the derived complex function Ω(z) via Eq. (6) is the seeked complex potential are the following. First, Ω(z) is analytic on C (as it is obtained via the composition of two analytic functions, f and Θ), implying that the real and imaginary parts of Ω(z) are related to each other via the Cauchy-Riemann equations and are both harmonic functions. Hence, they do satisfy all the necessary conditions for corresponding respectively to a potential function and a streamfunction of an incompressible and irrotational flow of an inviscid fluid. Second, the correspondence of ∂C and ∂D under the conformal mapping performed by f transposes the boundary conditions enforced by Θ(ζ) on ∂D to the boundary ∂C [35]. Finally, via conformal mappings, the flow induced by a vortex of circulation κ is indeed mapped to a vortex flow with the same circulation [45].
In the present work, we choose D to coincide with the upper half complex plane, i.e. D = {ζ ∈ C : m ζ > 0}. In this domain, the complex potential Θ(ζ) describing the flow induced by a vortex placed in ζ 0 ∈ D is obtained by the method of images, namely where sgn(ζ 0 ) is the sign of the vortex placed in ζ 0 (positive for anti-clockwise induced flow, negative for clockwise), ζ * 0 is the complex conjugate of ζ 0 where a vortex of opposite sign is placed (the image-vortex of ζ 0 ) and κ is the circulation of the flow generated by the vortex. The analytical function f transforming conformally the channel C = {z ∈ C : 0 < m z < 2D} onto D is as follows (see Fig. 1 for a schematic illustration) The Employing Eq. (6), the determination of the complex potential Ω(z) is straightforward, namely leading to the following complex velocity where The complex function χ(z, z 0 ) (and, correspondingly, χ(z, z * 0 )) can be physically interpreted as the complex velocity generated by an isolated vortex placed in z 0 (whose complex potential would be Ω(z, z 0 ) = −sgn(z 0 )iκ log(z−z 0 )/(2π)) and its infinite images with respect to the walls of the channel, m z = 0 and m z = 2D. The expression (10) for the complex potential w(z, z 0 ) can indeed be derived by considering two sets of infinite images of a vortex placed in z 0 and an anti-vortex in z * 0 [22]. If the channel is characterised by the presence of N vortices, the complex velocity w(z, z k {k=1,...,N } ) generated by the the set of N vortices is obtained via the superposition principle, i.e.
A crucial role in this N -vortex problem is played by the equations of motion of a generic j-th vortex. In order to derive such equations of motions, we define the position z j (t) := x j (t) + iy j (t) occupied by the vortex at time t in the channel C. Indicating with the superscript '˙' derivation with respect to time, we define the quantityż j (t) :=ẋ j (t) + iẏ j (t), where the real and imaginary part correspond to the x and y components of the j-th vortex velocity. As vortices are advected by the local fluid velocity, i.e.ẋ j (t) = v(x j (t), t), the following relation holdṡ where we have omitted the time dependence of z j and z k to ease notation and the complex conjugation on the r.h.s. arises from the definition (5) of complex velocity. In order to determine the complex velocity w(z j , z k {k=1,...,N } ), we employ Eq. (11) subtracting the term corresponding to the vortex placed in z j , obtaining the following relatioṅ which coincides with the equations of motion of the j-th vortex. The equations of motion of the whole N -vortex problem are hence a set of 2N coupled ordinary differential equations.

B. Gross-Pitaevskii equation model
The Gross-Pitaevskii model is a well-established theoretical framework for the investigation of the dynamics of BECs at temperatures much smaller than the critical transition temperature. The Gross-Pitaevskii (GP) equation describes the temporal evolution of the complex order parameter Ψ = Ψ(x, t) of the system, and reads as follows, where the dot is the time derivative, = h/(2π) is the reduced Planck's constant, m is the boson mass, is an externally applied potential, and g = 4π 2 a s /m models the two-body contact-like boson interaction, where a s is the s-wave scattering length for the collision of two bosons. The order parameter Ψ can be written in terms of its amplitude and its phase as where n = n(x, t) = |Ψ| 2 is the particle number density (number of bosons per unit volume) and θ = θ(x, t) is the phase. Without loss of generality, the order parameter Ψ can be written as Ψ(x, t) = e iµt/ Φ(x, t) where µ is called the chemical potential and Φ(x, t) obeys

Quantum vortices
In the context of BECs described by the Gross-Pitaevskii equation, quantum vortices are topological defects of the phase θ of the order parameter, at which Ψ = 0 (hence θ is undefined) and around which θ wraps by 2qπ with q ∈ Z \ {0}. In three dimensions, vortices take the form of one-dimensional curves which may form a vortex tangle, as observed both in BECs [69] and superfluid helium [67]. In two dimensions, vortices coincide with vortex points which have been observed extensively in oblate (pancake-like) BECs [39]. For the purpose of the present work, we will restrict our discussion to two dimensional systems.
The velocity field v(x, t) associated to a BECs whose dynamics is described by the order parameter Ψ, is defined from the phase θ via the relation v(x, t) = m ∇θ.
Employing the definition (17) of the velocity and the 2qπ phase wrapping existing around a vortex, it is straightforward to verify that the circulation Γ of the velocity field on any closed curve γ enclosing a vortex point is quantised in terms of the quantum of circulation κ = h/m, i.e.
Choosing γ to be a circle of radius r and assuming the flow around a vortex to be axisymmetric, the azimuthal component of the flow velocity around a vortex is given by the relation v φ = qκ/(2πr), coinciding with the expression for a classical point vortex. Hence, from a velocity point of view, quantum and classical vortices are identical. The important and dynamically significant distinction between classical and quantum vortices is that the latter are characterised by a finite core whose size is of the order of the so-called healing length ξ = / √ mgn. As we will very briefly illustrate in the next section, quantum fluids are indeed compressible fluids.

Fluid dynamical equations for a BEC
The Gross-Pitaevskii equation (14) may be rewritten via the Madelung transformation consisting in expressing Ψ in polar form (15) and separating the real and imaginary parts of (14). This procedure leads to the following equationṡ where p and p are respectively pressure and quantum pressure n∆ (ln (n)) .
Equation (19) coincides formally with the continuity equation of a classical fluid, while equation (20), exception made for the presence of the quantum pressure p , is formally identical to the momentum balance equation for a barotropic, compressible classical Euler (ideal) fluid. At length scales much larger than the healing length ξ (which is the typical length scale for density variations, associated e.g. to the presence of vortices or boundaries) p /p 1, implying that in this limit the BEC can indeed be considered as a barotropic, compressible classical inviscid fluid. Hence, at length scales ξ, the dynamics of quantum and classical point vortices only differ on the basis of compressible phenomena which may arise in BECs. In the other limit of ∼ ξ, the physics may be significantly different. For instance, if the relative distance between quantum vortices of opposite sign is of the order of ξ, the quantum pressure term would trigger the annihilation of the vortex pair, while in the classical point vortex model no loss of circulation is included in the model. Moreover, the behaviour of a co-rotating pair of quantum vortices of same sign also shows dissimilarities with respect to the classical case, in particular for the finite value of the rotation frequency ω τ as the distance tends to zero (in the classical model, the frequency diverges, ω τ ∼ 1/ 2 ).

A. Classical fluids
To make progress in understanding the impact of boundaries on the leapfrogging behaviour of classical point vortices in a two-dimensional channel, we consider the motion of four vortices, half with positive circulation κ, half with negative −κ. In Fig. 2 we show this initial condition. If we interpret our two-dimensional configuration as a model of a three-dimensional configuration of vortex rings, point vortices of same colour in the figure correspond to cross-sections of the same ring. Initially, the four vortices are vertically aligned on the y axis, i.e. x j (0) = 0 for j = 1, . . . , 4 and the vortex-anti vortex pairs are symmetrically positioned with respect to the channel mid axis y = D, namely y j (0) = D ± R for the first pair j = (1, 2) and y j (0) = D ± r for the second pair j = (3, 4), with the conditions R/D < 1 and r/R < 1. In order to characterise the dependence of vortex trajectories on the two non-dimensional parameters r/R and R/D which determine the flow, we numerically integrate the equations of motion (13) for the four vortices, j = 1, . . . , 4, varying r/R and R/D. In particular, we choose r/R = n/10 and R/D = m/10, with m, n = 1, . . . , 9. The time-advancement scheme employed in the numerical simulations is a second-order Adams-Bashforth method with a time step ∆t = T /1000 where T = 2π 2 δ 2 /κ is the rotation period of a pair of vortices of the same polarity placed at distance δ. In our numerical simulations δ is set to 10 −3 D.
For classical unbounded fluids, since the study performed by Love over a century ago [36], it is well known that vortices undergo leapfrogging motion only if r/R is larger than a critical value α c = 3 − 2 √ 2 ≈ 0.172. If r/R < α c , leapfrogging does not occur: the smaller, faster pair moves "too fast" for the larger ring to influence its dynamics in a significant way, and the vortices separate. More recently, Acheson [1] extended numerically the study performed by Love and established that leapfrogging motion is unstable when α c < r/R < α c , with α c = 0.382.
In our two-dimensional channel, the confinement of the flow leads to a richer dynamics than in an unbounded domain. In addition to distinction between leapfrogging and non-leapfrogging, which is already known, we also observe backward leapfrogging and periodic orbits. The phase diagram of the system resulting from the numerical simulations is illustrated in Fig. 3.
For values of R/D ≤ 1/2, the dynamics is very similar to what is observed in an unbounded fluid, the role of the boundaries being only marginal. For a given value of R/D ≤ 1/2, in fact, as we increase r/R, we first observe non leapfrogging motion (in black in Fig. 3), defined as the dynamics characterised byẏ j (t) = 0 for all j at late times; then we notice unstable leapfrogging motion (open red squares), and finally stable leapfrogging (filled red squares). These dynamical regimes therefore coincide with the scenario outlined by Acheson [1], the only significant and important difference being the dependence of α c on R/D: for small values of R/D, α c is very close to the constant value 0.172 for vortex leapfrogging in unbounded fluids (e.g. for R/D = 0.1, α c = 0.173), increasing for increasing values of R/D (e.g. α c = 0.216 for R/D = 0.5). This dependence of α c on R/D stems from the interaction of the outer vortices 1 and 2 in Fig. 2) with their corresponding images with respect to the closest channel wall; essentially, the interaction with image vortices is stronger compared to the interaction of the inner pair with the corresponding images. These images, of opposite sign, slow down the outer vortex pair, allowing the inner pair to escape towards infinity for values of r/R which would produce leapfrogging motion in an unbounded fluid; in order to recover leapfrogging, r/R would have to increase. As R/D increases, this effect is amplified as the outer pair is closer to the channel walls.
This increasing monotonous behaviour of α c with respect to R/D extends also for R/D > 1/2, where the role played by boundaries becomes significant, triggering a much richer dynamics. As R/D is larger than 1/2, for large values of r/R, we observe backward-leapfrogging, indicated by blue diamonds in Fig. 3. This dynamics, again, originates from the interaction of vortices with their images with respect to the closest channel wall. In particular, each vortex, paired to its image of opposite sign, forms a virtual vortex-anti vortex pair on its own. As a consequence, we observe two distinct leapfrogging motions, each involving two virtual vortex-anti vortex pairs. Due to the vortex polarity, the leapfrogging motion induces a net translation in the opposite direction with respect to standard (forward) leapfrogging. In the (R/D, r/R) plane, the forward-leapfrogging to backward-leapfrogging transition occurs via an intermediate regime in which vortices follow periodic orbits, indicated by green stars in Fig. 3. As shown in detail in the next section and in the analytical derivation presented in the Appendix, periodic orbits are observed when R+r = D, corresponding to the green dashed line in Fig. 3. For large values of R/D (R/D 3/4), the system crosses directly the no-leapfrogging to backward-leapfrogging boundary without passing through a forward-leapfrogging regime. Examples of all the different regimes observed in our system of classical point vortices are shown in Fig. 4. Note that in the three-dimensional coaxial vortex rings analog, vortices of the same colour correspond to cross-sections of the same vortex ring.

Derivation of periodic orbits
In this section we derive theoretically the existence of periodic orbits in the leapfrogging motion of four vortices in a channel using the classical point-vortex model. We show that under suitable conditions, namely when R + r = D, each pair of same signed vortices moves around a fixed point. Some analytic details are discussed in the Appendix.
We now consider the midpoint M between the vortex points P 1 and P 3 , namely z M (t) = x 0 (t) + x 1 (t) 2 + .
If we look for the conditions such that the velocity w(z M ) of the midpoint M is zero, we have Note that the same result Eq. (25) is found for the midpoint N between the two vortex points P 2 and P 4 . Since r, R, and D are positive real parameters, the only admissible values of k in (25) are k ∈ Z + . Moreover, we know that r < R < D, leading to r + R < 2D, which implies that the only admissible value for k is k = 1, i.e.
This is the most interesting result: it states that when the four vortices satisfy the condition (26) then the midpoints M and N are at rest: the two pairs of vortices (P 1 , P 3 ) and (P 2 , P 4 ) move hence symmetrically with respect to their correspondent midpoints, i.e.ẋ 0 (t) = −ẋ 1 (t) andṘ(t) = −ṙ(t). The last equality is fundamental as it expresses that if condition (26) is satisfied at a given t = t 0 , it will be satisfied for every t > t 0 . Thus, if the initial condition is prepared such that x 0 (0) = x 1 (0) = 0 and r(0) + R(0) = D, vortices will always move symmetrically with respect to their midpoints z M = i D 2 and z N = i 3D 2 .
The last step to demonstrate the existence of periodic orbits is to prove that the trajectories of the vortex points are closed curves rotating around the two midpoints M and N as, in principle more general trajectories with the restrictionṘ(t) = −ṙ(t) (for instance,Ṙ(t) =ṙ(t) = 0) could be possible, not leading to periodic orbits. We tackle this issue in the Appendix, to ease the readability of the manuscript.

B. Quantum fluids
The next step is to numerically probe the dynamical regimes of two quantum vortex-antivortex pairs interacting in a two-dimensional channel. We shall compare the results with the corresponding classical results outlined in the previous Section (III A).
We consider a two-dimensional BEC in a channel geometry, imprinting quantum vortices in the positions initially occupied by classical vortices. Note that, in addition to the parameters R, r and D already present in the classical point vortex formulation, in the Gross-Pitaevskii formulation of the problem we have an extra length scale -the healing length ξ -which plays a fundamental role in the dynamics. To assess the relevance of this extra length scale, we present numerical simulations of leapfrogging quantum vortices employing two distinct values of the channel half-width D: D 1 = 40ξ and D 2 = 20ξ. In order to model the channel confinement, we use the following potential V : corresponding to a channel of half-width D, where the density |Φ| 2 is constant everywhere with the exception of thin layer whose width is of the order of the healing length at the channel boundaries y = 0 and y = 2D. The trajectories of the quantum vortices are calculated as a function of time by numerically solving the equation of motion of the order parameter Φ, the dimensionless Gross-Pitaevskii equation Equation (28)  In these units the healing length and the bulk density in the channel are unity. The numerical integration of Eq. (28) is performed employing a fourth-order Runge-Kutta time advancement scheme and second-order finite differences to approximate spatial derivative operators. Time step ∆t/τ is set to 1.5 × 10 −2 and spatial discretization ∆x/ξ = ∆y/ξ is chosen to be equal to 0.25. In the set of simulations where D = D 1 = 40ξ, the numbers of grid-points in the x and y directions are N x = 6400 and N y = 400 respectively, leading to the computational box −800ξ ≤ x ≤ 800ξ and −10ξ ≤ y ≤ 90ξ. On the other hand, when D = D 2 = 20ξ, N x = 3200 and N y = 240 respectively, leading to the computational box −400ξ ≤ x ≤ 400ξ and −10ξ ≤ y ≤ 50ξ.
The initial imprinting of vortices is made by enforcing a uniform 2π phase wrapping around the positions employed as initial condition for the classical point vortex simulations and letting the system relax in imaginary time before starting the integration of Eq. (28) for t ∈ R. In Fig. 5 we report the density |Φ| 2 (x, y) (left) and the phase θ(x, y) (right) of the initial condition employed for R/D = 0.6 and r/R = 0.3 and D = D 1 = 40ξ. It can be easily observed that the density |Φ| 2 rapidly drops to zero at the vortex positions and outside the channel. Correspondingly, the four 2π phase wrappings can be distinguished in Fig. 5 (right).
To verify the existence in a BEC of all distinct regimes observed in the classical point vortex model (Section III A), we perform numerical simulations of quantum vortex leapfrogging along the vertical line R/D = 0.6 of the phasediagram reported in Fig. 3; we have chosen this value of R/D because along this line, as r/R varies from 1/10 to 9/10, all regimes which we have identified using the classical point vortex model are present.
The results are schematically outlined in Fig. 6, where classical vortex dynamics (left) is compared to quantum vortex dynamics at D = D 1 = 40ξ (middle) and D = D 2 = 20ξ (right). When D = D 1 , the boundaries of the phase diagram at R/D = 0.6 are at the same values of r/R in the classical and in the quantum case. When D = D 2 we observe two differences: first, periodic motion now occurs at (R/D, r/R) = (0.6 , 0.7) instead of (R/D, r/R) = (0.6 , 0.67); second, at (R/D, r/R) = (0.6 , 0.1) the internal vortex-anti vortex pair annihilates as their initial distance is only 2.4ξ. These differences at the smaller value of channel size are expected, as the healing length scale starts playing a role: only if D/ξ is sufficiently large we can expect classical and quantum dynamics to be the same. The exact matching of the observed dynamical regimes when comparing classical and quantum leapfrogging in a two-dimensional channel if D ≥ 40ξ is confirmed in Fig. 7, which shows the trajectories of quantum vortices for (R/D , r/R) pairs selected as for the classical trajectories illustrated in Fig. 4.
It is worth noting some minor differences between the quantum vortex trajectories and their classical counterparts reported in Fig. 4. Since the initial condition is not stationary with respect to any frame of reference, when we start integrating in time Eq. (28) for t ∈ R there is a sudden emission of sound waves, and as a result the entire vortex configuration is translated towards the positive x direction. The effect (which has been reported in the literature [16]), is visible in the top right, bottom left and bottom right panels of Fig. 7 when compared Fig. 4. In particular, this horizontal shift affects the periodic orbits reported in Fig. 7 (bottom right) whose center is slightly shifted towards positive x values. In addition, the number of periods observed in the x range [−2D, 2D] is different from the classical counterpart, possibly due to the compressible nature of a quantum Bose gas, in which incompressible kinetic energy may be transformed into compressible kinetic energy (sound) when vortices change their velocities (accelerate), as shown by Parker et al. [48], exactly as the accelerated motion of charged particles emit electro-magnetic radiation. The role played by this effective dissipation of kinetic energy into sound will be assessed in a future study.

IV. CONCLUSIONS
In conclusion, we have demonstrated that, in the confined space of a two dimensional channel, the classical problem of vortex leapfrogging acquires new aspects. Using the point vortex model we have found that, besides the known regimes of standard leapfrogging and absence of leapfrogging, there are two new regimes: backward leapfrogging and periodic motion. Using the Gross-Pitaevskii equation to model an atomic Bose-Einstein condensate (a compressible quantum fluid) confined within a channel, we have verified that all four regime also exist for quantum vortices. In large channels, the boundaries between these regimes are the same for classical and quantum vortices. Some differences appear if the channel size is reduced, and the finite-size nature of the quantum vortex core starts playing a role, or if the vortices are very close and sound radiation becomes important. The determination of a richer dynamics for the leapfrogging of vortices occuring in confined geometries will be particularly important for the interpretation and planning of ongoing and future experiments with atomic Bose-Einstein Condensates, where the dynamical regimes reported in the present work can be potentially observed.
Future work will address the problem in three dimensions, paying attention to the excitation of Kelvin waves along the vortex rings and the departure from axisymmetry.
By writingż 1 = −ρû ρ − ρθû θ , we then find the following equations forρ andθ: In order to prove that the trajectory of vortex P 1 is a closed curve, we need to show that the function ρ(θ) is a continuos and periodic function. However, the integration of equation (A5) is a hard task to achieve. Therefore, we choose to prove that ρ(θ) is a continuos and periodic function without finding the exact integral of (A5). In order to achieve this goal, we first need to recall a result from mathematical analysis, which states: Theorem 1 Given a continuos and periodic function f : R → R with period T such that Having recalled Theorem 1, we now need to prove the following theorem: Theorem 2 The primitive function ρ(θ) of ρ (θ) (as defined in (A5)) is C 1 (R) and periodic with period at least 2π.
Below the proof of each step: a) As stated in the previous sections, the complex velocity w(z) is an analytic function, and hence the curve describing the trajectory of the vortex point P 1 . This implies that the function ρ(θ) is C 1 (R). Moreover, we can assert that the denominator of ρ (θ) is = 0, or, better, it is easy to show that it is always positive for (ρ, θ) ∈ A. Indeed, the two terms in the denominator in (A5) are always positive (both for sin θ and cos θ positive, negative or null).
c) A sufficient condition to prove the last step is that the function ρ (θ) is an odd function in R. The proof follows directly from (A5) after substituting θ by −θ obtaining: Finally, we apply Theorem 1 to our function ρ (θ) and the theorem is proved. Theorem 2 leads hence to the conclusion that ρ(θ + 2π) = ρ and thus that the trajectory of vortex point P 1 is a closed curve.