N-Symmetric Interaction of N Hetons, II: Analysis of the Case of Arbitrary N

Abstract

This paper seeks and examines N-symmetric vortical solutions of the two-layer geostrophic model for the special case when the vortices (or eddies) have vanishing summed strength (circulation anomaly). This study is an extension [Sokolovskiy et al. Phys. Fluids 2020, 32, 09660], where the general formulation for arbitrary N was given, but the analysis was only carried out for $N=2$. Here, families of stationary solutions are obtained and their properties, including asymptotic ones, are investigated in detail. From the point of view of geophysical applications, the results may help interpret the propagation of thermal anomalies in the oceans.

1. Introduction

The concept of point-vortex hetons, proposed in [1,2] as the simplest model of vortex heat flow, is based on the equations of geophysical fluid dynamics [3]. It has found a wide range of applications in oceanic [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20], atmospheric [21,22,23,24] and laboratory [25,26,27,28] research. Hetons are most effectively used to model the propagation of temperature anomalies and, in particular, to explain the phenomenon of deep convection in the ocean [14,15,16,29,30,31,32,33,34,35,36,37]. Although hetons are relatively simple entities, the dynamics of a collection of hetons may exhibit both regular and chaotic motion [2,11,37,38,39,40,41,42,43].

The study of finite-core hetons has proved especially fruitful [44,45]. In particular, collections of finite-core hetons in mutual equilibrium have been found, and their stability addressed. Their nonlinear interactions are highly varied, including vortex merger, vortex splitting, filamentation, as well as vertical alignment [42,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60].

In this paper, we extend the results obtained in [61] (hereafter referred to as SKDR), for both point vortices and finite-core vortices, to the case of an arbitrary N. The paper is organised as follows. The formulation of the general problem for the motion of N hetons found in SKDR is restated in Section 2. The complete system of $2N$ differential equations is reduced to four equations for an ‘equivalent’ heton, and subsequently to just two equations for a ‘virtual’ heton. Next, in Section 3, the possible modes of behaviour for the heton system, following from the properties of the Hamiltonian, are described. Furthermore, vortex trajectories for various distinct types of motion are presented, and the class of so-called ‘choreographies’ is considered. The interesting problem of the direction of rotation of the vortices located at the vertices of the larger N-gon is discussed separately. Section 4 compares the evolution of the point-vortex hetons with the evolution of the more realistic finite-core hetons on selected examples. Section 5 discusses our main findings and their implications, as well as open problems for future research. Appendix A provides useful formulas for reference. Appendix B contains a detailed analysis of the asymptotic properties of the mean angular velocity of the vortex rotation.

2. Mathematical Model and Problem Formulation

We consider the motion of a system of $2N$ point vortices in a two-layer quasi-geostrophic flow. For simplicity, the layers are taken to have equal non-dimensional thickness $h_1=h_2=1/2$. In the conclusions, we discuss $h_1\ne h_2$ as a natural extension. Each layer contains N identical vortices, initially located at the vertices of a regular N-gon. Vortices in different layers have opposite strengths (volume-integrated potential vorticity), so that the net integrated potential vorticity is zero. The two regular N-gons have the same centre, but can have different size. Additionally, the two regular N-gons can be offset by an angle (see Figure 1) which shows two examples of initial conditions for $N=6$.

The equations are simplest if we introduce the complex coordinate $z_j^{\alpha }=x_j^{\alpha }+i{y}_j^{\alpha }$ for the $\alpha$^th vortex in layer j (with $j=1$ for the upper layer and $j=2$ for the lower layer). The vortex strengths ${\kappa }_j^{\alpha }$ are all equal to ${\kappa }_2^{\alpha }=\kappa >0$ in the lower layer and to ${\kappa }_1^{\alpha }=-\kappa <0$ in the upper layer, corresponding to ‘warm’ hetons, following [1].

The vortex motion in an unbounded horizontal domain is governed by the ordinary differential equations [57]

$$\begin{array}{ccc}\hfill {\dot{\overline{z}}}_j^{\alpha }& =& \frac{{(-1)}^j\kappa }{4\pi i}\{\sum _{\begin{array}{c}\beta =1\\ \beta \ne \alpha \end{array}}^N\frac{1}{z_j^{\alpha }-z_j^{\beta }}\left[1+\gamma |z_j^{\alpha }-z_j^{\beta }|{\mathrm{K}}_1\left(\gamma |z_j^{\alpha }-z_j^{\beta }|\right)\right]\hfill \\ & -& \sum _{\beta =1}^N\frac{1}{z_j^{\alpha }-z_{3-j}^{\beta }}\left[1-\gamma |z_j^{\alpha }-z_{3-j}^{\beta }|{\mathrm{K}}_1\left(\gamma |z_j^{\alpha }-z_{3-j}^{\beta }|\right)\right]\},\hfill \end{array}$$

(1)

where $j=1,2$ and $\alpha =1,2,\dots ,N$. Here, $\overline{q}$ denotes the complex conjugate of any quantity q, while ${\mathrm{K}}_1$ is the first-order modified Bessel function of the second kind, and $\gamma =1/\lambda$ is the inverse internal Rossby deformation length $\lambda =\sqrt{g\prime h_1{h}_2/(h_1+h_2)}/f$, with $g\prime =g\Delta \rho /{\rho }_0$ being the reduced gravity, g is the acceleration due to gravity, $\Delta \rho ={\rho }_2-{\rho }_1$ (${\rho }_1$ and ${\rho }_2$ are the densities of upper and lower layers, respectively, while ${\rho }_0$ is the mean density).

Note that Equation (1) is invariant under both translation and rotation. It follows that the two real components, $P_x$ and $P_y$, of the linear impulse and the angular impulse, M,

$$\begin{array}{c}\hfill P_x+i{P}_y=\frac{\kappa }{2}\sum _{j=1}^2{(-1)}^j\sum _{\alpha =1}^N{z}_j^{\alpha }\mathrm{and}M=\frac{\kappa }{2}\sum _{j=1}^2{(-1)}^j\sum _{\alpha =1}^N{\left|z_j^{\alpha }\right|}^2\end{array}$$

(2)

are invariants of motion. Moreover, the linear impulse is zero in each layer. A further invariant is the Hamiltonian, or interaction energy. It takes the form

$$\begin{array}{ccc}\hfill H=& -& \frac{{\kappa }^2}{8\pi }\sum _{j=1}^2\{\sum _{\begin{array}{c}\alpha ,\beta =1\\ \alpha \ne \beta \end{array}}^N\left[ln\left|z_j^{\alpha }-z_j^{\beta }\right|-{\mathrm{K}}_0\left(\gamma \left|z_j^{\alpha }-z_j^{\beta }\right|\right)\right]\hfill \\ & -& \sum _{\alpha ,\beta =1}^N\left[ln\left|z_{(3-j)}^{\alpha }-z_j^{\beta }\right|+{\mathrm{K}}_0\left(\gamma \left|z_{(3-j)}^{\alpha }-z_j^{\beta }\right|\right)\right]\}.\hfill \end{array}$$

(3)

Henceforth, we restrict attention to the N-fold symmetric situation. It follows from the invariants (2) that this symmetry is preserved for all time. Moreover, one can find a solution of (1) of the form

$$\begin{array}{c}\hfill z_j^{\alpha }=z_j\left(t\right)e^{i\frac{2\pi (\alpha -1)}{N}}=r_j\left(t\right)e^{i\left({\phi }_j\left(t\right)+\frac{2\pi (\alpha -1)}{N}\right)},j=1,2,\alpha =1,2,\dots ,N.\end{array}$$

(4)

For the initial configurations described in Figure 1, the initial angles are ${\phi }_1\left(0\right)=0$, ${\phi }_2\left(0\right)=0$ in panel (a) and ${\phi }_1\left(0\right)=0$, ${\phi }_2\left(0\right)=\pi /N$ in panel (b).

Remark 1.

Equation (2) for the angular impulse M simplifies to ${\left(r_1\right)}^2-{\left(r_2\right)}^2=-2M/\kappa N$. Since the latter is invariant, it follows that for any motion of the hetons, $r_1$ and $r_2$ must always increase or decrease simultaneously and, in particular, take extreme values synchronously. This fact is used in the analysis of the regimes of vortex motion.

Substituting a solution in the form (4) into (1) allows one to obtain evolution equations for the radial $r_j$ and angular coordinates ${\phi }_j$ of an ‘equivalent’ heton (see SKDR). Using angular momentum conservation, we can write the equations for the radial coordinates in the form:

$$\begin{array}{c}\hfill {\dot{r}}_j=\frac{R{\varrho }^{2-j}}{1-{\varrho }^2}\dot{,}j=1,2.\end{array}$$

(5)

Here, $\varrho =\frac{r_2}{r_1}$, $R\equiv r_1$, $R^2=-\frac{2\tilde{M}}{N\left(1-{\varrho }^2\right)}\left(\tilde{M}=\frac{M}{\kappa }\right)$, and

$$\begin{array}{ccc}\hfill \dot{\varrho }& =& -\frac{\kappa \left(1-{\varrho }^2\right)}{4\pi R^2}\left[\gamma R\sum _{n=0}^{N-1}\frac{sin\left(\Phi +\frac{2\pi n}{N}\right)}{G(\varrho ,\Phi ,n)}{\mathrm{K}}_1\left(\gamma RG(\varrho ,\Phi ,n)\right)\right.\hfill \\ & -& \left.\frac{N{\varrho }^{N-1}sinN\Phi }{1+{\varrho }^N({\varrho }^N-2cosN\Phi )}\right].\hfill \end{array}$$

(6)

The equations for the angular coordinates ${\phi }_j$ have the form

$$\begin{array}{ccc}\hfill {\dot{\phi }}_1& =& -\frac{\kappa }{4\pi R^2}\left[\gamma R\right(\sum _{n=1}^{N-1}sin\frac{\pi n}{N}{\mathrm{K}}_1\left(2\gamma Rsin\frac{\pi n}{N}\right)\hfill \\ & +& \sum _{n=0}^{N-1}\frac{1-\varrho cos\left(\Phi +\frac{2\pi n}{N}\right)}{G(\varrho ,\Phi ,n)}{\mathrm{K}}_1\left(\gamma RG(\varrho ,\Phi ,n)\right))\hfill \\ & -& \frac{1}{2}\left(\frac{N\left(1-{\varrho }^{2N}\right)}{1+{\varrho }^N\left({\varrho }^N-2cosN\Phi \right)}+1\right)],\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill {\dot{\phi }}_2& =& \frac{\kappa }{4\pi R^2{\varrho }^2}\left[\gamma R\varrho \right(\sum _{n=1}^{N-1}sin\frac{\pi n}{N}{\mathrm{K}}_1\left(2\gamma R\varrho sin\frac{\pi n}{N}\right)\hfill \\ & +& \sum _{n=0}^{N-1}\frac{\varrho -cos\left(\Phi +\frac{2\pi n}{N}\right)}{G(\varrho ,\Phi ,n)}{\mathrm{K}}_1\left(\gamma R\varrho G(\varrho ,\Phi ,n)\right))\hfill \\ & +& \frac{1}{2}\left(\frac{N\left(1-{\varrho }^{2N}\right)}{1+{\varrho }^N\left({\varrho }^N-2cosN\Phi \right)}-1\right)],\hfill \end{array}$$

(8)

where $\Phi \equiv {\phi }_2-{\phi }_1$. (Note that the argument $sin\frac{\pi n}{N}$ of ${\mathrm{K}}_1$ is positive for all $n=1,2,\dots ,$$N-1$). In Equations (7) and (8) above,

$$G\left(\varrho ,\Phi ,n\right)=\sqrt{1+\varrho \left[\varrho -2cos\left(\Phi +\frac{2\pi n}{N}\right)\right]}.$$

(9)

We have thus reduced the original system of $4N$ Equation (1) to a system of just three equations, (5), (7) and (8). The trajectories of the remaining $N-1$ hetons can be recovered by simple rotations through angles $2\pi n/N$, $n=1,2,\dots ,N-1$.

For a general qualitative analysis of this reduced dynamical system, it is convenient to further reduce (5)–(8) to a system of two equations for a ‘virtual’ vortex. The right-hand sides of (5) and (7)–(8) only depend on $\Phi$, the angular difference between the lower-layer and upper-layer vortices (not on ${\phi }_1$ and ${\phi }_2$). This fact allows us to reduce the system to just two equations. The first is Equation (6) above for $\varrho =r_2/r_1$, while the second is

$$\begin{array}{ccc}\hfill \dot{\Phi }& =& \frac{\kappa }{4\pi R^2{\varrho }^2}\{\frac{1}{2}\left(\frac{N\left(1-{\varrho }^{2N}\right)\left(1-{\varrho }^2\right)}{1+{\varrho }^N\left({\varrho }^N-2cosN\Phi \right)}-(1+{\varrho }^2)\right)\hfill \\ & +& \gamma R\varrho [\sum _{n=1}^{N-1}sin\frac{\pi n}{N}\left(\varrho {\mathrm{K}}_1\left(2\gamma Rsin\frac{\pi n}{N}\right)+{\mathrm{K}}_1\left(2\gamma R\varrho sin\frac{\pi n}{N}\right)\right)\hfill \\ & -& \sum _{n=0}^{N-1}\frac{\left(\left(1+{\varrho }^2\right)cos\left(\Phi +\frac{2\pi n}{N}\right)-2\varrho \right)}{G(\varrho ,\Phi ,n)}{\mathrm{K}}_1\left(\gamma RG(\varrho ,\Phi ,n)\right)\left]\right\}\hfill \end{array}$$

(10)

for $\Phi ={\phi }_2-{\phi }_1$.

The Hamiltonian (3), re-expressed in the variables ($\varrho ,\Phi$), takes the form

$$\begin{array}{ccc}\hfill H& =& \frac{{\kappa }^2}{4\pi }[ln\frac{R^2\left(1+{\varrho }^N({\varrho }^N-2cosN\Phi )\right)}{{\varrho }^{N-1}}+2\sum _{n=0}^{N-1}{\mathrm{K}}_0\left(\gamma RG(\varrho ,\Phi ,n)\right)\hfill \\ & +& \sum _{n=1}^{N-1}\left({\mathrm{K}}_0\left(2\gamma Rsin\frac{\pi n}{N}\right)+{\mathrm{K}}_0\left(2\gamma R\varrho sin\frac{\pi n}{N}\right)\right)].\hfill \end{array}$$

(11)

Here, we have used the properties of trigonometric products, see Equation (A2) of Appendix A. It is convenient to re-express Equations (6) and (10) in their Hamiltonian form

$$\begin{array}{c}\hfill \dot{\varrho }=\frac{1-{\varrho }^2}{2\kappa \varrho R^2}\frac{\partial H}{\partial \Phi },\dot{\Phi }=-\frac{1-{\varrho }^2}{2\kappa \varrho R^2}\frac{\partial H}{\partial \varrho },\end{array}$$

(12)

explicitly showing this is a (reduced) Hamiltonian system.

3. Analysis of Possible Motion Regimes

The possible motion regimes are completely determined by the structure of the Hamiltonian. The phase portrait of the Hamiltonian $H(\varrho ,\Phi )$ is shown in Figure 2 for $\tilde{M}=-1.2$ and $N=2,3$ and 4. The interval $\Phi \in [0;2\pi ]$ includes the entire periodically repeating structure of heton system.

Separatrices delimit regions of localised and unbounded motions. The intersections of the separatrices (red dashed lines) at the hyperbolic points occur at values of $\varrho$, for which the two vortex N-gons are in relative equilibrium for $\Phi =\pi /N$. Such a case is shown in Figure 1b for $N=6$. Scanning across in $\varrho$ for any fixed $\Phi$, bounded or localised motions occur only between the separatrices.

It should be noted that the system is invariant through a rotation of ${\displaystyle \frac{2\pi }{N}}$, so it suffices to consider the values $\Phi =0$ and $\Phi =\frac{\pi }{N}$ for the virtual vortex. In both cases, the equalities ${(\dots )}_i=-{(\dots )}_{N-i}$ hold, where $(\dots )$ are the expressions within the sums (see Appendix A). It follows that the right-hand side of (6) vanishes for these values of $\Phi$, and $\varrho$ reaches an extremum. This means that

$$\dot{\varrho }\left(\Phi =0,\frac{\pi }{N}\right)=0,$$

(13)

which is clearly seen in Figure 2. The same is true for $r_{1,2}$ due to Equation (5).

To analyse the extrema,

$${\varrho }^{-}=\varrho \left(\Phi =0\right),{\varrho }^{+}=\varrho \left(\Phi =\frac{\pi }{N}\right)$$

(14)

we evaluate the second derivative of $\varrho$,

$$\ddot{\varrho }{|}_{\Phi =0,\frac{\pi }{N}}=\frac{1-{\varrho }^2}{2\kappa \varrho R^2}\frac{{\partial }^{2}H}{\partial {\Phi }^2}\dot{\Phi }{|}_{\Phi =0,\frac{\pi }{N}}.$$

(15)

It is shown in Appendix B that $\frac{{\partial }^{2}H}{\partial {\Phi }^2}>0$ and $\dot{\Phi }>0$ for $\Phi =0$. Hence, the extremum ${\varrho }^{-}$ is always a minimum. It is also shown that $\frac{{\partial }^{2}H}{\partial {\Phi }^2}<0$ for $\Phi =\pi /N$. However, $\dot{\Phi }$ changes sign in this case. We have $\dot{\Phi }>0$ for ${\varrho }^{+}<{\varrho }^{*}$, $\dot{\Phi }=0$ at ${\varrho }^{+}={\varrho }^{*}$ and $\dot{\Phi }<0$ for ${\varrho }^{+}>{\varrho }^{*}$. This means that, for $0\le {\varrho }^{+}<{\varrho }^{*}$, the second derivative ${\ddot{\varrho }}^{+}<0$ and ${\varrho }^{+}$ is maximum on this interval. Then, the second derivative ${\ddot{\varrho }}^{+}>0$ for ${\varrho }^{*}<{\varrho }^{+}<1$ and ${\varrho }^{+}$ is minimum on this interval.

The extreme values of $\varrho$ can be obtained from the equation

$$H\left({\Phi }^{\pm },{\varrho }^{\pm }\right)=\overline{H},{\Phi }^{-}=0,{\Phi }^{+}=\frac{\pi }{N},$$

(16)

where $\overline{H}$ is a some fixed value of the Hamiltonian. It is, therefore, useful to analyse, in turn, the extreme values of the Hamiltonian (see also [61,62,63]):

$$H^{-}=H\left({\Phi }^{-},{\varrho }^{-}\right),H^{+}=H\left({\Phi }^{+},{\varrho }^{+}\right).$$

(17)

To this end, we examine the dependence of $H^{-}$ and $H^{+}$ on ${\varrho }^{\pm }$ in Figure 3. Here, we restrict attention to the interval $0<\varrho <1$ without loss of generality, since the situation $\varrho >1$ corresponds to swapping the two layers. It is shown in Appendix B that both $H^{-}$ and $H^{+}$ tend to infinity as $\varrho \to 0$, due to a logarithmic singularity.

For $\varrho \to 1$, the situation differs. $H^{+}$ tends to infinity (see (A24)), but $H^{-}$ tends to a finite constant (see (A25)). Since $\frac{\partial H^{-}}{\partial \varrho }$ is always negative (see Section 3) one concludes that this derivative is a monotonically decreasing function of $\varrho$. On the other hand, $\frac{\partial H^{+}}{\partial \varrho }$ changes sign at $\varrho ={\varrho }^{*}$, as shown in Appendix B. Thus, $H^{+}$ has a minimum at $\varrho ={\varrho }^{*}$.

It is also shown in Appendix B that the difference $H^{+}-H^{-}$ is always positive, but remains very small for small ${\varrho }^{+}$ and ${\varrho }^{-}$. The amplitude of the oscillation ${\varrho }^{+}-{\varrho }^{-}$ is also very small for $\varrho \ll {\varrho }^{*}$.

We next analyse the influence of the angular momentum and the number of vortices N on $H^{-}$ and $H^{+}$. We examine the possible localised (periodic) and non-localised (unbounded) types of motion. For this, it is convenient to introduce the following definitions, which are depicted graphically in Figure 3. We define ${\varrho }^{*}$ as the value of $\varrho$, for which the function $H^{+}$ has a minimum $H^{+*}$ (the dashed horizontal lines in Figure 3d). We also define ${\varrho }^{-}$ to be the intersection point between $H^{-}\left(\varrho \right)$ and either of the two levels of constant energy, $H=H_{\mathrm{d}}$ or $H_{\mathrm{u}}$, are defined so that $H_{\mathrm{d}}<H^{+*}$ and $H_{\mathrm{u}}>H^{+*}$. Finally, we define ${\varrho }^{+1}$ and ${\varrho }^{+2}$ to be the minimum and the maximum solutions to $H^{+}\left(\varrho \right)={\mathcal{H}}_{\mathrm{u}}$.

It is shown in Appendix B that ${\varrho }^{-}$ and ${\varrho }^{+2}$ always correspond to minima of $\varrho$, while ${\varrho }^{+1}$ always corresponds to a maximum of $\varrho$, for a given energy $H=\overline{H}$. It is worth mentioning that the limit $\varrho \to 1$ corresponds to the limit $r_1,r_2\to \infty$. If $\varrho =1$ at $t=0$, it follows from (12) that $\varrho =1$ for all times since $\dot{\varrho }=0$ exactly. For this special case, Equation (10) simplifies to

$$\begin{array}{ccc}\hfill \dot{\Phi }& =& \frac{\kappa }{4\pi R^2}\sum _{n=1}^{N-1}sin\frac{\pi n}{N}{\mathrm{K}}_1\left(2\gamma Rsin\frac{\pi n}{N}\right)\hfill \\ & +& 2\gamma R\sum _{n=0}^{N-1}\left|sin\left(\frac{\Phi }{2}+\frac{\pi n}{N}\right)\right|{\mathrm{K}}_1\left(2\gamma R\left|sin\left(\frac{\Phi }{2}+\frac{\pi n}{N}\right)\right|\right)-1,\hfill \end{array}$$

(18)

where $R=r_1=r_2$. This special case was studied in [1] (see also [11,38,57]). It was shown that there exists $R=R^{*}$ such that the solutions to (18) are localised and periodic, and the corresponding trajectories follow closed quasi-quadrangular patterns for $R<R^{*}$. The trajectories are unbounded for $R>R^{*}$. Thus, one deduces three possible regimes of motion (see also Figure 3 for illustration) for the virtual vortex:

(A)

Unbounded motion within the interval ${\varrho }^{-}\le \varrho \le 1$, when $H<H^{+*}$;

(B)

Localised (bounded and periodic) motion within the interval ${\varrho }^{-}\le \varrho \le {\varrho }^{+1}$, when $H>H^{+*}$;

(C)

Unbounded motion within the interval ${\varrho }^{+2}\le \varrho \le 1$, when $H>H^{+*}$.

For example, if we initialise the system with $\varrho \left(0\right)={\varrho }_0={\varrho }^{+}$ and $\Phi \left(0\right)={\Phi }_0=\pi /N$, and choose a value of the angular momentum $\tilde{M}$ to determine the initial values of the radial coordinates $r_{10}$ and $r_{20}$, the subsequent motion will be localised (bounded) if ${\varrho }^{+}<{\varrho }^{*}$, but unbounded if ${\varrho }^{+}>{\varrho }^{*}$. Both cases, in fact, correspond to the same condition: $H>H^{+*}$. It is worth mentioning that it is impossible to have $H<H^{+*}$ for $\varrho \left(0\right)={\varrho }_0={\varrho }^{+}$ and $\Phi \left(0\right)={\Phi }_0=\pi /N$, irrespective of $\tilde{M}$. We consider the case $H<H^{+*}$ separately below.

The region of unbounded motions [∞] is shaded in blue in Figure 4. The critical curve $\tilde{M}\left({\varrho }^{*}\right)$ separating bounded and unbounded motion is found by solving the equation

$$\frac{\partial H^{+}}{\partial \varrho }=\frac{\partial H(\Phi =\frac{\pi }{N},\varrho )}{\partial \varrho }=0$$

(19)

for each value of $\tilde{M}$. Notably, the asymptotes of this critical curve can be found for large values of $-\tilde{M}$, see Appendix B. We find

$${\varrho }^{*}=1-\frac{{\delta }^{*}}{N},R^{*}\approx \sqrt{-\frac{\tilde{M}}{{\delta }^{*}}}\left({\delta }^{*}\approx 1.5485,-\tilde{M}\gg 1\right)$$

(20)

and the asymptotes are indicated in Figure 4a–c by the vertical dashed lines; the boundary of the [∞] regions approaches the respective asymptotes as $-M\to \infty$.

An approximate form of the boundary of the region of parameter space corresponding to unbounded motion can also be found from Equation (A16) in Appendix B:

$$-\tilde{M}\approx 0.827\frac{N^3}{{\gamma }^2{\pi }^2}\left(1-{\varrho }^2\right)$$

(21)

for small $(1-\varrho )$. For sufficiently large N and $\gamma RG(\varrho ,\Phi ,n)\ge 1$, the value of ${\delta }^{*}$ in Equation (20) can be determined from Equation (A11) in Appendix B, and we have ${\delta }^{*}\sim 1.5485$. For $\gamma RG(\varrho ,\Phi ,n)\le 1$ the value of ${\delta }^{*}$ is smaller. It follows that ${\varrho }^{*}$ tends to the right boundary of the region when N grows and when $-\tilde{M}$ decreases. For $2<N<10$, we have checked this tendency by direct calculation. This tendency is also shown in Figure 3.

Hence, the boundary of the [∞]-region tends to the point $(\tilde{M},\varrho )=(0,1)$ along the curve (21), with a slope increasing with N.

3.1. Examples of Trajectories: Absolute and Relative Choreographies

Panels (a–f) of Figure 5 show early-time trajectories of the equivalent heton for $\tilde{M}/N=-0.28$, $N=3$, $\Phi \left(0\right)={\Phi }_0=\pi /3$ and ${\varrho }^{+}=0.864,0.878,0.8892,0.9,0.9144$ and $0.9145$. The parameter values $({\varrho }^{+},\tilde{M}/N)$ are indicated by white markers in Figure 4e. When (${\varrho }^{+},\tilde{M}/N$) = (0.864, −0.28), as shown on Figure 5a, the initial conditions belong to the red region of Figure 4e. This region, symbolised by $\left[-/+\right]$, corresponds to cases where the average angular velocity of the orbital motion of the upper layer vortex is negative, and the lower layer vortex is positive. In the general case, in addition to circular motion, the vortices also perform nutational oscillations relative to their mean position. Thus, their trajectories over long times fill concentric annular regions.

Panels (b–e) show examples of motions in the $\left[-0+/+\right]$ region (black area of Figure 4). Symbols ‘−’, ‘0’ and ‘+’, respectively, correspond to an anticyclonic displacement, the absence of any average displacement, and a cyclonic displacement of the upper layer vortex. The symbol ‘+’ after the forward slash corresponds to a cyclonic rotation in lower layer. In each panel of Figure 5, the corresponding symbol is highlighted in red. A characteristic feature of this class of motions is that the sign of the total angular velocity (not averaged) of the upper-layer vortices alternates (nutational trajectories are loop-shaped). In panel (c), the outer vortex only performs periodic nutational motions, while the inner vortex periodically rotates along a fixed quasi-triangular curve (the corresponding marker in Figure 4a lies on the white line). In reality, we have a six-vortex structure, in which three peripheral vortices in the upper layer rotate along individual closed trajectories shifted relative to each other by an angle $\pi /3$, and three vortices in the lower layer move along a common central quasi-triangular trajectory one after the other with the same angular displacement. For this reason, the trajectories in panel (c) are called ‘absolute complex choreographies’.

Remark 2.

Simó [64] introduced the concept of choreography applied to rigid body dynamics. Choreographies can be (a) absolute, if the trajectories are closed in a fixed coordinate system; (b) relative, if the trajectories are closed in a uniformly rotating or translating coordinate system; (c) simple, if all of the material points move along a common trajectory; and (d) complex, if at least one of the points moves along a separate trajectory. The concept of choreography in vortex theory was used for the first time in [65]. Modern studies of choreography as applied to 2D vortex dynamics are presented in [65,66,67,68,69], and for a two-layer rotating fluid in [56,57,61].

It should be noted that the values of the values of the parameter ${\varrho }^{+}$ in panels (b–d) differ only slightly and, therefore, the markers in Figure 4e are almost undistinguishable. The last pair of trajectories of this type in Figure 4e corresponds to initial conditions close to the border of the black region. There, the loops almost disappear, and the trajectory of the outer vortex approaches a circular form. An increase in ${\varrho }^{+}$ by just 0.0001 shifts the initial conditions into the region of unbounded motions, when three two-layer pairs (hetons) scatter in different directions—see Figure 5f. It should be noted that the example of absolute choreographies, shown in Figure 5c, is not unique. Trajectories of this type correspond to points belonging to the white lines in Figure 4a–c,e. The case $N=2$ was studied in detail in SKDR.

Figure 6 shows a gallery of absolute choreographies for $2\le N\le 6$, corresponding to the maximum values of ${\varrho }^{+}$ (the closed curves or orbits were obtained numerically). We call these periodic trajectories the ‘limiting’ ones. The parameters ${\varrho }^{+}$ and $\tilde{M}/N$ for the cases $N=2,3$ and 6 correspond to the white markers in the upper right corners of Figure 4a–c.

Figure 6 shows that the size of internal N-shaped (N-gons with smoothed corners) choreographies increase with increasing of N.

We next consider the topological properties of the choreographies for the case $N=6$ ($\Phi \left(0\right)={\Phi }_0=\pi /6$) in more detail in Figure 7. We show a family of choreographies over an interval for ${\varrho }^{+}\in [0.60;0.97]$ with an increment of 0.01. This interval is the green portion of the white line in Figure 4c. Numerical experiments show that as ${\varrho }^{+}$ decreases, the quasi-hexagonal choreographies decrease in size in the lower layer and the trajectories become nearly circular. In the upper layer, the situation is more complicated. For ${\varrho }^{+}\ge 0.97$, the vortices move along thin drop-shaped trajectories, whose innermost tips point towards the origin. For ${\varrho }^{+}$ = 0.96, figure-of-eight orbits begin to form due to the emergence of a second closed loop to the left of the tip. Such doubly connected patterns occur in the interval ${\varrho }^{+}\in [0.93,0.96]$. In this interval, as ${\varrho }^{+}$ decreases, the left part of the ‘eight’ grows, while the right part shrinks. Finally, for ${\varrho }^{+}$ = 0.93, the right part completely disappears, and now the drop-shaped trajectory has a tip pointing away from the origin. As we decrease ${\varrho }^{+}$ further, the size of the drop first increases to reach a local maximum at ${\varrho }^{+}\approx 0.90$, and then decreases, eventually shrinking to a point.

We next consider the evolution of the system close to the figure-of-eight stationary states for ${\varrho }^{+}=0.95$ (purple curves in Figure 7) and for ${\varrho }^{+}=0.97$ (black curves). The next three figures are devoted to ${\varrho }^{+}=0.95$, followed by another devoted to ${\varrho }^{+}=0.97$.

Figure 8 shows a series of vortex trajectories in the vicinity of an equivalent heton with $N=6$ and ${\varrho }^{+}=0.95$. Panel (b) shows the absolute choreography of the [-0+/+] type, for which the vortex trajectory of the upper layer has a figure-of-eight shape, and panels (a) and (c) show the absolute vortex trajectories of types [-0+/+] and [-0+/+],, respectively, when the upper layer vortex moves along a loop-like trajectory, unwinding in the anticyclonic and cyclonic directions, respectively. This behaviour is found in the black region of Figure 4c.

We construct the relative choreographies for these cases in Figure 8(a₁,c₁). The relative choreographies are the closed trajectories obtained in a reference frame rotating with the average angular velocity $\omega$ of the vortex in the upper layer. One may see that the relative choreographies have loop-like absolute upper-layer trajectories (found inside the black region of Figure 4c) also belong to a type of figure-of-eight trajectory. The figure-of-eight relative choreographies also occur in some parts of the red and yellow areas of the diagram, where the upper-layer vortices move along zigzag trajectories—see Figure 9.

Further away from the stationary state, it is possible to obtain drop-shaped relative choreographies, both in the red region (for $\tilde{M}/N=-0.2$) and in the yellow region (for $\tilde{M}/N=-0.5$) in Figure 4c, as shown in Figure 10.

We next consider the second class of motions in the vicinity of the stationary state, when the absolute choreography for the vortex of the upper layer has a drop-like shape, i.e., for ${\varrho }^{+}=0.97$. Results are presented in Figure 11, where all absolute trajectories of the upper-layer vortices are loop-shaped. The initial conditions are in the black region of Figure 4c, outside the yellow [+/+]-region. All relative choreographies are drop-shaped.

Let us point out another new type of absolute choreography in the N-symmetric dynamics of N hetons. In this choreography, trajectories are circular synchronous rotations of the lower-layer and upper-layer vortices around a common centre. Clearly, these movements are of the [+/+] type (the yellow areas in Figure 4). For a circular motion, the right-hand side of Equation (5) must be zero, i.e., $\varrho \equiv 0$, and we must have ${\omega }_1={\dot{\phi }}_1={\omega }_2={\dot{\phi }}_2$ (see Equations (7) and (8)). These conditions are satisfied for a discrete set of parameters in the plane (${\varrho }^{+},\tilde{M}/N$) in a narrow neighbourhood of the yellow region near the boundary of the region [∞]. Two examples of such circular choreographies for $N=5$ and 6 are shown in Figure 12. Markers in the figure only show the initial location of the two vortices of the equivalent heton. In reality, five (respectively, six) anticyclones/cyclones in the upper/lower layer move counter-clockwise, with the same angular velocity, thereby preserving the regular pentagonal or hexagonal shape.

3.2. Direction of Rotation of the External Vortices

Since we assume throughout that $r_1>r_2$ (as in Figure 1), the external N-gon is in the upper layer. Thus, ‘external/internal’ and ‘upper-layer/lower-layer’ are synonymous in this context. The internal vortices always maintain their direction of rotation, while the external ones can change it.

The direction of rotation of the external vortices depends on the region of the parameter space. The external vortices are expected to rotate in the opposite (anticyclonic) direction if the interaction between the two layers is sufficiently weak. In the parameter space $({\varrho }^{+};\tilde{M})$ in Figure 4, this anticyclonic motion occurs in the red [$-/+$] region. To determine the boundary of this region, we determine the value of the angular impulse $\tilde{M}\left({\varrho }^{+}\right)$ at which the angular velocity of the external vortices vanishes. This leads to the implicit equation:

$$\begin{array}{ccc}\hfill {{\dot{\phi }}^{+}}_1& =& \frac{\kappa }{4\pi R^2}[\frac{N\left(1-{\varrho }^N\right)+1+{\varrho }^N}{2\left(1+{\varrho }^N\right)}-\sum _{n=1}^{N-1}\gamma Rsin\frac{\pi n}{N}{\mathrm{K}}_1\left(2\gamma Rsin\frac{\pi n}{N}\right)\hfill \\ & -& \sum _{n=0}^{N-1}\frac{1-\varrho cos\frac{\pi (1+2n)}{N}}{G(\varrho ,\frac{\pi }{N},n)}\gamma R{\mathrm{K}}_1\left(\gamma RG(\varrho ,\frac{\pi }{N},n)\right)]=0.\hfill \end{array}$$

(22)

Analytical progress can be made assuming $\varrho \ll 1$. Then, Equation (22) reduces to

$${{\dot{\phi }}^{+}}_1=\frac{\kappa }{4\pi R^2}\left[\frac{N+1}{2}-\frac{1}{2}\sum _{n=1}^{N-1}2\gamma Rsin\frac{\pi n}{N}{\mathrm{K}}_1\left(2\gamma Rsin\frac{\pi n}{N}\right)-N\gamma R{\mathrm{K}}_1\left(\gamma R\right)\right]=0.$$

(23)

We can estimate the solution by manipulating the sum in the equation. On the one hand, if we neglect the sum, we obtain the approximate value of $R=1.257$, which is smaller than the actual solution to Equation (23). On the other hand, we can find an upper bound for the sum in question:

$$\begin{array}{ccc}& & \sum _{j=1}^{N-1}\sum _{n=n_j}^{n_{j+1}}2\gamma Rsin\frac{\pi n}{N}{\mathrm{K}}_1\left(2\gamma Rsin\frac{\pi n}{N}\right)\hfill \\ & <& \sum _{j=1}^{N-1}\left(n_{j+1}-n_j\right)2\gamma Rsin\frac{\pi n_j}{N}{\mathrm{K}}_1\left(2\gamma Rsin\frac{\pi n_j}{N}\right),n_1=0,n_2=1,\dots ,n_N=N-1,\hfill \end{array}$$

(24)

from which we obtain the upper bound $R=2$ for the solution to (23). Thus, we obtain the following bounds for the solution:

$$1.257<\gamma R\left({\varrho }^{+}=0\right)=\gamma \sqrt{\frac{-2\tilde{M}}{N}}<2,\tilde{M}=-\frac{N}{2}{R}^2\left({\varrho }^{+}=0\right).$$

(25)

Numerical solutions to Equation (23) for $50\le N\le 100$ suggest that $R\left({\varrho }^{+}=0\right)$ tends to $R_0^{\infty }\approx 1.72$.

From this result, we conclude that the boundary of the [$-/+$] region corresponds to

$$\tilde{M}=-\frac{N}{2}{\left(R_0^{\infty }\right)}^2\approx -1.4792N,$$

(26)

when $\varrho \to 0$, and decreases with N. The limiting lower boundary $\tilde{M}/N\approx -1.4792$ of the red regions is shown in Figure 4a–c by the dashed horizontal lines.

We next consider the direction of rotation of the internal (lower-layer) vortices. Their angular velocity at the minimum of $r_2=R\varrho$ is found from

$$\begin{array}{ccc}\hfill {{\dot{\phi }}_2}^{-}& =& \frac{\kappa }{4\pi R^2{\varrho }^2}[\sum _{n=1}^{N-1}\gamma R\varrho sin\frac{\pi n}{N}{\mathrm{K}}_1\left(2\gamma R\varrho sin\frac{\pi n}{N}\right)\hfill \\ & +& \gamma R\varrho \sum _{n=0}^{N-1}\frac{\varrho -cos\frac{2\pi n}{N}}{G(\varrho ,0,n)}{\mathrm{K}}_1\left(\gamma RG(\varrho ,0,n)\right)+\frac{N-1+\left(N+1\right){\varrho }^N}{2\left(1-{\varrho }^N\right)}].\hfill \end{array}$$

(27)

We can see that ${{\dot{\phi }}_2}^{-}$ is always positive. In Equation (27), only the second sum can be negative, but at sufficiently large R this sum is exponentially small. For small R one can use (A6) and approximate ${\mathrm{K}}_0$ as a natural logarithm and use (A2) to obtain the estimate $\frac{N{\varrho }^N}{\left(1-{\varrho }^N\right)}>0$ for this sum.

At the maximum distance of $r_2$, the value of ${\dot{\Phi }}^{+}$ vanishes only at the edge of the blue [∞] region. It follows that ${{\dot{\phi }}^{+}}_2$ must be opposite in sign to ${{\dot{\phi }}^{+}}_1$ at this distance, i.e., it is also positive at maximum distance. This confirms that the internal vortices always rotate in the same direction.

We next analyse the boundary of the yellow [$+/+$] region shown in Figure 4, where the vortices in the two layers rotate in the same direction. Along this boundary, given by the curve $\tilde{M}\left({\varrho }^{+}\right)$, the angular velocity of the external vortices vanishes at the minimum of $r_1=R$. It thus cannot change sign for any $r_1$. The equation for $\tilde{M}\left({\varrho }^{+}\right)$ is

$$\begin{array}{ccc}\hfill {{\dot{\phi }}_1}^{-}=& -& \frac{\kappa }{4\pi R^2}[\sum _{n=1}^{N-1}\gamma Rsin\frac{\pi n}{N}{\mathrm{K}}_1\left(2\gamma Rsin\frac{\pi n}{N}\right)\hfill \\ & +& \gamma R\sum _{n=0}^{N-1}\frac{1-\varrho cos\frac{2\pi n}{N}}{G(\varrho ,0,n)}{\mathrm{K}}_1\left(\gamma RG(\varrho ,0,n)\right)-\frac{N+1+\left(N-1\right){\varrho }^N}{2\left(1-{\varrho }^N\right)}].\hfill \end{array}$$

(28)

This curve coincides with the boundary of the red $[-/+]$ region shown in Figure 4, in the limit $\varrho \ll 1$, and reaches the same limiting point given by Equation (26). There are no solutions for $\tilde{M}\to -\infty$ or $\tilde{M}\to -0$. It follows that the boundaries of the yellow $[+/+]$ and the blue $[\infty ]$ regions intersect at a finite value of $\tilde{M}$ (as has been confirmed numerically).

There is a small region between the $[-/+]$ and $[+/+]$ regions, shown in black in Figure 4, where angular velocity of external vortices changes sign over the period $T_2$ of oscillation. The white curve lying inside this region indicates where the time-averaged angular velocity, ${\langle {\dot{\phi }}_1\rangle }_{T_2}$, vanishes. Along this curve, the external vortices only nutate, with no net rotation over a period.

We next consider the characteristics of the localised, bounded vortex motions. First, we examine angular velocity profiles in a few examples. Figure 13 shows the initial upper-layer and lower-layer profiles ${\omega }_{1,2}^{+}={\dot{\phi }}_{1,2}^{+}$ (the right-hand sides of Equations (7) and (8) at $t=0$), for $\Phi =\pi /N$ and $N=3$. In Figure 13a the profiles vary with ${\varrho }^{+}=r_2/r_1$ for two values of $R=r_1$, while in (b) the profiles vary with R for two values of ${\varrho }^{+}$. In panel (a), the angular velocity profile does not change sign over the entire interval of the parameter ${\varrho }^{+}$, but the signs are opposite for the upper-layer and lower-layer vortices when $R=1$ (they move in the direction corresponding to their strength). However, when $R=4$, both angular velocities are positive, due to the dominant role of the cyclonic lower layer. In Figure 13b in both cases (${\varrho }^{+}=0.6$ and ${\varrho }^{+}=0.8$), ${\omega }_1^{+}$ changes sign as R increases, when the upper-layer vortices fall under the prevailing influence the lower-layer vortices.

Following SKDR, we introduce a characteristic of the direction of the upper-layer vortex rotation:

$${\delta }_{max}=\frac{T_1}{T_0}=\frac{{\langle {\omega }_2-{\omega }_1\rangle }_{T_0}}{{\langle {\omega }_1\rangle }_{T_1}},$$

(29)

where $T_0$ is the time it takes a vortex to move from its maximum radius to its minimum radius and back, $T_1$ is the upper-layer vortex rotation period around the origin, and ${\langle \dots \rangle }_T$ denotes the time average over the interval $[0,T]$. For further insight, the parameter $\delta =T_1/T_2$, which characterises the ratio of the rotational periods of the vortices in the upper and lower layers, is also useful.

Figure 14 shows examples of the distributions ${\delta }_{max}\left({\varrho }^{+}\right)$, for several values of the angular impulse, $\tilde{M}$, and for $N=6$. Here, the red line curve shows ${\delta }_{max}\left({\varrho }^{+}\right)$ when $\tilde{M}/N=-0.1326$, corresponding to the white marker in the upper right corner in Figure 4c, as well as the red vortex trajectories (absolute choreography) in Figure 6, in which ${\varrho }^{+}=0.9854$. The nonmonotonic character of ${\delta }_{max}\left({\varrho }^{+}\right)$ is consistent with the [$-/+$]-type of motion, but the monotonically decreasing branch of this distribution corresponds to the [$+/+$]-type of motion.

According to Figure 4c, when $\tilde{M}/N=-1$ only the [$+/+$] type of bounded motion occurs (the green line on Figure 14). The two cases $\tilde{M}/N=-0.3$ and $-0.4$ by contrast contain both bounded and unbounded types of motion. The black markers in Figure 14 indicate the initial conditions for the periodic motions shown in Figure 15 and Figure 16. We see that ${\delta }_{max}$ varies continuously in certain intervals of the parameters ${\varrho }^{+}$ and $\tilde{M}/N$, with the exception of the resonance values of ${\varrho }^{+}$ when ${\omega }_1=0$. Of particular interest are the cases with integer values of ${\delta }_{max}$; they correspond to closed periodic trajectories of the vortex of the upper layer with an integer number of maxima. In view of Remark 1, in this case, the vortex in the lower layer must also move along a closed trajectory, which must also have an integer number of maxima. Thus, the value of $\delta$ for such cases is always rational or an integer. Figure 15 shows four examples of stationary periodic trajectories for an equivalent heton. In all cases, the closed shape of the vortex trajectory of the upper layer is such that, as noted above, the number of maxima coincides with value of ${\delta }_{max}$: at ${\delta }_{max}$ = 2 it resembles an ellipse symmetrically flattened above and below. In the other cases, the trajectories are quasi ${\delta }_{max}$-gonal in structure. In the lower layer, the situation as a whole is more complicated. For ${\delta }_{max}$ = 2, the trajectory is closed with three maxima of the radial coordinate $r_2$ and three self-intersection points. For ${\delta }_{max}$ = 3, the trajectory has an oval shape. The number of maximum increases as ${\delta }_{max}$ increases.

Figure 16 shows the periodic trajectories for ${\delta }_{max}$ = 20, $\tilde{M}/N=-0.1326,{\varrho }^{+}=0.96550$ (left panel) and $\tilde{M}/N=-0.3,{\varrho }^{+}=0.96447$ (right panel). In the first case, the corresponding marker in Figure 14 lies on the red non-monotonic curve, and in the second it lies on the blue monotonic one. Since the values ${\varrho }^{+}$ are very close for both cases, the markers in Figure 14 are visually indistinguishable. It is important to note that in the first case the solution belongs to the [$-/+$] type of motion; that is, when the upper-layer vortices rotate, on average, anticyclonically and the lower-layer vortices rotate cyclonically. In the second case, all vortices rotate cyclonically (region [$+/+$]).

Remark 3.

Recall that Figure 15 and Figure 16 show the motion of an equivalent heton (i.e., one vortex per layer). The motion of the complete vortex structure is obtained by N-fold symmetry.

Figure 17 shows the trajectories of all the vortices for the cases shown in Figure 15 with $N=6$. The combined trajectories are 6-fold symmetric, unlike that of the equivalent heton, but are also highly tangled. We discuss the cases with different integer values of ${\delta }_{max}$ in turn.

For ${\delta }_{max}=2$, there are three distinct trajectories in the upper layer (plotted in black, red and blue) along which vortices (1,4), (2,5) and (4,6) move, respectively. In the lower layer, we observe two distinct trajectories. Vortices 1, 3 and 5 move along the black trajectory, while vortices 2, 4 and 6 move along the red one.

For ${\delta }_{max}$ = 3, now there are only two distinct trajectories in the upper layer, with vortices 1, 3 and 5 moving along the trajectory plotted in black, and vortices 2, 4 and 6 along the trajectory plotted in red. In the lower layer, there are three distinct trajectories along which vortices (1,4), (2,5) and (4,6) move, respectively. This behaviour is exactly the reverse of what was found for ${\delta }_{max}=2$. The fact that the vortices move along independent closed curves is characteristic of the vortex motion for low-mode periodic solutions (low integer values of ${\delta }_{max}$) for $N>2$.

For ${\delta }_{max}$ = 4, as in the first case, three distinct trajectories are observed in the upper layer, while now all vortices of the lower layer move along a single trajectory. During the full period of revolution around the centre, the radial coordinate of the vortices reaches a maximum twelve times.

The case ${\delta }_{max}=5$ stands out from the previous ones, because now, in the upper layer, each of the six vortices moves along its own individual trajectory. This is because $N=6$ and ${\delta }_{max}=5$ are co-prime. Each trajectory is plotted in a different colour in Figure 17. In the lower layer, as before, all the vortices follow one another along one common trajectory, and now the the radial coordinate reaches a maximum thirty times per period.

Thus, in all cases with integer values of ${\delta }_{max}$, all trajectories are absolute choreographies. In the first two cases, the choreographies are convoluted in both layers, and in last two cases, the choreographies are convoluted in the upper layer, but simple in the lower layer.

Figure 18 shows three more examples of the trajectories of all six two-layer pairs of vortices for $\tilde{M}/N=-0.1396$. The trajectories of the vortices of the upper layer are plotted in red, while those of the lower layer are plotted in blue. The central part of the figure shows the trajectories of all the vortices for ${\varrho }^{+}=0.96550$ corresponding to $R=1.97761$ and $\delta =7/4$, as previously shown in Figure 16. In this case, vortices (1,6), (2,5) and (3,6) of the upper layer move along three distinct trajectories, while the vortices of the lower layer move along a single trajectory. In the second case, for ${\varrho }^{+}=0.8954$ corresponding to $R=3.024729$ and $\delta \to \infty$, the vortices of the lower layer follow a single trajectory resembling an hexagon with smoothed edges. The vortices of the upper layer move along distinct, closed, peripheral trajectories, as previously shown in Figure 6. In the third case, for ${\varrho }^{+}=0.8960$ corresponding to $R=3.088396$ (i.e., a small increase in ${\varrho }^{+}$ from the second case), the vortices scatter as six heton pairs. This is justified by the fact that the second case lies on the boundary between the regions of bounded and unbounded motion shown in the regime diagram in Figure 4.

4. Finite-Core Hetons

In this section, we present some results for the interaction of finite-core hetons (i.e., hetons comprised vortices with finite area). The hetons have uniform PV, $q=\pm 2\pi$, without loss of generality. Considering finite-core hetons introduces a new length scale to the problem, namely the size of the vortices. Vortices have initially a circular horizontal cross-section of radius $R_v$. Finite-core vortices have a shape and can deform. The deformation affects the velocity field the vortices induce on each other. Vortices can also strongly interact and merge when close together. No symmetry is imposed in the nonlinear evolution, and perturbations (including asymmetric ones) are allowed to develop. Simulations are run using Contour Surgery with the method’s standard set-up [70]. Equations are marched in time using a fourth-order Runge–Kutta scheme with a time step set by the vortices PV: $\Delta t=\pi /\left(10q\right)$.

In the first set of numerical experiments, we revisit the point-vortex case for $N=3$, $\tilde{M}/N=-0.28$ and ${\varrho }^{+}=0.8892$ shown in Figure 5c, but now for finite-core vortices. We first consider vortices of small to moderate radius, specifically $R_v=0.1,0.2,0.3,0.4$ and $0.5$. We keep $\tilde{M}/N=-0.28$ fixed and we determine the value of ${\varrho }^{+}\left(R_v\right)$ corresponding to quasi-closed trajectories (absolute choreographies) for the vortex centres, by trial and error. The values of ${\varrho }^{+}\left(R_v\right)$ thereby obtained are listed in

Table 1. Values of ${\varrho }^{+}$ for an absolute choreography for $N=3$ and $\tilde{M}/N=-0.28$ for various initial radii $R_v$.

$R_v$	0.05	0.1	0.2	0.3	0.4	0.5
${\varrho }^{+}$	0.8897	0.8904	0.8935	0.8978	0.9025	0.9078

.

The values of ${\varrho }^{+}\left(R_v\right)$ are remarkably close to the value obtained for point vortex hetons, and the maximum relative departure for the largest hetons considered, $R_v=0.5$, is only about $2\%$.

Figure 19 shows both the long-time trajectories of the centre of the first vortex in each layer (in black), as well as a the shape of the bounding contour for all six vortices at various times (in blue for the anticyclones in the upper layer, and in red for the cyclones in the lower layer).

We next measure the vortex deformation by evaluating the semi-axis lengths a and b, with $a>b$ without loss of generality. We perform this by fitting ellipses having the same second-order spatial moments to the actual vortex bounding contours (see, e.g., [63]). As $R_v$ is increased while keeping the other parameters fixed, as expected, the vortex deformation becomes more significant, as shown in Figure 20 for the first vortex of the upper layer. Instead of moving along a unique trajectory, the vortices move along criss-crossing trajectories without deviating far from a mean line overall. This is a consequence of the deformation of the vortices shown in Figure 20. The apparent thickness of the black lines (increasing with $R_v$) showing the vortex trajectories, in fact, indicates the spread of the trajectories.

When $R_v$ is increased further, the relative distance $\ell /R_v$, where ℓ is a typical scale of the distance between the vortex centres in a given layer, decreases, and vortices are likely to strongly interact. Figure 21 shows such a strong interaction for $N=3$, $\tilde{M}/N=-0.28$ and $\varrho =0.9073$, i.e., close to the corresponding absolute choreography for larger vortices $R_v=0.8$. The three cyclonic vortices of the lower layer deform strongly, and merge first in the centre of the domain. Three secondary satellite vortices form and are ejected away from the centre. This is a consequence of the invariance of the angular impulse M, and the migration of positive PV towards the centre of the lower layer.

We next consider the interaction between six hetons close to the regime of absolute choreography for $\tilde{M}/N=-0.3$. For $R_v=0.2$, we recover the absolute choreography for ${\varrho }^{+}=0.965$ (see Figure 22a). With increasing $R_v$, obtaining absolute choreographies for long times becomes difficult: numerical results indicate that the configuration is or becomes unstable. For $R_v=0.5$ and ${\varrho }^{+}=0.9678$, when $t\le 280$, the vortices move along quasi-closed trajectories, as seen in Figure 22b. They however depart from them at later times, as shown in Figure 22c.

5. Discussion and Concluding Remarks

In SKDR, we analysed the motion of an N-symmetric two-layer system of $2N$ vortices for the special case N = 2. In the present paper, we have extended these results to the case of arbitrary N. Diagrams of all possible regimes of motion have been established and analysed. Numerous examples of vortex trajectories have been considered, mainly for N = 3 and N = 6. Particular attention has been paid to the study of the properties of relative and absolute choreographies, both simple and complex [64]. Note that all the results of the calculations and the analysis presented in this article for $r_1<r_2$ are new, and have not previously been considered by other authors. We can only note the work [2], where test calculations were carried out for the cases N = 2 and N = 6, but at $r_1=r_2$.

The concept of a two-layer quasi-geostrophic heton is attractive, first of all, because it is the simplest model for the horizontal transfer of temperature anomalies. The vortex structures considered in this work, with anticyclones in the upper layer and cyclones in the lower layer, depress the interface between the layers. Thus, any vortex motion is accompanied by the transport of a fluid column warmer than the surrounding medium. So, the following ‘thermal’ interpretation can be given to the numerical experiments shown in Figure 18.

• ${\varrho }^{+}$ = 0.9655: A regular mixing of warm water in both layers occurs inside an annular region, somewhat wider than the area occupied by the trajectories. The total transfer of warm fluid occurs in the cyclonic direction (type [+/+]).
• ${\varrho }^{+}$ = 0.8954: All the vortices of the lower layer move move along a single trajectory in the cyclonic direction, and all the vortices of the upper layer move in the anticyclonic direction along another peripheral closed curves (type [0/+]). Both sets of vortices carry with them columns of fluid that are warmer than their surroundings.
• ${\varrho }^{+}$ = 0.8960: the twelve vortices form six two-layer pairs with tilted axes, which scatter radially while preserving the N-fold symmetry of the system, and carrying warm water (motion type [∞]).

In the general case, bounded motions of N-symmetric warm hetons are accompanied either by cyclonic transport of warm water column from surface to bottom, or cyclonic transport in the lower layer and anticyclonic in the upper layer inside an annular area. Unbounded motions expel the warmer fluid radially, away from its initial location. A similar transport of cold water would occur if the hetons were cold instead.

Finally, we have shown a few selected examples that the idealised point vortex trajectories are reproduced for more realistic finite core hetons, so long as the vortices are not too large. Large vortices may strongly interact in the nonlinear regime.

In this paper, we have not addressed the stability of stationary configurations. This will be addressed in a future paper. We expect that the motions of unstable configurations may also be irregular, breaking the initial symmetry of the heton.