1. 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 (Meleshko Reference Meleshko2010) which dates back to von Helmholtz (Reference von Helmholtz1858). 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 (Love Reference Love1894; Hicks Reference Hicks1922; Acheson Reference Acheson2000; Tophøj & Aref Reference Tophøj and Aref2013) to the deformation of the vortex cores and to the effects of viscosity (Shariff & Leonard Reference Shariff and Leonard1992) using numerical (Riley & Stevens Reference Riley and Stevens1993Reference Samson, Wilson, Newman and AndersonReference Wright, Blakestad, Lobb, Phillips and Campbell; Borisov Reference Borisov2014; Cheng & Lim Reference Cheng and Lim2015) as well as experimental methods (Maxworthy Reference Maxworthy1972; Yamada & Matsui Reference Yamada and Matsui1978; Lim Reference Lim1997; Qin, Liu & Xiang Reference Qin, Liu and Xiang2018). The most recent developments concern leapfrogging of vortex bundles (Wacks, Baggaley & Barenghi Reference Wacks, Baggaley and Barenghi2014) and helical waves (Hietala et al. Reference Hietala, Hänninen, Salman and Barenghi2016; Selçuk, Delbende & Rossi Reference Selçuk, Delbende and Rossi2018; Quaranta et al. Reference Quaranta, Brynjell-Rakhola, Leweke and Henningson2019).
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 (White, Anderson & Bagnato Reference White, Anderson and Bagnato2014). 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 (Barenghi & Parker Reference Aioi, Kadokura, Kishimoto and SaitoReference Barenghi and Parker2016) 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 as an ordinary fluid and the total energy is constant, although the incompressible kinetic energy fraction may vary with time being transformed either into acoustic, interaction or potential energy, or vice versa); it is compressible; its vorticity is concentrated to thin hollow vortex lines with fixed width 
$a_0$ and fixed circulation 
$\pm h/m$ where 
$h$ is Planck's constant and 
$m$ is the mass of a boson (while vortices with larger quanta of circulation, 
$\pm 2 h/m, \pm 3 h/m, \ldots ,$ 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 tuneable, 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 (Inouye et al. Reference Inouye, Andrews, Stenger, Miesner, Stamper-Kurn and Ketterle1998)); 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 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 (Serafini et al. Reference Serafini, Galantucci, Iseni, Bienaimé, Bisset, Barenghi, Dalfovo, Lamporesi and Ferrari2017) and vortex points (Seo et al. Reference Seo, Ko, Kim and Shin2017) in real time, including inference of their individual circulations.
Atomic BECs have been employed as a context to study a range of fundamental vortex phenomena, including vortex nucleation from moving obstacles (Frisch, Pomeau & Rica Reference Frisch, Pomeau and Rica1992; Rica Reference Rica2001; Neely et al. Reference Neely, Samson, Bradley, Davis and Anderson2010; Nore, Brachet & Fauve Reference Nore, Brachet and Fauve2012Reference Nowak, Schole, Sexty and GasenzerReference Piazza, Collins and Smerzi; Stagg, Parker & Barenghi Reference Stagg, Parker and Barenghi2014; Kwon, Seo & Shin Reference Kwon, Seo and Shin2015; Musser et al. Reference Musser, Proment, Onorato and Irvine2020) and flow constriction (Valtolina et al. Reference Valtolina2015; Burchianti et al. Reference Burchianti, Scazza, Amico, Valtolina, Seman, Fort, Zaccanti, Inguscio and Roati2018; Xhani et al. Reference Xhani2020), von Kármán vortex streets (Sasaki, Suzuki & Saito Reference Sasaki, Suzuki and Saito2010; Kwon et al. Reference Kwon, Kim, Seo and Shin2016), vortex–antivortex annihilations (Seo et al. Reference Seo, Ko, Kim and Shin2017), vortex line reconnections (Serafini et al. Reference Serafini, Galantucci, Iseni, Bienaimé, Bisset, Barenghi, Dalfovo, Lamporesi and Ferrari2017; Galantucci et al. Reference Galantucci, Baggaley, Parker and Barenghi2019), vortex chaos (Navarro et al. Reference Navarro, Carretero-González, Torres, Kevrekidis, Franzeskakis, Ray, Altuntas and Hall2013), vortex scattering (Barenghi et al. Reference Barenghi, Parker, Proukakis and Adams2005; Caplan et al. Reference Caplan, Talley, Carretero-Gonzales and Kevrekedis2014; Griffin et al. Reference Griffin, Stagg, Proukakis and Barenghi2017), quantum turbulence (Henn et al. Reference Henn, Seman, Roati, Magalhaes and Bagnato2009; Neely et al. Reference Neely, Bradley, Samson, Rooney, Wright, Law, Carretero-Gonzáles, Kevrekidis, Davis and Anderson2013; Kwon et al. Reference Kwon, Moon, Choi, Seo and Shin2014; White et al. Reference White, Anderson and Bagnato2014; Stagg et al. Reference Stagg, Allen, Parker and Barenghi2015; Tsatsos et al. Reference Tsatsos, Tavares, Cidrim, Fritsch, Caracanhas, dos Santos, Barenghi and Bagnato2016; García-Orozco et al. Reference García-Orozco, Madeira, Galantucci, Barenghi and Bagnato2020Reference Gallucci and Proukakis), and self-organisation and clustering of vortices (Billam et al. Reference Billam, Reeves, Anderson and Bradley2014Reference Bisset, Wang, Ticknor, Carretero-Gonzales, Frantzeskakis, Collins and Kevrekidis; Simula, Davis & Helmerson Reference Simula, Davis and Helmerson2014; Gauthier et al. Reference Gauthier, Reeves, Yu, Bradley, Baker, Bell, Rubinsztein- Dunlop, Davis and Neely2019; Johnstone et al. Reference Johnstone, Groszek, Starkey, Billington, Simula and Helmerson2019). With regards to vortex leapfrogging, this has been considered theoretically in idealised unconfined condensates (Ikuta, Sugano & Saito Reference Ikuta, Sugano and Saito2019), including spinor condensates (Kaneda & Saito Reference Kaneda and Saito2014).
Atomic BECs, however, are characterised 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 (Valtolina et al. Reference Valtolina2015; Xhani et al. Reference Xhani2020). The aim of the present work is to provide insight into 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 are 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 the literature neither for classical nor for quantum fluids (Kaneda & Saito (Reference Kaneda and Saito2014) and Ikuta et al. (Reference Ikuta, Sugano and Saito2019) 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 BECs. 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 BECs at temperatures 
$T \ll 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 (Brewczyk, Gajda & Rzazewski Reference Brewczyk, Gajda and Rzazewski2007; Blakie et al. Reference Blakie, Bradley, Davis, Ballagh and Gardiner2008; Proukakis & Jackson Reference Proukakis and Jackson2008; Berloff, Brachet & Proukakis Reference Berloff, Brachet and Proukakis2014). We also remark that, on the one hand, the two-dimensional nature of the system that we consider is an idealisation, but, on the other hand, where atomic BECs 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 § 2, 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 zero-temperature limit. In § 3, 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 § 4, we summarise our findings and illustrate their importance in the future of quantum vortex experiments.
2. Models
2.1. 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 analogue of three-dimensional coaxial vortex rings), each of circulation 
$\pm \kappa$. 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 (Dalfovo et al. Reference Dalfovo, Giorgini, Pitaevskii and Stringari1999Reference Eckel, Lee, Jendrzejewski, Muray, Clark, Lobb, Phillips, Edwards and Campbell), current experimental techniques (Gaunt et al. Reference Gaunt, Schmidutz, Gotlibovych, Smith and Hadzibabic2013) 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 BECs, a vortex is a topological defect of the phase of the governing complex wavefunction (or order parameter), as we shall describe with more detail in § 2.2.1. Therefore, the vortex core is a thin tubular region around the vortex axis which is depleted of atoms: as 
$r \rightarrow 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 
$\xi$ (see § 2.2.1). A similar difference between the classical point vortex model and BECs occurs near a hard boundary: the classical model assumes that the fluid's density is constant up to the boundary; in a BEC a thin boundary region (again of the order of 
$\xi$) 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 BECs provided that vortex–vortex and vortex–boundary distances are larger than the healing length 
$\xi$.
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 (Barenghi et al. Reference Barenghi, Parker, Proukakis and Adams2005). 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 \ll 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 § 3.
Despite these approximations, we believe that the model captures the essential ingredient of our problem: the motion of quantised irrotational vortices in the 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. (Reference Simula, Davis and Helmerson2014). It must also be noticed that Mason, Berloff & Fetter (Reference Mason, Berloff and Fetter2006) 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 
$\xi$ 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.
2.1.1. Equations of motion
Our physical domain under investigation is a two-dimensional infinite strip 
$\mathcal {C}_{\mathbb {R}^{2}}$ defined as 
$\mathcal {C}_{\mathbb {R}^{2}} =\{ (x,y) \in \mathbb {R}^{2}: (x,y) \in ( -\infty , \infty ) \times ( 0 , 2 D ) \}$, which hereafter we will refer to as the channel. We assume the flow to be two-dimensional, 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
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v} = 0. \end{equation}
The velocity field 
$\boldsymbol {v}$ can hence be expressed as the curl of vector field 
$\boldsymbol {\varPsi }$ which, given the two-dimensionality of the flow, has non-vanishing components only in the 
$z$ direction, 
$\boldsymbol {\varPsi } = (0,0,\psi (x,y,t))$. The velocity components have hence the following expressions in terms of the function 
$\psi$ which is often called the streamfunction: 
$v_x=\partial _y \psi$ and 
$v_y=-\partial _x \psi$, where 
$\partial _i$ indicates spatial derivatives in the 
$i$ direction.
The irrotationality of the flow implies that the velocity field can be expressed via a potential function 
$\varphi$, i.e.
\begin{equation} \boldsymbol{v} = \boldsymbol{\nabla} \varphi,\end{equation}
leading to the following relations for the components 
$v_i=\partial _i \varphi$. Equations (2.1) and (2.2) imply that both 
$\varphi$ and 
$\psi$ satisfy the Laplace equation, 
${\rm \Delta} \varphi = {\rm \Delta} \psi = 0$, and the following equalities between their spatial derivatives:
\begin{gather} \partial_x \varphi = \partial_y \psi, \end{gather}
\begin{gather}\partial_y \varphi = -\partial_x \psi. \end{gather}
Equations (2.3) and (2.4) coincide with the well-known Cauchy–Riemann relations for the complex function 
$\varOmega (z) := \varphi + \rm i \psi$, where 
$z=x + \textrm {i} y$. Hence, following basic complex analysis, the function 
$\varOmega (z)$, called the complex potential, is an analytical complex function on the simply connected open domain 
$\mathcal {C}=\{ z\in \mathbb {C}: 0 < \mbox {Im}\, z < 2D \} \subsetneq \mathbb {C}$. As a consequence, 
$\varOmega (z)$ is differentiable and its derivative
\begin{equation} w(z) := \frac{\textrm{d}\varOmega}{\textrm{d} z} = v_x - \textrm{i}v_y \end{equation}
is the so-called complex velocity. In the framework of complex potentials, the impermeable boundary conditions for ideal fluids correspond in our channel 
$\mathcal {C}$ to the following constraint: 
$\mbox {Im}\, \varOmega (z)_{\arrowvert _{z\in \partial \mathcal {C}}} = \alpha (t)$, with 
$\alpha (t)\in \mathbb {R}$ depending only on time 
$t$.
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 
$\varOmega (z)$ describing the velocity field induced by a point vortex in our channel 
$\mathcal {C}$. The essential steps for this derivation are as follows. The necessary ingredients are mainly twofold: (a) the knowledge of the complex potential 
$\varTheta (\zeta )$ describing the flow induced by a point vortex in a simply connected open subset 
$\mathcal {D}$ of the complex plane, with 
$\zeta \in \mathcal {D}\subsetneq \mathbb {C}$; and (b) the construction of a conformal map 
$\zeta =f(z)$ transforming our channel 
$\mathcal {C}$ onto the domain 
$\mathcal {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)\,{\neq}\, 0$ for all 
$z\in \mathcal {C}$. The requirement 
$\mathcal {D}$ not coincide with the entire complex plane 
$\mathbb {C}$, is fundamental in order to exploit the Riemann mapping theorem which ensures the existence of the conformal map 
$f$ mapping 
$\mathcal {C}$ onto 
$\mathcal {D}$. Once 
$\varTheta (\zeta )$ and 
$f(z)$ are determined, the complex potential 
$\varOmega (z)$ for a vortex flow in 
$\mathcal {C}$ is obtained by transforming the potential 
$\varTheta (\zeta )$ via the conformal map 
$f^{-1}(\zeta )$, i.e.
\begin{equation} \varOmega(z)=\varTheta(\,f(z)). \end{equation}
The reasons why the complex function 
$\varOmega (z)$ derived via (2.6) is the sought complex potential are the following. First, 
$\varOmega (z)$ is analytic on 
$\mathcal {C}$ (as it is obtained via the composition of two analytic functions, 
$f$ and 
$\varTheta$), implying that the real and imaginary parts of 
$\varOmega (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 
$\partial \mathcal {C}$ and 
$\partial \mathcal {D}$ under the conformal mapping performed by 
$f$ transposes the boundary conditions enforced by 
$\varTheta (\zeta )$ on 
$\partial \mathcal {D}$ to the boundary 
$\partial \mathcal {C}$ (Lavrentiev & Chabat Reference Lavrentiev and Chabat1972). Finally, via conformal mappings, the flow induced by a vortex of circulation 
$\kappa$ is indeed mapped to a vortex flow with the same circulation (Newton Reference Newton2001).
In the present work, we choose 
$\mathcal {D}$ to coincide with the upper half-complex plane, i.e. 
$\mathcal {D}=\{ \zeta \in \mathbb {C} : \mbox {Im}\, \zeta > 0 \}$. In this domain, the complex potential 
$\varTheta (\zeta )$ describing the flow induced by a vortex placed in 
$\zeta _0\in \mathcal {D}$ is obtained by the method of images, namely
\begin{equation} { \varTheta(\zeta,\zeta_0) = -\text{sgn}(\zeta_0)\frac{\textrm{i}\kappa }{2 {\rm \pi}} \log{\left( \frac{\zeta -\zeta_0}{\zeta - \zeta_0^{*}}\right)}, } \end{equation}
where 
$\text {sgn}(\zeta _0)$ is the sign of the vortex placed in 
$\zeta _0$ (positive for anti-clockwise induced flow, negative for clockwise), 
$\zeta _0^{*}$ is the complex conjugate of 
$\zeta _0$ where a vortex of opposite sign is placed (the image vortex of 
$\zeta _0$) and 
$\kappa$ is the circulation of the flow generated by the vortex. The analytical function 
$f$ transforming conformally the channel 
$\mathcal {C}=\{ z\in \mathbb {C} : 0 < \mbox {Im} \,z < 2D \}$ into 
$\mathcal {D}$ is as follows (see figure 1 for a schematic illustration)
\begin{equation} \zeta=f(z)= \textrm{e}^{{{{\rm \pi} z}/{2D}}}. \end{equation}
The conformal map 
$f$ transforms 
$\partial \mathcal {C}$ into 
$\partial \mathcal {D}$, with 
$f(\{ z\in \mathbb {C} : \mbox {Im} \,z= 0\})=\mathbb {R}^{+}$ and 
$f(\{ z\in \mathbb {C} : \mbox {Im} \,z= 2D\})=\mathbb {R}^{-}$. Employing (2.6), the determination of the complex potential 
$\varOmega (z)$ is straightforward, namely
\begin{equation} {\varOmega(z,z_0)=-\text{sgn}(z_0)\frac{\textrm{i}\kappa}{2{\rm \pi}} \log{\left(\frac{1-\exp\left({{-\dfrac{\rm \pi}{2D}(z-z_{0})}}\right)}{1-\exp\left({{-\dfrac{\rm \pi}{2D}(z-z_{0}^{*})}}\right)}\right)}, \quad z_0=f^{-1}(\zeta_0) } \end{equation}leading to the following complex velocity
\begin{align} w(z,z_0)&=-\text{sgn}(z_0)\frac{\textrm{i}{\rm \pi}}{4}\frac{\kappa}{2{\rm \pi} D}\left \{ \coth \left [ \frac{\rm \pi}{4D}(z-z_0)\right ]-\coth \left [ \frac{\rm \pi}{4D}(z-z_0^{*}) \right ] \right \} \nonumber\\ &= \chi(z,z_0) + \chi(z,z_0^{*}), \end{align}
where 
$\chi (z,z_0)=-\text {sgn}(z_0)({\textrm {i}{\rm \pi} }/{4})({\kappa }/{2{\rm \pi} D}) \coth [({{\rm \pi} }/{4D})(z-z_0)]$ and 
$\text {sgn}(z_0^{*})=-\text {sgn}(z_0)$.
Schematic illustration of the conformal map 
$\zeta =f(z)= \textrm {e}^{{{{\rm \pi} z}/{2D}}}$ transforming 
$\mathcal {C}$ into 
$\mathcal {D}$ and a vortex placed in 
$z_0$ into a vortex in 
$\zeta _0$, 
$\zeta _0=f(z_0)$.
The complex function 
$\chi (z,z_0)$ (and, correspondingly, 
$\chi (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 
$\varOmega (z,z_0)=-\text {sgn}(z_0)\textrm {i}\kappa \log (z-z_0)/(2{\rm \pi} )$) and its infinite images with respect to the walls of the channel, 
$\mbox {Im} \,z=0$ and 
$\mbox {Im} \,z=2D$. The expression (2.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 antivortex in 
$z_0^{*}$ (Greengard Reference Greengard1990).
If the channel is characterised by the presence of 
$N$ vortices, the complex velocity 
$w(z,z_{k_{\{k=1,\ldots ,N\}}})$ generated by the the set of 
$N$ vortices is obtained via the superposition principle, i.e.
\begin{equation} w(z,z_{k_{\{k=1,\ldots,N\}}})=\sum_{k=1}^{N} w(z,z_k)=\sum_{k=1}^{N} \left [ \chi(z,z_k) + \chi(z,z_k^{*}) \right].\end{equation}
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) + \textrm {i} y_j(t)$ occupied by the vortex at time 
$t$ in the channel 
$\mathcal {C}$. Indicating with the superscript ‘
$\,\dot {}\,$’ derivation with respect to time, we define the quantity 
$\dot {z}_j(t):=\dot {x}_j(t) + \textrm {i} \dot {y}_j(t)$, where the real and imaginary parts correspond to the 
$x$ and 
$y$ components of the 
$j$th vortex velocity. As vortices are advected by the local fluid velocity, i.e. 
$\dot {\boldsymbol {x}_j}(t)= \boldsymbol {v}(\boldsymbol {x}_j(t),t)$, the following relation holds:
\begin{equation} \dot{z}_j=w^{*}(z_j,z_{k_{\{k=1,\ldots,N\}}}),\end{equation}
where we have omitted the time dependence of 
$z_j$ and 
$z_k$ to ease notation and the complex conjugation on the right-hand side arises from the definition (2.5) of the complex velocity. In order to determine the complex velocity 
$w(z_j,z_{k_{\{k=1,\ldots ,N\}}})$, we employ (2.11) subtracting the term corresponding to the vortex placed in 
$z_j$, obtaining the following relation:
\begin{align} \dot{z}_j &= w^{*}(z_j,z_{k_{\{k=1,\ldots,N;\ k\!{\neq}j\}}}) + \chi^{*}(z_j,z_j^{*}) \nonumber\\ &= \sum_{k\!{\neq}j} w^{*}(z_j,z_k) + \chi^{*}(z_j,z_j^{*}) \nonumber\\ &= \sum_{k\!{\neq}j} \left [ \chi^{*}(z_j,z_k) + \chi^{*}(z_j,z_k^{*}) \right ] + \chi^{*}(z_j,z_j^{*}), \end{align}
which coincides with the equations of motion of the 
$j$th vortex. The equations of motion for the whole 
$N$-vortex problem are hence a set of 
$2N$ coupled ordinary differential equations.
2.1.2. Hamiltonian formulation and conserved quantities
The motion of 
$N$ vortices located at 
$\boldsymbol {x}_{j_{\{\,j=1,\ldots ,N\}}}$ can be described via a Hamiltonian formalism, where the equations of motion (2.13) may be re-written as follows:
\begin{equation} \left.\begin{aligned} \dot{x}_j = \frac{1}{\text{sgn}(z_j)}\frac{1}{\kappa}\frac{\partial \mathcal{H}}{\partial y_j}\\ \dot{y}_j = - \frac{1}{\text{sgn}(z_j)}\frac{1}{\kappa}\frac{\partial \mathcal{H}}{\partial x_j} \end{aligned}\right\}, \end{equation}
where 
$\mathcal {H}$ is the regular part of the (incompressible) kinetic energy, i.e. the kinetic energy without the vortex self-interaction terms (which we have consistently neglected in the derivation of (2.13) by subtracting 
$\chi ^{*}(z_j,z_j)$ from (2.11)). In an 
$N$-vortex system with circulation 
$\pm \kappa$ the expression of 
$\mathcal {H}$ is as follows (Newton Reference Newton2001)
\begin{equation} \mathcal{H}(t)=\frac{\kappa}{2}\sum_{k=1}^{N} \text{sgn}(z_k) \psi_{\mathcal{C}}(\boldsymbol{x}_{k},t) =\frac{\kappa}{2}\sum_{k=1}^{N} \text{sgn}(z_k) \mbox{Im} \left [\varOmega_{\mathcal{C}} (z_k,t) \right], \end{equation}
where with 
$\psi _{\mathcal {C}}$ and 
$\varOmega _{\mathcal {C}}$ we indicate respectively the overall streamfunction and the overall complex potential arising from the contribution of all vortices, i.e.
\begin{align} \varOmega_{\mathcal{C}} (z_k,t) = \sum_{l{\neq} k}\varOmega(z_k(t),z_l(t)) + \text{sgn}(z_k)\frac{\textrm{i}\kappa}{2{\rm \pi}}\log \left [ 1-\exp\left({{-\frac{\rm \pi}{2D}(z_k(t)-z_{k}^{*}(t))}}\right)\right], \end{align}
where 
$\varOmega (z_k,z_l)$ is given by (2.9) and the last term is obtained by removing the singular part from 
$\varOmega (z_k,z_k)$. In an 
$N$-vortex problem, 
$\mathcal {H}(t)$ is a conserved quantity, given the absence of dissipative phenomena in an ideal fluid. In the circumstance where 
$\mathcal {H}$ is invariant with respect to translations in the 
$x$ or 
$y$ directions or with respect to rotations, other constants of motions exist, namely
\begin{gather} I_x(t) = \kappa \sum_{k=1}^{N} \text{sgn}(z_k) y_k(t), \end{gather}
\begin{gather}I_y(t) = -\kappa \sum_{k=1}^{N} \text{sgn}(z_k) x_k(t), \end{gather}
\begin{gather}L(t) = -\frac{\kappa}{2}\sum_{k=1}^{N} \text{sgn}(z_k) [x_k^{2}(t) + y_k^{2}(t)]. \end{gather}
The 
$x$ and 
$y$ components of the linear impulse 
$\boldsymbol {I}=(I_x, I_y)$ are conserved if the domain is invariant with respect to translations in the 
$x$ and 
$y$ directions, respectively, while the angular momentum 
$L$ is conserved if the geometry of the system is invariant under rotations. Thus, in our two-dimensional channel, 
$I_x$ is always conserved, while in general 
$I_y$ and 
$L$ depend upon time.
2.2. 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 
$\varPsi = \varPsi (\boldsymbol {x},t)$ of the system, and reads as follows:
\begin{equation} \textrm{i}\hbar \dot{\varPsi} = -\frac{\hbar^{2}}{2m}{\rm \Delta} \varPsi + V \varPsi + g |\varPsi|^{2} \varPsi, \end{equation}
where the dot is the time derivative, 
$\hbar =h/(2{\rm \pi} )$ is the reduced Planck's constant, 
$m$ is the boson mass, 
$V=V(\boldsymbol {x},t)$ is an externally applied potential and 
$g=4{\rm \pi} \hbar ^{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 
$\varPsi$ can be written in terms of its amplitude and its phase as
\begin{equation} \varPsi = \sqrt{n} \,\textrm{e}^{\textrm{i}\theta}, \end{equation}
where 
$n=n(\boldsymbol {x},t)=|\varPsi |^{2}$ is the particle number density (number of bosons per unit volume) and 
$\theta =\theta (\boldsymbol {x},t)$ is the phase. Without loss of generality, the order parameter 
$\varPsi$ can be written as 
$\varPsi (\boldsymbol {x},t)=\textrm {e}^{\textrm {i} \mu t/\hbar } \varPhi (\boldsymbol {x},t)$, where 
$\mu$ is called the chemical potential and 
$\varPhi (\boldsymbol {x},t)$ obeys
\begin{equation} {\textrm{i}\hbar \dot{\varPhi} = -\frac{\hbar^{2}}{2m}{\rm \Delta} \varPhi + V \varPhi + g |\varPhi|^{2} \varPhi -\mu \varPhi.} \end{equation}2.2.1. Quantum vortices
In the context of BECs described by the GP equation, quantum vortices are topological defects of the phase 
$\theta$ of the order parameter, at which 
$\varPsi =0$ (hence 
$\theta$ is undefined) and around which 
$\theta$ wraps by 
$2q{\rm \pi}$ with 
$q\in \mathbb {Z}\setminus \{ 0 \}$. In three dimensions, vortices take the form of one-dimensional curves which may form a vortex tangle, as observed both in BECs (White et al. Reference White, Barenghi, Proukakis, Youd and Wacks2010) and superfluid helium (Vinen Reference Vinen1957). In two dimensions, vortices coincide with vortex points which have been observed extensively in oblate (pancake-like) BECs (Matthews et al. Reference Matthews, Anderson, Haljan, Hall, Wieman and Cornell1999). For the purpose of the present work, we will restrict our discussion to two-dimensional systems.
The velocity field 
$\boldsymbol {v}(\boldsymbol {x},t)$ associated with a BEC whose dynamics is described by the order parameter 
$\varPsi$, is defined from the phase 
$\theta$ via the relation
\begin{equation} \boldsymbol{v}(\boldsymbol{x},t)=\frac{\hbar}{m}\boldsymbol{\nabla} \theta. \end{equation}
Employing the definition (2.23) of the velocity and the 
$2q{\rm \pi}$ phase wrapping existing around a vortex, it is straightforward to verify that the circulation 
$\varGamma$ of the velocity field on any closed curve 
$\gamma$ enclosing a vortex point is quantised in terms of the quantum of circulation 
$\kappa =h/m$, i.e.
\begin{equation} \varGamma = \oint_\gamma \boldsymbol{v}\cdot \boldsymbol{dl}=q \kappa,\quad q\in\mathbb{Z}\setminus \{ 0 \}. \end{equation}
Choosing 
$\gamma$ 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_\phi =q\kappa /(2{\rm \pi} 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 
$\xi =\hbar /\sqrt {mgn}$. As we will very briefly illustrate in the next section, quantum fluids are indeed compressible fluids.
2.2.2. Fluid dynamical equations for a BEC
The GP equation (2.20) may be rewritten via the Madelung transformation consisting in expressing 
$\varPsi$ in polar form (2.21) and separating the real and imaginary parts of (2.20). This procedure leads to the following equations:
\begin{gather} \dot{n} + \boldsymbol{\nabla}\boldsymbol{\cdot} \left ( n\boldsymbol{v}\right ) = 0, \end{gather}
\begin{gather}m n \left [ \dot{\boldsymbol{v}} + \left ( \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} \right ) \boldsymbol{v} \right ] = -\boldsymbol{\nabla} \left ( p + p'\right ) -n\boldsymbol{\nabla} V, \end{gather}
where 
$p$ and 
$p'$ are respectively pressure and quantum pressure
\begin{equation} p = \frac{g n^{2}}{2}, \end{equation}
\begin{equation}p'= -\frac{\hbar^{2}}{4m} n {\rm \Delta} \left ( \ln{(n)} \right ). \end{equation}
Equation (2.25) coincides formally with the continuity equation of a classical fluid, while (2.26), 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 
$\ell$ much larger than the healing length 
$\xi$ (which is the typical length scale for density variations, associated e.g. with the presence of vortices or boundaries) 
$p'/p \ll 1$, implying that in this limit the BEC can indeed be considered as a barotropic, compressible classical inviscid fluid. Hence, at length scales 
$\ell \gg \xi$, 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 
$\ell \sim \xi$, the physics may be significantly different. For instance, if the relative distance between quantum vortices of opposite sign is of the order of 
$\xi$, 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 
$\omega _\tau$ as the distance 
$\ell$ tends to zero (in the classical model, the frequency diverges, 
$\omega _\tau \sim 1/\ell ^{2}$).
3. Results
3.1. 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 
$\kappa$, half with negative 
$-\kappa$. In figure 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,\ldots ,4$ and the vortex–antivortex pairs are symmetrically positioned with respect to the channel mid-axis 
$y=D$, namely 
$y_j(0)=D \pm R$ for the first pair 
$j=(1,2)$ and 
$y_j(0)=D \pm r$ for the second pair 
$j=(3,4)$, with the conditions 
$R/D < 1$ and 
$r/R < 1$. We have verified that the initial vortex symmetry with respect to the channel mid-axis 
$y=D$ is conserved along the motion of the vortices. This symmetry plays an important role in the dynamics of vortices as it implies that, in addition to the horizontal impulse 
$I_x$ ((2.17), always conserved as a result of the channel geometry), also the vertical impulse 
$I_y$ (2.18) and the angular momentum 
$L$ (2.19) are constants of motion. In more detail, in our channel, with the mentioned symmetry with respect to 
$y=D$, we have 
$I_x = 2\kappa (R + r)$, 
$I_y=0$ and 
$L=-2D\kappa (R + r)$.
Initial vortex configuration for the classical point vortex numerical simulations: filled (open) circles correspond to vortices with positive (negative) circulation. Numerical labels close to vortices indicate the vortex numeration employed.
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 (2.13) for the four vortices, 
$j=1,\ldots ,4$, varying 
$r/R$ and 
$R/D$. In particular, we choose 
$r/R=n/10$ and 
$R/D=m/10$, with 
$m, n =1, \ldots , 9$. The time-advancement scheme employed in the numerical simulations is a second-order Adams–Bashforth method with a time step 
${\rm \Delta} t=T/1000$ where 
$T=2{\rm \pi} ^{2}\delta ^{2}/\kappa$ is the rotation period of a pair of vortices of the same polarity placed at distance 
$\delta$. In our numerical simulations 
$\delta$ is set to 
$10^{-3}D$.
For classical unbounded fluids, since the study performed by Love over a century ago (Love Reference Love1894), it is well known that vortices undergo leapfrogging motion only if 
$r/R$ is larger than a critical value 
$\alpha _c = 3 - 2\sqrt {2} \approx 0.172$. If 
$r/R < \alpha _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 (Reference Acheson2000) extended numerically the study performed by Love and established that leapfrogging motion is unstable when 
$\alpha _c < r/R < \alpha _c'$, with 
$\alpha _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 the distinction between leapfrogging and non-leapfrogging, which is already known, we also observe image-driven leapfrogging and periodic orbits. The phase diagram of the system resulting from the numerical simulations is illustrated in figure 3.
Phase diagram of the classical motion of two vortex–antivortex pairs in a two-dimensional plane channel. All symbols refer to performed numerical simulations. Black circles indicate no leapfrogging motion; red filled (open) squares stand for stable (unstable) forward, standard leapfrogging; blue filled (open) diamonds correspond to stable (unstable) image-driven leapfrogging; green stars stand for periodic orbits. The dashed green lines indicate the analytical solution for periodic orbits (see § 3.1.1 and appendix A). Dark and light green colours indicate stable and unstable periodic orbits, as discussed in § 3.3. The dashed violet curve is the numerically computed 
$\alpha _c$ dependence on 
$(R/D)$: for each value of 
$(R/D)$ considered, we increase 
$r/D$ by 0.001 until we observe the onset of leapfrogging motion, identifying the corresponding critical value 
$(r/D)_c$. The value of 
$\alpha _c$ is then determined as 
$\alpha _c=(r/D)_c/(R/D)$ and indicated with small violet circles.
For values of 
$R/D \le 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 \le 1/2$, in fact, as we increase 
$r/R$, we first observe non-leapfrogging motion (in black in figure 3), defined as the dynamics characterised by 
$\dot {y}_j(t)=0$ for all 
$j$ at late times; then we notice unstable leapfrogging motion (open red squares, we refer the reader to § 3.3 for the definition of stable/unstable trajectories), and finally stable leapfrogging (filled red squares). These dynamical regimes therefore coincide with the scenario outlined by Acheson (Reference Acheson2000), the only significant and important difference being the dependence of 
$\alpha _c$ on 
$R/D$: for small values of 
$R/D$, 
$\alpha _c$ is very close to the constant value 
$0.172$ for vortex leapfrogging in unbounded fluids (e.g. for 
$R/D=0.1$, 
$\alpha _c=0.173$), increasing for increasing values of 
$R/D$ (e.g. 
$\alpha _c=0.216$ for 
$R/D=0.5$). This dependence of 
$\alpha _c$ on 
$R/D$ stems from the interaction of the outer vortices 1 and 2 in figure 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 monotonic behaviour of 
$\alpha _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 image-driven leapfrogging, indicated by blue diamonds in figure 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–antivortex pair on its own. As a consequence, we observe two distinct leapfrogging motions, each involving two virtual vortex–antivortex 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 image-driven leapfrogging transition occurs via an intermediate regime in which vortices follow periodic orbits, indicated by dark and light green stars in figure 3 (dark and light green correspond to stable and unstable periodic orbits respectively, as described in § 3.3). As shown in detail in the next section and in the analytical derivation presented in appendix A, periodic orbits are observed when 
$R + r =D$, corresponding to the green dashed line in figure 3. For large values of 
$R/D$ (
$R/D\gtrsim 3/4$), the system crosses directly the no-leapfrogging to image-driven 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 figure 4. Note that in the three-dimensional coaxial vortex ring analogue, vortices of the same colour correspond to cross-sections of the same vortex ring.
Examples of dynamical regimes and trajectories for classical 4-vortex motion in a two-dimensional channel. Filled (open) symbols indicate positive (negative) vortices: (a) 
$R/D=5/10$, 
$r/R=1/10$, no leapfrogging (vortices moving to the right); (b) 
$R/D=4/10$, 
$r/R=4/10$, forward (standard) leapfrogging (vortices moving to the right); (c) 
$R/D=8/10$, 
$r/R=6/10$, image-driven leapfrogging (vortices moving to the left); (d) 
$R/D=6/10$, 
$r/R=67/100$, periodic orbits, showing more than 30 overlapping closed orbits.
3.1.1. 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 appendix A.
With reference to figure 2, we consider the pair of vortices 
$P_1 =(x_0(t),D-R(t))$, with negative circulation 
$-\kappa$, and 
$P_2= (x_0(t),D+R(t))$, with positive circulation 
$\kappa$, and the pair of vortices 
$P_3 =(x_1(t),D-r(t))$, with negative circulation 
$-\kappa$, and 
$P_4= (x_1(t),D+r(t))$, with positive circulation 
$\kappa$, where 
$t$ is time. In the complex domain, omitting the time dependence to ease notation, these vortices are located in 
$z_1 =x_0 +\textrm {i} (D-R)$ for 
$P_1$, 
$z_2 =x_0 +\textrm {i} (D+R)$ for 
$P_2$, 
$z_3 =x_1 +\textrm {i} (D-r)$ for 
$P_3$ and 
$z_4 =x_1 +\textrm {i} (D+r)$ for 
$P_4$, and they generate the following complex velocity in the point 
$z$, as given by (2.11):
\begin{equation} w(z)=w(z,z_1)+w(z,z_2)+w(z,z_3)+w(z,z_4). \end{equation}
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}+\textrm {i}(D-{(r(t)+R(t))}/{2})$ and the complex velocity generated by vortices in 
$z_M$ which we indicate with 
$w(z_M)$
\begin{equation} w(z_M)= \frac{\dfrac{\textrm{i} \kappa}{2 D} \left(-1+\exp\left({\dfrac{2 \textrm{i} {\rm \pi}(r+R)}{D}}\right)\right) \exp\left({\dfrac{{\rm \pi} ({x}_0+{x}_1)}{2 D}}\right)}{\begin{array}{l} \left(\exp\left({\dfrac{{\rm \pi} (4 \textrm{i} r+4 \textrm{i} R+{x}_0+{x}_1)}{2 D}}\right)- \exp\left({\dfrac{{\rm \pi} (2 {x}_0+\textrm{i} (r+3 R))}{2 D}}\right) \right.\\ \quad \left.-\exp\left({\dfrac{{\rm \pi} (2 {x}_1+\textrm{i} (3 r+R))}{2 D}}\right)+ \exp\left({\dfrac{{\rm \pi} ({x}_0+{x}_1)}{2 D}}\right)\right) \end{array}}. \end{equation}
If we look for the conditions such that the velocity 
$w(z_M)$ of the midpoint 
$M$ is zero, we have
\begin{align} w(z_M)= 0 &\Longleftrightarrow \exp\left({\frac{2 \textrm{i} {\rm \pi}(r+R)}{D}}\right)-1=0 \nonumber\\ &\Longleftrightarrow \frac{2{\rm \pi} (r+R)}{D} = 2k {\rm \pi},\quad k\in \mathbb{Z}. \end{align}
Note that the same result equation (3.3) 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 (3.3) are 
$k \in \mathbb {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.
\begin{equation} r(t)+R(t)=D. \end{equation}
This is the most interesting result: it states that when the four vortices satisfy the condition (3.4) then the midpoints 
$M$ and 
$N$ are at rest; the two pairs of vortices 
$(P_1,\ P_3)$ and 
$(P_2,\ P_4 )$ hence move symmetrically with respect to their corresponding midpoints, i.e. 
$\dot {x}_0(t)=-\dot {x}_1(t)$ and 
$\dot {R}(t)=-\dot {r}(t)$. The last equality is fundamental as it expresses that if condition (3.4) 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=\textrm {i} ({D}/{2})$ and 
$z_{N}= \textrm {i} ({3D}/{2})$.
The last step to demonstrating 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 
$\dot {R}(t)=-\dot {r}(t)$ (for instance, 
$\dot {R}(t)=\dot {r}(t)=0$) could be possible, not leading to periodic orbits. We tackle this issue in appendix A, to ease the readability of the manuscript.
3.2. 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 § 3.1.
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 GP formulation of the problem we have an extra length scale – the healing length 
$\xi$ – 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\xi$ and 
$D_2=20\xi$. In order to model the channel confinement, we use the following potential 
$V$:
\begin{equation} V = V(\,y) = \begin{cases} 0 & \text{if } 0 < y < 2D \\ 10\mu & \text{if } y \le 0 \text{ or } y \ge 2D, \end{cases}\end{equation}
corresponding to a channel of half-width 
$D$, where the density 
$|\varPhi |^{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 
$\varPhi$, the dimensionless GP equation
\begin{equation} \textrm{i}\dot{\varPhi} = -\frac{1}{2}{\rm \Delta} \varPhi + \frac{V}{\mu}\varPhi + |\varPhi|^{2} \varPhi - \varPhi. \end{equation}
Equation (3.6) is obtained from (2.22) after introducing characteristic units of length, time and energy: 
$\xi =\hbar /\sqrt {m\mu }$ (the healing length), 
$\tau =\xi /c$ (where 
$c=\sqrt {\mu /m}$ is the speed of sound) and 
$\mu$ (the chemical potential) respectively, and normalising the order parameter with respect to the unperturbed homogeneous solution 
$\varPhi _0=\sqrt {\mu /g}$ of (2.22). In these units the healing length and the bulk density in the channel are unity.
The numerical integration of (3.6) is performed employing a fourth-order Runge–Kutta time advancement scheme and second-order finite differences to approximate spatial derivative operators. Time step 
${\rm \Delta} t/\tau$ is set to 
$1.5\times 10^{-2}$ and spatial discretisation 
${\rm \Delta} x/\xi = {\rm \Delta} y/\xi$ is chosen to be equal to 
$0.25$. In the set of simulations where 
$D=D_1=40\xi$, 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\xi \le x \le 800\xi$ and 
$-10\xi \le y \le 90\xi$. On the other hand, when 
$D=D_2=20\xi$, 
$N_x=3200$ and 
$N_y=240$, respectively, leading to the computational box 
$-400\xi \le x \le 400\xi$ and 
$-10\xi \le y \le 50\xi$. We employ periodic boundary conditions in the 
$x$ direction and reflective in the 
$y$ direction where the wave function 
$\varPhi \rightarrow 0$ well before reaching the boundary of the numerical domain, given the buffer layer of 
$10\xi$.
The initial imprinting of vortices is made by enforcing a uniform 
$2{\rm \pi}$ 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 (3.6) for 
$t\in \mathbb {R}$. In figure 5 we report the density 
$|\varPhi |^{2}(x,y)$ (a) and the phase 
$\theta (x,y)$ (b) of the initial condition employed for 
$R/D=0.6$ and 
$r/R=0.3$ and 
$D=D_1=40\xi$. It can be easily observed that the density 
$|\varPhi |^{2}$ rapidly drops to zero at the vortex positions and outside the channel. Correspondingly, the four 
$2{\rm \pi}$ phase wrappings can be distinguished in figure 5(b).
Initial condition for numerical simulation of leapfrogging of quantum vortices in a two-dimensional channel for 
$R/D=0.6$ and 
$r/R=0.3$ and 
$D=D_1=40\xi$. (a) The density of the BEC 
$|\varPhi (x,y)|^{2}$ (presented as a ratio of the bulk density 
$|\varPhi _0|^{2}$) is displayed: it is unity (yellow) in the bulk of the channel and vanishes (blue) in the vortex cores and at the channel's boundaries; (b) the phase 
$\theta (x,y)$ of the BEC is illustrated in the range 
$[-{\rm \pi} ,{\rm \pi} )$. We underline that these figures only show a tenth of the numerical domain which spans from 
$-800\xi$ to 
$800\xi$ in the 
$x$ direction.
To verify the existence in a BEC of all distinct regimes observed in the classical point vortex model (§ 3.1), we perform numerical simulations of quantum vortex leapfrogging along the vertical line 
$R/D=0.6$ of the phase diagram reported in figure 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. To reconstruct the vortex configuration, we print to file the non-dimensional wavefunction every 100 time steps; the position of each vortex is precisely determined by monitoring the density as well as the circulation.
The results are schematically outlined in figure 6, where classical vortex dynamics (left) is compared to quantum vortex dynamics at 
$D=D_1=40\xi$ (middle) and 
$D=D_2=20\xi$ (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: the two cuts in the phase diagram coincide. When 
$D=D_2$ we observe two differences: first, at 
$(R/D,r/R)=(0.6, 0.1)$ the internal vortex–antivortex pair annihilates as their initial distance is only 
$2.4\xi$; second, the no leapfrogging and the unstable leapfrogging regions are shifted towards larger values of 
$r/R$. For 
$r/R>0.5$ the phase diagram at 
$R/D=0.6$ is identical to the classical diagram. These differences at the smaller value of channel size 
$D$ are expected, as the healing length scale starts playing a role: only if 
$D/\xi$ is sufficiently large we can expect classical and quantum dynamics to be the same.
Cuts in the dynamical regimes phase diagram corresponding to 
$R/D=0.6$ for the classical point vortex model (left), the GP model with 
$D=D_1=40\xi$ (middle) and 
$D=D_2=20\xi$ (right). Symbols as in figure 3 except for the newly introduced up-pointing orange triangle corresponding to the annihilation of the inner vortex–antivortex pair.
The matching of the observed dynamical regimes when comparing classical and quantum leapfrogging in a two-dimensional channel if 
$D\ge 40\xi$ is confirmed in figure 7, which shows the trajectories of quantum vortices for 
$(R/D, r/R)$ pairs selected as for the classical trajectories illustrated in figure 4. We in fact observe that the same regimes (no leapfrogging/leapfrogging/image-driven leapfrogging/periodic orbits) occur in the quantum system for the corresponding parameters 
$(R/D, r/R)$ selected for the classical case.It is, however, worth noting some minor differences between the quantum vortex trajectories and their classical counterparts reported in figure 4. Since the initial condition is not stationary with respect to any frame of reference, when we start integrating in time (3.6) for 
$t\in \mathbb {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 Frisch et al. (Reference Frisch, Pomeau and Rica1992)), is visible in figure 7(b,c,d) when compared with figure 4. In particular, this horizontal shift affects the periodic orbits reported in figure 7(d) whose centre is slightly shifted towards positive 
$x$ values. The orbits illustrated in figure 7(d) are indeed periodic, as, after the initial sound-induced translation towards the right, the vortex trajectories do follow the same closed orbits (more than five closed orbits are shown in figure 7). For the sake of completeness, in appendix B we show the comparison between classical and quantum vortex trajectories for 
$(R/D, r/R) = (6/10, 67/100)$ (d of figures 7 and 4) and 
$(R/D, r/R) = (72/100, 39/100)$. For both choice of parameters, as 
$R+r=D$, classical vortices perform periodic orbits, while in the quantum case, for 
$(R/D, r/R) = (72/100, 39/100)$ vortices perform a very slow image-driven leapfrogging. We can hence conclude that while all regimes identified in the classical point vortex model occur in quantum fluids, not necessarily the same regime takes place at given parameters 
$(R/D, r/R)$.
Dynamical regimes observed in a quantum 4-vortex configuration in a two-dimensional channel with 
$D=D_1=40\xi$. Filled (open) symbols indicate positive (negative) vortices: (a) 
$R/D=5/10$, 
$r/R=1/10$, no leapfrogging motion is observed; (b) 
$R/D=4/10$, 
$r/R=4/10$, forward (standard) leapfrogging; (c) 
$R/D=8/10$, 
$r/R=6/10$, image-driven leapfrogging; (d) 
$R/D=6/10$, 
$r/R=67/100$, periodic orbits, showing more than five completed periods.
In addition, we observe the number of periods observed in the 
$x$ range 
$[-2D, 2D]$ is different from the classical counterpart. This difference in the average 
$x$ velocity of the vortices cannot be ascribed to the well-known additional average velocity 
$u_{GP}$ in the 
$x$ direction possessed by quantum vortices (with respect to the classical system) when periodic boundary conditions are used in GP numerical simulations (Griffin et al. Reference Griffin, Shukla, Brachet and Nazarenko2004). In fact, the periodic boundary conditions employed in the 
$x$ direction imply that the constant superposed velocity 
$u_{GP}$ is given by 
$u_{GP}=\kappa (R+r)/(2 D L_x)$, where 
$L_x$ is the length of the numerical domain in the 
$x$ direction. As the 
$r < R < D$ and in figure 7 
$L_x/D=40$, the magnitude of 
$u_{GP}$ is too small to account for the distinct average velocity of vortices in the 
$x$ direction observed while comparing classical and quantum trajectories. These dissimilarities between classical and quantum trajectories are possibly due to the compressible nature of a quantum Bose gas and hence to the interplay between vortices and sound. An example of vortex–sound interaction in BECs reported in the literature (Parker et al. Reference Parker, Proukakis, Barenghi and Adams2004) is the transformation of incompressible kinetic energy into compressible kinetic energy (sound) when vortex–vortex pairs of the same polarity rotate in a dimple trap: the vortex pair radiates sound similarly to charged particles radiating electromagnetic waves. The role played by this effective dissipation of kinetic energy into sound will be assessed in a future study.
3.3. Stability of periodic orbits
In this section we address the numerical stability of periodic orbits for the leapfrogging of vortices occurring in both classical and quantum fluids. The aim is to provide an optimal range of parameters 
$(R/D, r/R)$ for the experimental observation of such novel leapfrogging regime in confined BECs. In order to probe the stability of the orbits in the classical point vortex model, we follow the green line 
$R+r=D$ in figure 3 and for each simulation indicated with a green star (corresponding to a step increase of 0.02 in 
$R/D$) we perform an additional simulation where the top vortex in figure 2 is initially translated by a small amount 
$\epsilon$ in the 
$x$ direction, i.e. 
$x_2(0)=\epsilon$ (Acheson Reference Acheson2000). We choose 
$\epsilon =\eta R$, with 
$\eta =0.005$. We then numerically determine the time 
$t^{\ast }$ where for the first time 
$|y_2(t^{\ast }) - y_N| > 2\mathcal {A}$, where 
$\mathcal {A}=(R-r)/2$ is the amplitude of the unperturbed (
$\epsilon =0$) periodic orbits of vortices 2 and 4 (see figure 2) and 
$y_N=3D/2$ is the vertical coordinate of the centre 
$N$ of such periodic orbits (cf. § 3.1.1). Conceptually, the time 
$t^{\ast }$ indicates the amount of time which is necessary for vortex 2 to reach a vertical distance to the centre 
$N$ larger than twice the amplitude 
$\mathcal {A}$ of the corresponding unperturbed periodic orbit. As the period 
$T$ of periodic orbits depends upon 
$R/D$, we successively define the ratio 
$\mathcal {P}=t^{\ast }/T$ for each value of 
$R/D$, and consider a periodic orbit unstable if 
$\mathcal {P}<5$. With this criterion, we conclude that periodic orbits are stable for 
$R/D < 7/10$. This is due to the fact that when 
$R/D \rightarrow 1/2$, 
$r$ tends to the value of 
$R$, i.e. pairs of vortices of the same polarity are closer. As a consequence, the predominant dynamics is the relatively fast rotation of same signed vortices which is less affected by the initial displacement 
$\epsilon$. In the other limit of large 
$R/D$, the overall velocity of vortices is slower and hence more affected by the non-symmetrical initial vortex configuration. In figure 3 we indicate unstable (stable) periodic orbits with light (dark) green stars. This procedure based on the threshold value of 
$\mathcal {P}=5$ has been employed throughout the paper for distinguishing stable from unstable vortex trajectories.
In figure 8 we report vortex trajectories with 
$x_2(0)=\epsilon =\eta R$ for 
$R/D=6/10, 66/100, 7/10$ (classical point vortex model) and 
$R/D=7/10$ (GP model). Relatively to the classical point vortex model, we observe that as 
$R/D$ increases the number of periods 
$\mathcal {P}$ completed by the vortices before crossing the line 
$y = y_N - 2\mathcal {A}$ decreases, eventually falling below 5. We also observe that for 
$R/D=7/10$, quantum trajectories follow quite closely the corresponding classical trajectories. Together with the similarity between classical and quantum trajectories already described in the previous paragraph 3.2, this suggests that the conclusions concerning the stability of periodic orbits in the classical case may be extended to the quantum case.
Vortex trajectories with the initial position of vortex 2 (cf. figure 2) being translated by the small quantity 
$\epsilon =\eta R$, with 
$\eta =0.005$, in the positive 
$x$ direction. Symbols as in figures 4 and 7. Green dashed lines correspond to 
$y = y_N + 2\mathcal {A}$ and 
$y = y_N - 2\mathcal {A}$: (a) 
$R/D=6/10$, 
$r/R=67/100$, classical point vortex model; (b) 
$R/D=66/100$, 
$r/R=52/100$, classical point vortex model; (c) 
$R/D=7/10$, 
$r/R=43/100$, classical point vortex model; (d) 
$R/D=7/10$, 
$r/R=43/100$, quantum GP model. All trajectories correspond to the same interval in physical time.
In the last paragraph of this section we address the possibility of observing periodic orbits in experimental accessible BECs. We believe that in order to achieve this aim, three criteria have to be satisfied: (a) periodic orbits have to satisfy the stability requirement determined previously, i.e. 
$R/D < 7/10$; (b) vortices have to be distinguishable experimentally and therefore their distance has to greater than at least twice the size of a vortex core, i.e. 
$R-r > d_{core}=10\xi$ (as a vortex core is of the order of five healing lengths); (c) the period of the periodic orbits has to be larger than the time interval between consecutive images of the BEC. As we are interested in periodic orbits, condition (b) can be written as follows 
${R}/{D} > \frac {1}{2} + {d_{core}}/{D}$, which becomes 
${R}/{D} > 0.625$ if we consider a width 
$D=40\xi$, easily achievable in current experimental facilities. As for condition (c), we first approximate for simplicity the period 
$T$ of periodic orbits as the period of two co-rotating vortices in an unbounded BEC, i.e. 
$T={2{\rm \pi} ^{2}(R-r)^{2}}/{\kappa }$. We then take into account the experimental setting employed in previous work (Serafini et al. Reference Serafini, Galantucci, Iseni, Bienaimé, Bisset, Barenghi, Dalfovo, Lamporesi and Ferrari2017) where the imaging time interval 
${\rm \Delta} t_{img} \sim 10\ \text {ms}$ and the atomic species condensed is 
${}^{23}\text {Na}$. With these parameters and requiring 
$T > 2 {\rm \Delta} t_{img}$ we have that 
$R/D > 0.64$. Thus, combining conditions (a)–(c) we have that the optimal range of parameters for observing experimentally periodic orbits is 
$0.64 < R/D < 0.7$.
4. 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: image-driven leapfrogging and periodic motion. Using the GP equation to model an atomic Bose–Einstein condensate (a compressible quantum fluid) confined within a channel, we have verified that all four regimes also exist for quantum vortices. In large channels, the boundaries between these regimes are the same for classical and quantum vortices although some differences arise while comparing the trajectories of individual vortices. Further 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 occurring in confined geometries will be particularly important for the interpretation and planning of ongoing and future experiments with atomic BECs, 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.
Funding
L.G., N.G.P. and C.F.B. acknowledge the support of the Engineering and Physical Sciences Research Council (grant no. EP/R005192/1). M.S. acknowledges the support of MIUR-Italy through the project PRIN ‘Multiscale phenomena in Continuum Mechanics: singular limits, off-equilibrium and transitions’ (grant no. PRIN2017 2017YBKNCE) and the support of the Dipartimento di Ingegneria of Universitá di Palermo under grant no. FFR_D26_PREMIO_SINGOLO_RICERCA_2020.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of periodic orbits
In order to show the existence of periodic orbits, we have to prove that, if condition (3.4) is satisfied, the trajectories of the vortex points are closed curves with vortices rotating around the two midpoints 
$M$ and 
$N$ defined in § 3.1.1.
For the sake of simplicity, and with reference to § 3.1.1, we prove the closedness of the trajectory only for the vortex point 
$P_1$, as the proof for the other vortex points is an iterative procedure.
We consider (2.12) for the vortex point 
$z_1$ with the complex velocity given by the expression (3.1) evaluated on the vortex point 
$z_1$. Since the middle point 
$z_M$ is at rest for 
$r+R=D$, we rewrite the dynamic equation of 
$z_1$, namely 
$\dot {z_1}= w^{*}(z_1)$, in the polar coordinate system 
$(\rho ,\theta )$ centred on 
$z_M$. The middle point 
$z_M$, under the condition 
$r+R=D$, becomes 
$z_M={(x_0+x_1)}/{2}+\textrm {i}({D}/{2})$, which requires the condition 
$\rho < {D}/{2}$ to ensure that vortices 
$P_1$ and 
$P_3$ are in 
$\{ z\in \mathbb {C} : 0 < \mbox {Im} \,z < D \}$.
Thus, in the new reference system the vortex points correspond to
\begin{equation} \left.\begin{aligned} & z_1= z_M-\rho \cos(\theta)-\textrm{i} \rho \sin(\theta), \\ & z_2= z_M+iD-\rho \cos(\theta)+\textrm{i} \rho \sin(\theta), \\ & z_3= z_M+\rho \cos(\theta)+\textrm{i} \rho \sin(\theta), \\ & z_4= z_M+iD+\rho \cos(\theta)-\textrm{i} \rho \sin(\theta), \end{aligned}\right\} \end{equation}
where 
$z_M$ is now the origin of the new frame of reference, which can be set 
$z_M=0 + 0 i$. Note that the condition (3.4) is automatically satisfied by construction; indeed, 
${(z_2-z_1)}/{2}=\textrm {i}({D}/{2}+\rho \sin (\theta ))$ and 
${(z_4-z_3)}/{2}=\textrm {i}({D}/{2}-\rho \sin (\theta ))$, implying 
$R\equiv {D}/{2}+\rho \sin (\theta )$ and 
$r\equiv {D}/{2}-\rho \sin (\theta )$ and, hence, condition (3.4). We now substitute the coordinates (A1) into the equation
\begin{equation} \dot{z}_1=w^{*}(z_1,z_{k_{\{k=1,\ldots,4\}}}), \end{equation}
according to (2.12), and change the vectorial basis from 
$(\hat {\boldsymbol {x}},\ \hat {\boldsymbol {y}})$ to 
$(\hat {\boldsymbol {u}}_\rho ,\ \hat {\boldsymbol {u}}_\theta )$ by means of the following rotation:
\begin{equation} \hat{\boldsymbol{u}}_\rho =\cos(\theta) \hat{\boldsymbol{x}} + \sin(\theta) \hat{\boldsymbol{y}},\quad \hat{\boldsymbol{u}}_\theta= -\sin(\theta) \hat{\boldsymbol{x}} + \cos(\theta) \hat{\boldsymbol{y}}. \end{equation}
By writing 
$\dot {z_1}=-\dot {\rho } \hat {\boldsymbol {u}}_\rho -\rho \dot {\theta } \hat {\boldsymbol {u}}_\theta$, we then find the following equations for 
$\dot {\rho }$ and 
$\dot {\theta }$:
\begin{align} \dot{\rho}= \,f_1(\rho, \theta) &= \frac{k}{4 D} \text{csch}\left(\frac{{\rm \pi} \,\textrm{e}^{-\textrm{i} \theta} \rho}{D}\right) \text{csch}\left(\frac{{\rm \pi} \,\textrm{e}^{\textrm{i} \theta} \rho}{D}\right) \nonumber\\ &\quad \times\left[\cos(\theta ) \tan \left(\frac{{\rm \pi} \rho \sin (\theta )}{D}\right) \cosh^{2}\left(\frac{{\rm \pi} \rho \cos (\theta )}{D}\right) \right. \nonumber\\ &\quad \left. -\sin (\theta ) \cos ^{2}\left(\frac{{\rm \pi} \rho \sin (\theta )}{D}\right) \tanh \left(\frac{{\rm \pi} \rho \cos (\theta )}{D}\right)\right] \end{align}
\begin{align} \dot{\theta}= \,f_2(\rho, \theta) &= -\frac{k}{4 D \rho} \text{csch}\left(\frac{{\rm \pi} \,\textrm{e}^{-\textrm{i} \theta} \rho}{D}\right) \text{csch}\left(\frac{{\rm \pi} \,\textrm{e}^{\textrm{i} \theta} \rho}{D}\right) \nonumber\\ &\quad \times\left[\cos(\theta ) \cos ^{2}\left(\frac{{\rm \pi} \rho \sin (\theta )}{D}\right) \tanh \left(\frac{{\rm \pi} \rho \cos (\theta )}{D}\right)\right. \nonumber\\ &\quad \left. +\sin (\theta )\tan \left(\frac{{\rm \pi} \rho \sin (\theta )}{D}\right) \cosh^{2}\left(\frac{{\rm \pi} \rho \cos (\theta )}{D}\right)\right]. \end{align} From (A4) and (A5), we finally derive the equation for 
$\rho '=\text {d}\rho /\text {d} \theta$ as follows:
\begin{equation} \rho'= \frac{\dot{\rho}}{\dot{\theta}}= \frac{\begin{array}{l} \rho \sin (\theta ) \cos^{2}\left(\dfrac{{\rm \pi} \rho \sin (\theta )}{D}\right) \tanh \left(\dfrac{{\rm \pi} \rho \cos (\theta )}{D}\right)\\ \quad -\rho \cos (\theta ) \tan \left(\dfrac{{\rm \pi} \rho \sin (\theta )}{D}\right) \cosh ^{2}\left(\dfrac{{\rm \pi} \rho \cos (\theta )}{D}\right) \end{array}}{\begin{array}{l} \cos (\theta )\cos ^{2}\left(\dfrac{{\rm \pi} \rho \sin (\theta )}{D}\right) \tanh \left(\dfrac{{\rm \pi} \rho \cos (\theta )}{D}\right)\\ \quad +\sin (\theta )\tan \left(\dfrac{{\rm \pi} \rho \sin (\theta )}{D}\right) \cosh ^{2}\left(\dfrac{{\rm \pi} \rho \cos (\theta )}{D}\right) \end{array}}, \end{equation}
which is well defined in 
$\mathcal {A}=\{(\rho , \theta )\in \mathbb {R}^{+} \times \mathbb {R} : 0 < \rho < D/2\}$ because (a) all the elementary functions are well defined (included the function 
$\tan (\cdots )$ through the condition 
$0 < \rho < D/2$); (b) the denominator is positive (in the first term 
$\cos (\theta ) \times \tanh ({{\rm \pi} \rho \cos (\theta )}/{D}) \geq 0$ and in the second term 
$\sin (\theta) \tan ({{\rm \pi} \rho \sin (\theta )}/{D})\geq 0$) and never zero (both terms are never zero in 
$\mathcal {A}$).
In order to prove that the trajectory of vortex 
$P_1$ is a closed curve, we need to show that the function 
$\rho (\theta )$ is a continuous and periodic function. However, the integration of (A6) is a hard task to achieve. Therefore, we choose to prove that 
$\rho (\theta )$ is a continuous and periodic function without finding the exact integral of (A6). In order to achieve this goal, we first need to recall a result from mathematical analysis, which states
Theorem 1 Given a continuous and periodic function 
${f:\mathbb {R} \to \mathbb {R}}$ with period 
$T$ such that 
$\int _0^{T} f(x)\, \textrm {d}x=0$, then the primitive function of 
$f(x)$ is periodic with period 
$T$.
Having recalled Theorem 1, we now need to prove the following theorem:
Theorem 2 The primitive function 
$\rho (\theta )$ of 
$\rho '(\theta )$ (as defined in (A6)) is 
$C^{1}(\mathbb {R})$ and periodic with period at least 
$2 {\rm \pi}$.
Proof. The proof consists in three steps:
(a)
$\rho (\theta )$ is a $C^{1}(\mathbb {R})$ function;(b)
$\rho '(\theta )$ is a periodic function, at least of period $T=2{\rm \pi}$; and(c)
$\int _0^{2{\rm \pi} } \rho '(\theta )\,\textrm {d} \theta =0$.
Below we give 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 $\rho (\theta )$ is $C^{1}(\mathbb {R})$. Moreover, we can assert that the denominator of $\rho '(\theta )$ is ${\neq} 0$, or, better, it is easy to show that it is always positive for $(\rho , \theta )\in \mathcal {A}$. Indeed, the two terms in the denominator in (A6) are always positive (both for $\sin \theta$ and $\cos \theta$ positive, negative or null).(b) That
$\rho '(\theta )$ is a periodic function in fact follows directly from (A6) as
(A7)\begin{equation} \rho'(\theta+2{\rm \pi})=\rho'(\theta). \end{equation}(c) A sufficient condition to prove the last step is that the function
$\rho '(\theta )$ is an odd function in $\mathbb {R}$. The proof follows directly from (A6) after substituting $\theta$ by $-\theta$ obtaining
(A8)\begin{equation} \rho'(-\theta)=-\rho'(\theta). \end{equation}
Finally, we apply Theorem 1 to our function 
$\rho '(\theta )$ and the theorem is proved.
Theorem 2 leads hence to the conclusion that 
$\rho (\theta +2{\rm \pi} )=\rho$ and thus that the trajectory of vortex point 
$P_1$ is a closed curve.
Appendix B. Comparison of classical and quantum trajectories for periodic orbits
In this appendix we compare the trajectories of vortices 1 and 2 (cf. figure 2) determined with both the point vortex model and the GP model for two cases where classically we expect periodic orbits: (a) 
$(R/D, r/R) = (6/10, 67/100)$ (figure 9a) and (b) 
$(R/D, r/R) = (72/100, 39/100)$ (figure 9b). For 
$(R/D, r/R) = (6/10, 67/100)$ we observe for both models a periodic orbit regime, although the quantum trajectories undergo an initial translation towards the positive 
$x$ direction arising from the sudden emission of sound once we start integrating the GP equation (3.6) for 
$t\in \mathbb {R}$ as the initial condition is not stationary with respect to any frame of reference. On the contrary, for 
$(R/D, r/R) = (72/100, 39/100)$ quantum trajectories do not coincide with periodic orbits: they perform a very slow image-driven leapfrogging (each period, orbits translate towards the negative 
$x$ direction by a distance of the order of 
$2\xi$). These results show that while all regimes observed classically can effectively occur in quantum fluids, in certain circumstances some differences exists between the trajectories of corresponding vortices, as also noted in §§ 3.2 and 3.3. This may imply that, for given set of parameters 
$(R/D, r/R)$, regimes observed in classical and quantum systems may slightly differ.
Comparison of classical (red) and quantum (blue) trajectories of vortices 1 and 2 for: (a) 
$(R/D, r/R) = (6/10, 67/100)$ and (b) 
$(R/D, r/R) = (72/100, 39/100)$. For each separate set of parameters, the same physical time is shown. For 
$(R/D, r/R) = (6/10, 67/100)$, quantum trajectories are shown for more than five periods.
Open access
