1. Introduction
Planetary fluids on large scales, often dominated by rapid rotation and stable stratification, ubiquitously self-organize into coherent vortices. In appropriately scaled vertical coordinates the vortices take the form of roughly spherical patches of potential vorticity (McWilliams Reference McWilliams1984; McWilliams & Weiss Reference McWilliams and Weiss1994; McWilliams, Weiss & Yavneh Reference McWilliams, Weiss and Yavneh1994; Dritschel, de la Torre Juárez & Ambaum Reference Dritschel, de la Torre Juárez and Ambaum1999; Carton Reference Carton2001; Viúdez Reference Viúdez2016). A population of physically well-separated coherent vortices, a ‘vortex gas’, can be approximated using a hierarchy of idealized models (Kida Reference Kida1981; Melander, Zabusky & Styczek Reference Melander, Zabusky and Styczek1986; Meacham Reference Meacham1992; Meacham, Morrison & Flierl Reference Meacham, Morrison and Flierl1998; Miyazaki, Furuichi & Takahashi Reference Miyazaki, Furuichi and Takahashi2001; Dritschel, Reinaud & McKiver Reference Dritschel, Reinaud and McKiver2004). The simplest member of the hierarchy is the point-vortex model (Helmholtz Reference Helmholtz1867; Kirchhoff Reference Kirchhoff1876; Aref Reference Aref2007) in which the vorticity is considered to be concentrated at a collection of isolated points in space. Point-vortex models are useful idealizations for geophysical and astrophysical turbulence and are used to model atmospheres and oceans on the Earth, planets and exoplanets and astrophysical disks (Aref & Stremler Reference Aref and Stremler2001; Carton Reference Carton2001; Lucarini et al. Reference Lucarini, Blender, Herbert, Ragone, Pascale and Wouters2014; Hirt et al. Reference Hirt, Schielicke, Múller and Névir2018; Abrahamyan Reference Abrahamyan2020).
The rotation in planetary fluids is characterized by a non-dimensional Rossby number, 
$\textit {Ro} = U/fL$, the ratio of the planetary rotation time scale 
$1/f$ to the advective time scale 
$L/U$, where 
$U$ is the advective velocity scale, 
$L$ is a characteristic length scale, 
$f = 2\varOmega \sin (\theta _{lat})$ is the Coriolis parameter, 
$\varOmega$ is the planetary rotation rate and 
$\theta _{lat}$ is the latitude. A flow whose advective time scale is much longer than its rotation time scale has 
$\textit {Ro} \ll 1$. Large scale midlatitude flows in the Earth's atmosphere typically have 
$\textit {Ro} \sim 0.1$, flows in the Earth's oceans typically have 
$\textit {Ro} \sim 0.01$, while flows on planets and exoplanets often have similarly small Rossby numbers (Showman, Cho & Menou Reference Showman, Cho and Menou2010). The small Rossby number renders geostrophic balance the dominant physics: Coriolis forces approximately balance horizontal pressure gradient forces. The leading-order asymptotic theory for rapidly rotating stratified flows is quasigeostrophy (QG), which results from an asymptotic expansion of the primitive equations for small Rossby number. QG can be derived for flows with many different vertical density structures. Most well known is QG with a stably stratified background density, either in layers, or with a continuous stratification. Two notable features of stably stratified QG are that there is no vertical geostrophic velocity and there is an antisymmetry between cyclones and anticyclones.
Balance models go beyond QG and include ageostrophic (AG) dynamics as an 
$O(\textit {Ro})$ correction to QG, retaining the dominant QG balance as the lowest-order dynamics. There is a large literature on balance models and many different balance models have been developed, (e.g. McWilliams & Gent Reference McWilliams and Gent1980; Allen Reference Allen1993; Holm Reference Holm1996; Muraki, Snyder & Rotunno Reference Muraki, Snyder and Rotunno1999, and references therein). These models have 
$O(\textit {Ro}^2)$ differences and thus are, formally, equivalently correct through 
$O(\textit {Ro})$. Different approaches to balance models choose different quantities as a distinguished variable. For example, Allen (Reference Allen1993) chooses pressure. Here, we investigate the balance model of Muraki et al. (Reference Muraki, Snyder and Rotunno1999) (MSR1999), QG
$^{+1}$ , which chooses potential vorticity 
$q$ as its distinguished variable. This makes it a natural choice for formulating balanced point vortices. Balance models differ qualitatively from QG in that the AG dynamics includes a non-zero 
$O(\textit {Ro})$ vertical velocity and breaks the cyclone–anticyclone antisymmetry. Despite being small, both of these effects are important for many aspects of planetary fluid dynamics. There have been a number of studies of finite-sized coherent vortices in balance models as well as in primitive equation models with varying Rossby number (e.g. Tsang & Dritschel Reference Tsang and Dritschel2015; McKiver & Dritschel Reference McKiver and Dritschel2016; Mahdinia et al. Reference Mahdinia, Hassanzadeh, Marcus and Jiang2017; Reinaud & Dritschel Reference Reinaud and Dritschel2018; McKiver Reference McKiver2020; Sokolovskiy, Carton & Filyushkin Reference Sokolovskiy, Carton and Filyushkin2020, and references therein).
Point-vortex models have a long history in two-dimensional (2-D) fluid dynamics (see Aref (Reference Aref2007) for a review). They are also well developed in layered and continuously stratified QG dynamics (Gryanik Reference Gryanik1983, Reference Gryanik1991; Gryanik et al. Reference Gryanik, Doronina, Olbers and Warncke2000; Gryanik, Sokolovskiy & Verron Reference Gryanik, Sokolovskiy and Verron2006; Reznik & Kizner Reference Reznik and Kizner2007). Due to the constraints of QG dynamics, QG point-vortex models are limited to modelling physical situations where the vertical velocity and cyclone–anticyclone antisymmetry are unimportant.
Here, we find point-vortex solutions to the nonlinear balance model of MSR1999, QG
$^{+1}$ . The solutions have many similarities to 2-D and QG point vortices, but go beyond them by inheriting from balance models 
$O(\textit {Ro})$ vertical velocities and broken cyclone–anticyclone antisymmetry. In § 2, we review the aspects of QG
$^{+1}$ theory relevant for the development of QG
$^{+1}$ point vortex dynamics. In § 3 we explore vortex-gas solutions of QG
$^{+1}$ . Section 4.1 finds the equations of motion for QG
$^{+1}$ point vortices advected by the far field of the other vortices, § 4.2 shows that the self-advection is zero in the point-vortex limit and § 4.3 summarizes the QG
$^{+1}$ point vortex point-vortex equations of motion. Section 5 describes some of the properties of QG
$^{+1}$ point vortex point-vortex dynamics obtained by analysing the equations of motion. Section 6 briefly describes a few numerical simulations of the QG
$^{+1}$ point-vortex equations, highlighting some of the new features of QG
$^{+1}$ point vortices.
2. QG
$^{+1}$ theory
Following MSR1999, we begin with an inviscid, adiabatic, 3-D Boussinesq hydrostatic fluid on an 
$f$-plane (i.e. constant Coriolis parameter 
$f$) described by the primitive equations,
$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}$$
$$\begin{gather}\frac{{\rm D} \boldsymbol{u}_h}{{\rm D}t} + f \hat{\boldsymbol{z}} \times \boldsymbol{u}_h ={-}\boldsymbol{\nabla}_h \phi^T, \end{gather}$$
$$\begin{gather}\frac{{\rm D}\theta^T}{{\rm D}t} = 0, \end{gather}$$
$$\begin{gather}\partial_z\phi^T = g \theta^T/\theta_0, \end{gather}$$
$$\begin{gather}\frac{{\rm D}}{{\rm D} t} = \partial_t + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}, \end{gather}$$
where 
$\boldsymbol {u} = u\hat{\boldsymbol{x}}+ v\hat{\boldsymbol{y}} + w\hat{\boldsymbol{z}}$ is the 3-D fluid velocity, 
$\boldsymbol {u}_h = u\hat{\boldsymbol{x}}+ v\hat{\boldsymbol{y}}$ is the horizontal velocity, 
$\boldsymbol {\nabla }$ and 
$\boldsymbol {\nabla }_h$ are the 3-D and horizontal gradient operators, respectively, 
$\theta ^T$ is the total potential temperature and 
$\theta _0$ is a reference potential temperature. We follow MSR1999 and refer to 
$\phi ^T$ and 
$z$ as the total pressure and vertical height, although strictly speaking they are more properly a modified pressure and the geopotential height (Hoskins & Bretherton Reference Hoskins and Bretherton1972). Equations (2.1)–(2.5) represent, in order, incompressibility, conservation of horizontal momentum, conservation of energy, hydrostatic balance and the definition of the 3-D material derivative.
We restrict ourselves to flows with a constant Brunt–Väisälä frequency 
$N$. The pressure and potential temperature are separated into horizontally uniform reference fields 
$\phi ^{ref}$ and 
$\theta ^{ref}$, and departures from the reference fields
$$\begin{gather} \phi^{ref}= \tfrac{1}{2} N^2 z^2, \end{gather}$$
$$\begin{gather}g \theta^{ref}/\theta_0 = N^2 z, \end{gather}$$
$$\begin{gather}\phi^T = \phi + \phi^{ref}, \end{gather}$$
$$\begin{gather}\theta^T = \theta + \theta^{ref}. \end{gather}$$ The dynamics is non-dimensionalized as in MSR1999 in terms of a vertical distance 
$H$, a horizontal distance 
$L$ and a velocity 
$U$. Two non-dimensional parameters arise: the Burger number which is taken to be unity, 
$B = (NH/fL)^2=1$, and the Rossby number which is assumed to be small, 
$Ro= U/fL \ll 1$. These scalings give the non-dimensional primitive equations. Three scalar fields have explicit advective dynamics: the two components of the horizontal velocity vector and the potential temperature. The remaining variables, the pressure and vertical velocity, are determined diagnostically from incompressibility and hydrostatic balance, respectively. MSR1999 write the dynamics in terms of a 3-D vector potential 
$\boldsymbol {A}$, with 
$\boldsymbol {u} = \boldsymbol {\nabla }\times \boldsymbol {A}$ and 
$\theta = -\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {A}$. The curl relation for the velocity guarantees that the fluid is incompressible in three dimensions.
The relevant potential vorticity is the Ertel potential vorticity, 
$Q = (\hat{\boldsymbol{z}} + Ro \boldsymbol {\nabla }\times \boldsymbol {u}_h)\boldsymbol {\cdot }\boldsymbol {\nabla }\theta ^T$. The dynamically active disturbance potential vorticity 
$q$ is the departure from a uniform potential vorticity which, in this scaling, is unity, 
$Q = 1 + Ro \, q$. Following MSR1999, we refer to this disturbance Ertel potential vorticity 
$q$ as just the potential vorticity (PV).
Technically, QG appears at 
$O(\textit {Ro})$ in the full asymptotic expansion. As seen, for example, in the relation between 
$Q$ and 
$q$ above, the 
$O(1)$ term is uninteresting and is typically removed, following which the QG quantities are scaled by 
$\textit {Ro}$ and are rendered 
$O(\textit {Ro}^0)$. At next order, AG effects enter. QG and AG fields will be denoted, respectively, by the subscripts ‘0’ and ‘1.’ The asymptotic expansions then take the form, using the vector potential as an example, of
\begin{equation} \boldsymbol{A} = \boldsymbol{A}_0 + Ro \,\boldsymbol{A}_1 + O(Ro^2). \end{equation}The familiar QG relations between horizontal velocity, PV and streamfunction
\begin{equation} \boldsymbol{u} = \begin{pmatrix} -\partial_y\psi_0\\ \partial_x \psi_0 \\ 0 \end{pmatrix}, \quad q_0=\nabla^2\psi_0, \end{equation}
are rewritten in terms of the vector potential which is the solution to a Poisson equation 
$\nabla ^2 \boldsymbol {A}_0 = \boldsymbol {S}_0$ where 
$\boldsymbol {S}_0$ is the QG source,
\begin{equation} \boldsymbol{A}_0 = \begin{pmatrix} 0\\0\\-\psi_0 \end{pmatrix}, \quad \boldsymbol{S}_0 = \begin{pmatrix} 0\\0\\-q_0 \end{pmatrix}. \end{equation}
This structure carries over to higher order 
$\boldsymbol {u}_n=\boldsymbol {\nabla }\times \boldsymbol {A}_n$, 
$\nabla ^2 \boldsymbol {A}_n = \boldsymbol {S}_n$, where, unlike the QG terms, the higher order 
$\boldsymbol {A}_n$ 
$\boldsymbol {S}_n$ are generally non-zero in all three vector components.
The QG
$^{+1}$ model is a specific instance of an iterated balance model (Allen Reference Allen1993). The iteration procedure is based on treating one distinguished physical variable as exact and iteratively computing corrections to other variables at each instant of time. For example, Allen (Reference Allen1993) chose pressure as the distinguished variable. QG
$^{+1}$ chooses PV which makes it natural to seek balanced point-vortex solutions in QG
$^{+1}$ . Other balance models may also have point-vortex solutions but we defer that question to future work.
From a vorticity perspective, QG flow is completely described by the QG PV, 
$q_0(\boldsymbol {x},t)$. The QG streamfunction is obtained by solving a Poisson equation with 
$q_0$ as the source, which then determines the velocity and potential temperature. Conservation of QG PV under QG flow is given by
\begin{equation} \frac{{\rm D}_0 q_0}{{\rm D} t} = \partial_t q_0 + (\boldsymbol{u}_0\boldsymbol{\cdot}\boldsymbol{\nabla}) q_0 = 0, \end{equation}
where 
${\rm D}_0/{\rm D}t$ indicates the material derivative of fluid parcels advected with the QG velocity.
QG
$^{+1}$ dynamics is governed by PV conservation, where the full PV 
$q$ is advected by the approximate velocity
\begin{equation} \frac{{\rm D} q}{{\rm D}t} = \frac{\partial q}{\partial t} + ( (\boldsymbol{u}_0 + \textit{Ro}\, \boldsymbol{u}_1 + O(\textit{Ro}^2))\boldsymbol{\cdot}\boldsymbol{\nabla}) q = 0. \end{equation}
The PV can be partitioned into QG and AG contributions, 
$q = q_0 + \textit {Ro}\, q_1 + O(\textit {Ro}^2)$, however, as discussed in MSR1999, this partition has some subtleties. For our purposes, we will show that for point vortices, 
$q_1 = 0$ is consistent with QG
$^{+1}$ and proceed assuming 
$q = q_0$.
MSR1999 showed that the AG source 
$\boldsymbol {S}_1$ is an operator acting on the QG streamfunction 
$\psi _0$ plus the vector 
$(0,0,-q_1)$. With our assumption that 
$q_1 = 0$, the AG source becomes
\begin{equation} \boldsymbol{S}_1=\begin{pmatrix} -2 J(\partial_z\psi_0,\partial_y\psi_0)\\ 2 J(\partial_z\psi_0,\partial_x\psi_0) \\ (\nabla^2 \psi_0)\partial_{zz} \psi_0 - | \boldsymbol{\nabla} (\partial_z \psi_0)|^2 \end{pmatrix}, \end{equation}
where 
$J$ is the usual horizontal Jacobian operator, 
$J(\,f,g) = (\partial _x f) (\partial _y g) - (\partial _y f)(\partial _x g)$. Depending on the context, it is convenient to consider 
$\boldsymbol {S}_1$ as either a function of position 
$\boldsymbol {S}_1(\boldsymbol {x})$ obtained from considering 
$\psi _0(\boldsymbol {x})$ as an explicit function of position, or as a differential operator acting on a function 
$\psi _0$. Since 
$\psi _0$ appears quadratically in 
$\boldsymbol {S}_1$, we can, in the operator view, consider 
$\boldsymbol {S}_1$ to an operator acting on two potentially independent scalar functions 
$f$ and 
$g$,
\begin{equation} \boldsymbol{S}_1 =\begin{pmatrix} -2 J(\partial_z f,\partial_y g)\\ \hphantom{-}2 J(\partial_z f, \partial_x g) \\ (\nabla^2 f)\partial_{zz} g - (\boldsymbol{\nabla} \partial_z f)\boldsymbol{\cdot} (\boldsymbol{\nabla} \partial_z g) \end{pmatrix}. \end{equation}
In this perspective, 
$\boldsymbol {S}_1$ is a bilinear non-symmetric vector differential operator with the properties
$$\begin{gather} \boldsymbol{S}_1(c_1 f_1+c_2 f_2,g) = c_1\boldsymbol{S}_1(\,f_1,g)+c_2 \boldsymbol{S}_1(\,f_2,g), \end{gather}$$
$$\begin{gather}\boldsymbol{S}_1(\,f, c_1 g_1 + c_2 g_2) = c_1\boldsymbol{S}_1(\,f,g_1)+c_2 \boldsymbol{S}_1(\,f,g_2), \end{gather}$$
$$\begin{gather}\boldsymbol{S}_1(\,f,g) \neq \boldsymbol{S}_1(g,f), \end{gather}$$
where 
$f$ and 
$g$ are scalar functions of 
$\boldsymbol {x}$ and 
$c_1$ and 
$c_2$ are constants. The bilinearity of 
$\boldsymbol {S}_1$ plays a large role in the subsequent development of QG
$^{+1}$ vortex gas vortex-gas dynamics where 
$\psi _0$ is decomposed into a sum of contributions from individual vortices. Note that since the AG vector potential, through the AG source, depends quadratically on the QG streamfunction, it manifestly breaks the cyclone–anticyclone (anti)symmetry of QG,
$$\begin{gather} \boldsymbol{A}_0 \xrightarrow[q_0 \to -q_0]{} -\boldsymbol{A}_0, \end{gather}$$
$$\begin{gather}\boldsymbol{A}_1 \xrightarrow[q_0 \to -q_0]{} +\boldsymbol{A}_1, \end{gather}$$
and so 
$\boldsymbol {A} = \boldsymbol {A}_0 + \textit {Ro} \,\boldsymbol {A}_1$ is non-symmetric under sign changes of the PV.
MSR1999 consider a fluid with doubly periodic horizontal boundaries and solid vertical boundaries. It is usual to consider point-vortex dynamics in an infinite domain with boundary conditions requiring physical fields to decay to zero as one goes to infinity. Then, point-vortex dynamics in closed and periodic domains are obtained through the method of images. Here, we follow the point vortex viewpoint and consider QG
$^{+1}$ in a 3-D infinite domain and require that the velocity and potential temperature go to zero at infinity.
The QG
$^{+1}$ equations are completed by solvability constraints, which, in the infinite domain considered here, are
$$\begin{gather} \iiint q_0 \,{\rm d}^3\boldsymbol{x} = \lim_{z\to\infty}\iint \theta_0|^{{+}z}_{{-}z} \,{{\rm d}\kern0.06em x} \,{{\rm d} y}, \end{gather}$$
$$\begin{gather}\iiint q_1 \,{\rm }^3\boldsymbol{x} = \lim_{z\to\infty}\iint [\theta_1 + \partial_z \psi_0(\partial_{xx}\psi_0 + \partial_{yy}\psi_0)]|^{{+}z}_{{-}z}\, {{\rm d}\kern0.06em x} \,{{\rm d} y}, \end{gather}$$
where (2.22) guarantees the invertibility of the QG PV–streamfunction relation, and (2.23) guarantees the invertibility of the 
$z$-component of (2.15).
3. The vortex-gas approximation
We begin by assuming the fluid takes the form of a so-called ‘vortex gas’, where the PV consists of a collection of 
$N$ physically separated coherent vortices. Each vortex has a location, 
$\boldsymbol {x}_i(t)$, and a circulation 
$\varGamma _i$. We assume each coherent vortex has the same finite-sized axisymmetric shape in PV, 
$\hat {q}$. The only time dependence in the flow is in the location of the vortices. The number of vortices, their shape and their circulations are assumed to be constant. The PV field is then
\begin{equation} q(\boldsymbol{x}) = \sum_{i=1}^N \varGamma_i\hat{q}(\boldsymbol{x} - \boldsymbol{x}_i(t)), \end{equation}
where 
$\hat {q}$ is normalized to have unit circulation, 
$\int \hat {q} \, {\rm d}^3\boldsymbol {x} = 1$. Potential vorticity conservation, (2.14), then requires the vortex locations move with the local velocity,
\begin{equation} \frac{{\rm d} \boldsymbol{x}_i(t)}{{\rm d}t} = \boldsymbol{u}(\boldsymbol{x}_i) = \boldsymbol{u}_0(\boldsymbol{x}_i) + \textit{Ro}\, \boldsymbol{u}_1(\boldsymbol{x}_i) + O(\textit{Ro}^2), \end{equation}where we have used
\begin{equation} \frac{\partial \hat{q}(\boldsymbol{x} - \boldsymbol{x}_i(t))}{\partial t} ={-} \boldsymbol{\nabla}\hat{q}(\boldsymbol{x})|_{\boldsymbol{x} = \boldsymbol{x} - \boldsymbol{x}_i}\boldsymbol{\cdot} \frac{{\rm d} \boldsymbol{x}_i(t)}{{\rm d}t}. \end{equation}
Thus, while conservation of the full 
$q$ in QG
$^{+1}$ gives AG corrections to the motion of the vortices, the vortex itself is advected coherently.
The logic of QG
$^{+1}$ , reflecting the iterated nature of balance models generally, is that one begins with the QG PV 
$q_0$, uses the QG relations to obtain the QG streamfunction 
$\psi _0$, and then uses 
$\psi _0$ to obtain the AG fields. Following this logic for a vortex gas, the fields decompose into sums over shape fields, denoted by hats, which are identical for all vortices. Denoting the asymptotic order by the subscript 
$n=0$ for QG and 
$n=1$ for AG, and letting 
$F$ represent any of the fields 
$\boldsymbol {S}$, 
$\boldsymbol {A}$, 
$\boldsymbol {u}$ and 
$\theta$, vortex-gas fields take the form
\begin{gather} F_0(\boldsymbol{x}, t) = \sum_{i=1}^N \varGamma_i \hat{F}_0(\boldsymbol{x} - \boldsymbol{x}_i), \end{gather}
\begin{gather} F_1(\boldsymbol{x}, t) = \sum_{i, j =1}^N \varGamma_i \varGamma_j \hat{F}_1(\boldsymbol{x} - \boldsymbol{x}_i, \boldsymbol{x} - \boldsymbol{x}_j). \end{gather}
The QG shape fields 
$\hat {F}_0$ are single-vortex fields and depend on the location of a single vortex, while the AG shape fields 
$\hat {F}_1$ are pair-vortex fields and depend on the locations of a pair of vortices. The pair-vortex form of the AG fields is a consequence of both the vortex-gas approximation and the bilinear nature of the AG source operator. This vortex-pair interaction is fundamentally different from the vortex-pair dynamics of 2-D and QG point vortices. In 2-D and QG dynamics, all fields are single-vortex fields and the interaction is solely through the sum of the advection induced by each vortex separately. In QG
$^{+1}$ , vortex pairs interact in creating the AG source, 
$\boldsymbol {S}_1$, for the AG vector potential, 
$\boldsymbol {A}_1$ and the advecting velocity is the sum of velocities induced by vortex pairs. At next order in 
$\textit {Ro}$, three-vortex and four-vortex contributions would enter.
The AG fields can be rewritten in terms of single-vortex nonlinear contributions 
$\hat {F}_1^{s}$, and symmetrized pair-vortex nonlinear contributions 
$\hat {F}_1^p$
\begin{gather} \hat{F}_1^{s}(\boldsymbol{x} - \boldsymbol{x}_i) = \hat{F}_1(\boldsymbol{x} - \boldsymbol{x}_i, \boldsymbol{x} - \boldsymbol{x}_i), \end{gather}
\begin{gather} \hat{F}_1^{p}(\boldsymbol{x} - \boldsymbol{x}_i, \boldsymbol{x} - \boldsymbol{x}_j) = \hat{F}_1(\boldsymbol{x} - \boldsymbol{x}_i, \boldsymbol{x} - \boldsymbol{x}_j) + \hat{F}_1(\boldsymbol{x} - \boldsymbol{x}_j, \boldsymbol{x} - \boldsymbol{x}_i), \end{gather}resulting in
\begin{equation} F_1(\boldsymbol{x},t) = \sum_{i=1}^N \varGamma_i^2 \hat{F}_1^{s}(\boldsymbol{x} - \boldsymbol{x}_i) + \sum_{\substack{i=1\\j=i+1}}^N \varGamma_i \varGamma_j\hat{F}_1^{p}(\boldsymbol{x} - \boldsymbol{x}_i, \boldsymbol{x} - \boldsymbol{x}_j). \end{equation}
In all sums we use the convention that terms in a sum that contradict the bounds, such as the term in the second sum above with 
$i=N$ and 
$j=i+1 = N+1$, are zero.
Following the logic of a QG
$^{+1}$ vortex gas, all the shape fields are determined by the potential vorticity shape field 
$\hat {q}$ through solving Poisson equations and applying differential operators
$$\begin{gather} \hat{\boldsymbol{S}}_{0}(\boldsymbol{x}) = \begin{pmatrix} 0\\ 0\\ -\boldsymbol{\nabla}^2\hat{\psi}_0(\boldsymbol{x}) \end{pmatrix}=\begin{pmatrix} 0\\ 0\\ -\hat{q}(\boldsymbol{x}) \end{pmatrix}, \end{gather}$$
$$\begin{gather}\hat{\boldsymbol{S}}_1^s(\boldsymbol{x} - \boldsymbol{x}_i) = \boldsymbol{S}_1(\hat{\psi}_0(\boldsymbol{x} - \boldsymbol{x}_i),\hat{\psi}_0(\boldsymbol{x} - \boldsymbol{x}_i)), \end{gather}$$
$$\begin{gather}\hat{\boldsymbol{S}}_1^p(\boldsymbol{x} - \boldsymbol{x}_i, \boldsymbol{x} - \boldsymbol{x}_j) = \boldsymbol{S}_1(\hat{\psi}_0(\boldsymbol{x} - \boldsymbol{x}_i),\hat{\psi}_0(\boldsymbol{x} - \boldsymbol{x}_j)) + i \longleftrightarrow j, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}^2 \hat{\boldsymbol{A}}_n^{\{s,p\}} = \hat{\boldsymbol{S}}_n^{\{s,p\}}, \end{gather}$$
$$\begin{gather}\hat{\boldsymbol{u}}_n^{\{s,p\}} = \boldsymbol{\nabla}\times \hat{\boldsymbol{A}}_n^{\{s,p\}}, \end{gather}$$
where the operator 
$\boldsymbol {S}_1$ is given by (2.15).
Vortex-gas advection, (3.2), becomes
\begin{align} \frac{{\rm d} \boldsymbol{x}_i(t)}{{\rm d}t} &= \sum_{j=1}^N \varGamma_j \hat{\boldsymbol{u}}_0(\boldsymbol{x}_i- \boldsymbol{x}_j) \nonumber\\ &\quad + \textit{Ro} \left(\sum_{j=1}^N \varGamma_j^2 \hat{\boldsymbol{u}}_1^s(\boldsymbol{x}_i-\boldsymbol{x}_j) + \sum_{\substack{j=1\\k=j+1}}^N \varGamma_j \varGamma_k \hat{\boldsymbol{u}}_1^p(\boldsymbol{x}_i-\boldsymbol{x}_j, \boldsymbol{x}_i - \boldsymbol{x}_k)\right). \end{align}
The advecting velocities fall into two categories: far-field advection, where the single vortex 
$j$ and the pair of vortices 
$j, k$ are different than the advected vortex 
$i$, and self-advection, where either 
$i = j$ or 
$i = k$. It is well known that 2-D and QG vortices do not self-advect on the 
$f$-plane considered here, 
$\hat {\boldsymbol {u}}_0(0) = 0$ (vortices on the 
$\beta$-plane do, however, self-advect). Below, we will show that in addition there is no AG self-advection, 
$\hat {\boldsymbol {u}}_1^s(0) = \hat {\boldsymbol {u}}_1^p(0,\boldsymbol {x}) = \hat {\boldsymbol {u}}_1^p(\boldsymbol {x},0) = 0$.
4. QG$^{+1}$ point vortex point-vortex dynamics
4.1. Far-field advection
The far-field vector potentials can be determined by taking the vortex size to zero before performing the calculations, 
$\hat {q} = \delta (\boldsymbol {x})$. We assume for now that 
$q_1 = 0$, i.e. the QG point vortex carries all the PV. We will show below that this choice is consistent with the solvability constraint equation (2.23). To solve for single-vortex fields we can assume the vortex is at the origin.
As the QG vortex shape is a 
$\delta$-function, the streamfunction shape, solving 
$\nabla ^2\hat {\psi }_0 = \delta (\boldsymbol {x})$, is the Green's function, 
$\hat {\psi }_0 = G(\boldsymbol {x}) = -1/4 {\rm \pi}|\boldsymbol {x}|$, and the QG vector potential shape is 
$\hat {\boldsymbol {A}}_0 = (0,0, -\hat {\psi }_0)$. The single-vortex AG Poisson source 
$\hat {\boldsymbol {S}}_1^s$ from a vortex at the origin can be directly calculated, leading to
\begin{equation} \hat{\boldsymbol{S}}_1^s = \frac{1}{8 {\rm \pi}^2 |\boldsymbol{x}|^8} \begin{pmatrix} 3 x z \\ 3 y z \\ -(x^2+y^2+4 z^2)/2 \end{pmatrix}. \end{equation} We solve the Poisson problem by first seeking a solution that matches the denominator of the source. We note that 
$\nabla ^2 |\boldsymbol {x}|^{-m} \sim |\boldsymbol {x}|^{-m-2}$. Since the denominator of 
$\hat {\boldsymbol {S}}_1^s$ is proportional to 
$|\boldsymbol {x}|^{8}$, we expect the denominator of 
$\hat {\boldsymbol {A}}_1^s$ to be proportional to 
$|\boldsymbol {x}|^{6}$. Dimensionally, 
$\hat {\boldsymbol {S}}_1^s \sim 1/L^6$, so 
$\hat {\boldsymbol {A}}_1^s \sim 1/L^4$. This means the numerator of 
$\hat {\boldsymbol {A}}_1^s$ must be quadratic in 
$(x, y, z)$. We thus set the numerator of each component of 
$\hat {\boldsymbol {A}}_1^s$ to a general quadratic 
$a x^2 + b y^2 + c z^2 + {{\rm d}\kern0.06em x} y + e x z + f y z$, and take the Laplacian of this general 
$\hat {\boldsymbol {A}}_s^2$. Requiring that 
$\hat {\boldsymbol {A}}_1^s$ solves the Poisson equation determines the values of the constants in the quadratic. The result is
\begin{equation} \hat{\boldsymbol{A}}_1^s(\boldsymbol{x}) = \frac{z}{16 {\rm \pi}^2 |\boldsymbol{x}|^6} \begin{pmatrix} x\\ y\\ -z/2 \end{pmatrix}. \end{equation}The next step is to solve the Poisson equation for the far-field pair-vortex potential. We notice that the components of the far-field single-vortex potential (4.2) can each be written as the product of two functions
\begin{equation} \hat{\boldsymbol{A}}_1^s = \begin{pmatrix} h_x(\boldsymbol{x}) h_z(\boldsymbol{x}) \\ h_y(\boldsymbol{x}) h_z(\boldsymbol{x}) \\ - h_z(\boldsymbol{x})^2/2 \end{pmatrix}, \end{equation}
where 
$\boldsymbol {h} = \boldsymbol {x}/4 {\rm \pi}|\boldsymbol {x}|^3 = \boldsymbol {\nabla } G(\boldsymbol {x})$ and its subscript represents the corresponding vector component. We guess that 
$\hat {\boldsymbol {A}}_1^p$ takes a similar symmetrized form
\begin{equation} \hat{\boldsymbol{A}}_1^p(\boldsymbol{x}-\boldsymbol{x}_i, \boldsymbol{x}-\boldsymbol{x}_j) =\begin{pmatrix} h_x(\boldsymbol{x}-\boldsymbol{x}_i) h_z(\boldsymbol{x}-\boldsymbol{x}_j) + h_x(\boldsymbol{x}-\boldsymbol{x}_j) h_z(\boldsymbol{x}-\boldsymbol{x}_i) \\ h_y(\boldsymbol{x}-\boldsymbol{x}_i) h_z(\boldsymbol{x}-\boldsymbol{x}_j) + h_y(\boldsymbol{x}-\boldsymbol{x}_j) h_z(\boldsymbol{x}-\boldsymbol{x}_i) \\ -h_z(\boldsymbol{x}-\boldsymbol{x}_i) h_z(\boldsymbol{x}-\boldsymbol{x}_j) \end{pmatrix}.\end{equation}
Direct calculation of 
$\boldsymbol {\nabla }^2 \hat {\boldsymbol {A}}_1^p$ verifies that this is indeed the solution. These solutions manifestly satisfy the boundary conditions, going to zero at infinity.
The right-hand side of the integral constraints involve integrals evaluated on the vertical boundaries, which we take here to be at 
$z = \pm \infty$. We assume all vortices are at finite heights, and then the integral constraints only depend on the far fields. The QG integral constraint equation (2.22) can be integrated analytically and verified. We have not found an analytic form for the integral on the right-hand side of (2.23). However, the scaling behaviour of the fields lets us deduce that the integral is zero. Since 
$\hat {\boldsymbol {A}}_1$ decays as 
$|\boldsymbol {x}|^{-4}$, 
$\hat \theta _1$ decays as 
$|\boldsymbol {x}|^{-5}$. The value of 
$\hat \psi _0$ decays as 
$|\boldsymbol {x}|^{-1}$, so the second term inside the integral also decays as 
$|\boldsymbol {x}|^{-5}$. The right-hand side thus goes to zero as 
$z\to \pm \infty$, indicating that the integral of 
$q_1$ is zero, consistent with 
$q_1 = 0$.
4.2. Self-advection
Point-vortex self-advection is the limit of the finite-sized coherent vortex self-advection as the vortex size goes to zero. This requires some consideration of a finite vortex. Here, we consider all vortices to have the same time-independent, spherically symmetric shape with size 
$r_0$.
In 2-D and QG point-vortex dynamics, the self-advection velocity of vortex 
$i$ is the velocity at 
$\boldsymbol {x}$ induced by a single finite-sized coherent vortex in the ordered limit, 
$\lim _{r_0\to 0}\lim _{\boldsymbol {x}\to \boldsymbol {x}_i}$. In QG
$^{+1}$ point vortex point-vortex dynamics, there are three self-advection velocities: 
$\hat {\boldsymbol {u}}_0(0)$ the QG velocity at a vortex centre induced by its own vorticity distribution; 
$\hat {\boldsymbol {u}}_1^s(0)$, the AG velocity at a vortex centre induced by the nonlinear interaction between the vortex and itself; and 
$\hat {\boldsymbol {u}}_1^p(0,\Delta \boldsymbol {x})$, the AG velocity at a vortex centre induced by the nonlinear interaction between the advected vortex and a different vortex separated by 
$\Delta \boldsymbol {x}$.
The consequences of the assumed vortex shape symmetry are leveraged using the 3-D infinite-domain Green's function, written as
\begin{equation} G(\boldsymbol{x},\boldsymbol{x}') = \frac{-1}{4{\rm \pi} |\boldsymbol{x} - \boldsymbol{x}'|}, \end{equation}and the Green's function solution of the Poisson equation,
\begin{equation} \boldsymbol{\nabla}^2 \boldsymbol{A}_n(\boldsymbol{x}) = \boldsymbol{S}_n(\boldsymbol{x}) \,\Longleftrightarrow \,\boldsymbol{A}_n(\boldsymbol{x}) = \int G(\boldsymbol{x}, \boldsymbol{x}') \boldsymbol{S}_n(\boldsymbol{x}') \,{\rm d}V'. \end{equation}The velocity is then
\begin{equation} \boldsymbol{u}_n(\boldsymbol{x}) = \boldsymbol{\nabla} \times \boldsymbol{A}_n(\boldsymbol{x}) = \int \boldsymbol{\nabla}_{\boldsymbol{x}} G(\boldsymbol{x},\boldsymbol{x}') \times \boldsymbol{S}_n(\boldsymbol{x}') \,{\rm d}^3\boldsymbol{x}', \end{equation}
where the subscript on the gradient operator indicates derivatives are with respect to 
$\boldsymbol {x}$.
We first recall the QG point-vortex self-advection, 
$\hat {\boldsymbol {u}}_0(0)$, which can be calculated by considering a single unit-circulation vortex at the origin. Then
\begin{equation} \hat{\boldsymbol{u}}_0(0) = \lim_{r_0\to0}\int [\boldsymbol{\nabla}_{\boldsymbol{x}} G(\boldsymbol{x},\boldsymbol{x}')]_{\boldsymbol{x} = 0} \times \hat{\boldsymbol{S}}_{0}(\boldsymbol{x}') \,{\rm d}^3\boldsymbol{x}'. \end{equation}
The symmetries in the vortex shape and the Green's function render the argument of the integral in (4.8) odd in one direction and thus 
$\hat {\boldsymbol {u}}_0(0) = 0$, there is no QG point vortex self-advection.
The AG single-vortex self-advection 
$\hat {\boldsymbol {u}}_1^s(0)$ is also zero by symmetry but the argument is a bit more complicated. The relevant source is 
$\hat {\boldsymbol {S}}_1^s(\boldsymbol {x}')$. The derivatives in (2.15), allow one to determine the symmetry of the components of 
$\hat {\boldsymbol {S}}_1^s(\boldsymbol {x}')$ in the different directions. One can show that each component has at least one direction where the integrand is odd and integrates to zero. Thus, 
$\hat {\boldsymbol {u}}_1^s(0) = 0$ and there is no AG single-vortex self-advection.
Physically, we know that single coherent vortices on the 
$f$-plane do not self-advect and so it is satisfying that symmetry-based arguments show no single-vortex self-advection. The pair self-advection is different. Since a pair of QG vortices do advect each other, there is no physical reason why the AG pair self-advection must be zero. A finite AG pair self-advection would just provide a small correction to the QG mutual advection and is not a priori ruled out on physical grounds.
Because the pair self-advection relies on a second vortex at location 
$\boldsymbol {x}_2 \ne 0$, one cannot use symmetry arguments to determine its value. Instead we explicitly consider a specific finite vortex shape. One common finite vortex shape is a Gaussian monopole
$$\begin{gather} \hat{q}^{r_0} = \frac{{\rm e}^{{-}r^2/2 r_0^2}}{(2{\rm \pi} r_0^2)^{3/2}}, \end{gather}$$
$$\begin{gather}\hat{\psi}_0^{r_0} = G(r) \,\text{Erf}(r/(\sqrt{2} r_0)), \end{gather}$$
where 
$r = |\boldsymbol {x}|$, 
$\text {Erf}(s) = ({2}/{\sqrt {{\rm \pi} }})\int _0^s {\rm e}^{-t^2}\,{\rm d}t$ is the error function, and the superscript 
$r_0$ denotes a finite vortex.
To calculate the AG pair self-advection of a vortex, we can place that vortex at the origin, 
$\boldsymbol {x}_1 = 0$. We are interested in the velocity near the vortex centre, the origin. In the vortex-gas approximation, the second vortex is far away 
$|\boldsymbol {x}_2| \gg r_0$, and we can use its far-field streamfunction 
$\hat {\psi }_0 = G(\boldsymbol {x} - \boldsymbol {x}_2)$ in calculating 
$\hat {\boldsymbol {S}}_1^p$
\begin{equation} \hat{\boldsymbol{S}}_1^p = \boldsymbol{S}_1(\hat{\psi}_0^{r_0}(\boldsymbol{x}), \hat{\psi}_0(\boldsymbol{x} - \boldsymbol{x}_2)) + \boldsymbol{S}_1( \hat{\psi}_0(\boldsymbol{x} - \boldsymbol{x}_2),\hat{\psi}_0^{r_0}(\boldsymbol{x})). \end{equation}
Expanding 
$\hat {\boldsymbol {S}}_1^p$ in a Taylor series near the origin, 
$r \ll r_0$, gives 
$\hat {\boldsymbol {S}}_1^p\sim O(r^0)$. This then implies that 
$\hat {\boldsymbol {A}}_1^p \sim O(r^2)$ and 
$\hat {\boldsymbol {u}}_1^p \sim O(r)$ near the origin. Thus 
$\hat {\boldsymbol {u}}_1^p \to 0$ as 
$r \to 0$ for finite Gaussian vortices and remains zero in the point vortex limit 
$r_0\to 0$. Thus, at least for the point-vortex limit of Gaussian vortices, the AG pair self-advection is zero.
4.3. Summary equations
As shown above, QG
$^{+1}$ point vortices, like 2-D and QG point vortices, have no self-advection. Thus, their dynamics is completely determined by their far-field advection. Putting together the pieces described above gives the equations of motion for a set of 
$N$ QG
$^{+1}$ point vortices
$$\begin{gather} \frac{{\rm d} \boldsymbol{x}_i}{{\rm d}t} = \boldsymbol{u}(\boldsymbol{x}_i),\notag\\ \boldsymbol{u}(\boldsymbol{x}_i) = \sum_{j=1}^N{}^{'} ( \varGamma_j\hat{\boldsymbol{u}}_0(\boldsymbol{x}_i - \boldsymbol{x}_j) + \textit{Ro}\, \varGamma_j^2 \hat{\boldsymbol{u}}_1^s(\boldsymbol{x}_i - \boldsymbol{x}_j)) \end{gather}$$
$$\begin{gather}\quad + \,\textit{Ro} \sum_{\substack{j=1\\k=j+1}}^N{}^{'} \varGamma_j \varGamma_k \hat{\boldsymbol{u}}_1^p(\boldsymbol{x}_i - \boldsymbol{x}_j, \boldsymbol{x}_i - \boldsymbol{x}_k), \end{gather}$$
where 
$\sum {}^{'}$ denotes skipping terms with 
$j = i$ and 
$k = i$. The single-vortex and pair-vortex velocities are
$$\begin{gather} \hat{\boldsymbol{u}}_0(\boldsymbol{x}) = \frac{1}{4 {\rm \pi}|\boldsymbol{x}|^3} \begin{pmatrix} -y\\ x\\ \,0 \end{pmatrix}, \end{gather}$$
$$\begin{gather}\hat{\boldsymbol{u}}_1^s(\boldsymbol{x}) = \frac{x^2+y^2-8z^2}{16 {\rm \pi}^2 |\boldsymbol{x}|^8} \begin{pmatrix} -y\\ x\\ \,0 \end{pmatrix}, \end{gather}$$
$$\begin{gather}\hat{\boldsymbol{u}}_1^p(\boldsymbol{x}_1,\boldsymbol{x}_2) = \frac{\tilde{\boldsymbol{u}}}{16 {\rm \pi}^2 |\boldsymbol{x}_1|^5 |\boldsymbol{x}_2|^5}, \end{gather}$$where
\begin{align} \tilde{\boldsymbol{u}} =\begin{pmatrix} 3 |\boldsymbol{x}_1|^2 (y_1 z_2^2 + 2 y_2 z_1 z_2) + 3|\boldsymbol{x}_2|^2(y_2 z_1^2 + 2 y_1 z_1 z_2) - |\boldsymbol{x}_1|^2 |\boldsymbol{x}_2|^2 (y_1+y_2) \\ |\boldsymbol{x}_1|^2 |\boldsymbol{x}_2|^2 (x_1+x_2)- 3 |\boldsymbol{x}_1|^2 (x_1 z_2^2 + 2 x_2 z_1 z_2) - 3|\boldsymbol{x}_2|^2(x_2 z_1^2 + 2 x_1 z_1 z_2)\\ 3 (x_2 y_1 - x_1 y_2)( |\boldsymbol{x}_2|^2 z_1 - |\boldsymbol{x}_1|^2 z_2). \end{pmatrix}. \end{align}We note that, as is the case for 2-D and QG point vortices, it is straightforward to study passive scalar transport by considering passive tracers as vortices with zero circulation.
Asymptotic expansions break down if the higher-order terms become larger than lower-order terms. Here, asymptotic consistency requires 
$|\hat {\boldsymbol {u}}_1^s| /|\hat {\boldsymbol {u}}_0| \sim O(1)$ and 
$|\hat {\boldsymbol {u}}_1^p| /|\hat {\boldsymbol {u}}_0| \sim O(1)$, which requires 
$\varGamma /4{\rm \pi} |\boldsymbol {x}|^3 \sim O(1)$. For single-vortex advection, 
$\boldsymbol {x}$ is the distance from the advected particle, either a vortex or passive, to the advecting vortex and 
$\varGamma$ is the circulation of the advecting vortex. For pair-vortex advection, 
$\boldsymbol {x}$ is the distance to the nearest member of the advecting pair and 
$\varGamma$ is its circulation. For a collection of point vortices with similar circulations, 
$\varGamma _i \sim O(\varGamma )$, there is thus a minimum distance between QG
$^{+1}$ vortices and between passives and vortices, an ‘asymptotic horizon’, if you will, 
$r_a \sim (\varGamma /4{\rm \pi} )^{1/3}$. The asymptotic horizon is a boundary in phase space and inside the boundary the asymptotic expansion becomes misordered. Systems where all vortices remain further apart than 
$r_a$ along their entire trajectories are outside the horizon and are asymptotically valid QG
$^{+1}$ solutions. Similarly, a passive particle which remains further than 
$r_a$ from all vortices is an asymptotically valid passive trajectory. There is no restriction on how close passives can approach each other as they do not induce any advective velocity. For systems with a wide range of 
$\varGamma$ values, the asymptotic horizon becomes more complex.
5. Properties of QG$^{+1}$ point vortex point-vortex dynamics
5.1. Vortex pairs
The qualitative behaviour of a system composed of only two vortices, 
$N=2$ can be understood from the dynamical equation, (4.12)–(4.17). We begin by noting that for two vortices, it is not possible to have 
$i$, 
$j$, 
$k$ which satisfy 
$i\ne j$, 
$i\ne k$ and 
$k>j$. As a result, the second sum in (4.13) has no terms and the pair-vortex velocity is zero.
Both 2-D and QG vortex pairs co-rotate in horizontal planes. Since 
$\hat {\boldsymbol {u}}_0$ and 
$\hat {\boldsymbol {u}}_1^s$ have zero vertical velocity and 
$\hat {\boldsymbol {u}}_1^p=0$, QG
$^{+1}$ vortex pairs also move in horizontal planes with no vertical velocity. In addition 
$\hat {\boldsymbol {u}}_1^s$ shares with 
$\hat {\boldsymbol {u}}_0$ the property that it is proportional to 
$(-y, x, 0)$ and is thus directed perpendicular to the line connecting the vortices. The AG velocities induced at vortices 
$1$ and 
$2$ are, as in two dimensions and QG, in opposite directions. Thus, QG
$^{+1}$ vortex pairs also co-rotate.
The AG correction can either speed up or slow down the QG rotation. The magnitude of the AG velocity scales as 
$\varGamma ^2$ and so is positive for both same-sign and opposite-sign vortex pairs. Furthermore, the numerator of 
$\hat {\boldsymbol {u}}_1^s$ changes sign on the surfaces 
$z = \pm \sqrt {(x^2+y^2)/8}$. As QG
$^{+1}$ vortex pairs rotate with constant vertical and horizontal separations, the sign of the numerator is constant. Thus, whether QG
$^{+1}$ vortex pairs rotate faster or slower than QG vortex pairs depends on both the signs of the circulations and the specific vector vortex separation. The difference in rotation speeds with the sign of the vortex circulation is one example of QG
$^{+1}$ breaking the cyclone–anticyclone anti-symmetry.
The relative magnitudes of the AG and QG vortex speeds, 
$|\boldsymbol {u}_1^s|/|\boldsymbol {u}_0|$ scales as 
$\textit {Ro}\varGamma /{|\boldsymbol {x}_1-\boldsymbol {x}_2|^3}$. Thus, as the vortex separation grows, the AG velocity decays much faster than the QG velocity.
Opposite-sign 2-D and QG vortices with 
$\varGamma _1 = - \varGamma _2$, called hetons in the oceanographic context, propagate along straight lines. Due to the dependence of the AG velocity on 
$\varGamma ^2$, the two vortices no longer have the same horizontal velocity and they travel in circles, similar to opposite-sign unequal magnitude vortices in two dimensions and QG, with a radius of curvature that scales as 
$1/\textit {Ro}$ (figure 1). The heton curvature breaks the cyclone–anticyclone anti-symmetry. If the symmetry held, switching the signs of the circulations of the heton component vortices would cause the heton to rotate in the opposite direction. However, as the curvature is caused only by the AG contribution, the direction of curvature is the same regardless of which component vortex is positive and which is negative,
Figure 1. QG (dashed) and QG
$^{+1}$ (solid) heton trajectories viewed from above for initial vortex positions 
$\pm (1, 0, 1)$ with 
$\varGamma = \mp 4{\rm \pi}$ and 
$\textit {Ro} = 0.2$. The height and separations of the vortices comprising QG and QG
$^{+1}$ hetons are constant. The main figure shows the horizontal projection of the trajectories for 
$0 \le t \le 500$ and the inset shows the projection for the QG
$^{+1}$ heton for 
$0 \le t \le 2000$.
5.2. Solid horizontal boundary
In the oceanographic context, idealized vortex dynamics is typically studied by placing a solid boundary, representing the ocean surface, at 
$z=0$, where the vertical velocity is zero. This boundary condition can be satisfied for QG
$^{+1}$ point vortices with arbitrary 
$N$ using the method of images. Each physical vortex below the surface is supplemented by an image vortex at the same horizontal location and equally spaced above the surface, i.e. vortex 
$i$ at 
$\boldsymbol {x}_i = (x,y,z)$ has an image vortex at 
$\boldsymbol {x}_{i_{image}} = (x,y,-z)$. Then it can be seen from the equations of motion that the pair-vortex vertical velocity induced on the surface by a vortex and its image is zero. Furthermore, the pair-vortex vertical velocity at the surface induced by vortex 
$i$ and the image of vortex 
$j$ is non-zero but is cancelled by the vertical velocity induced by the image of vortex 
$i$ and vortex 
$j$. To keep the images vertically aligned with their physical vortices as the vortices evolve, their horizontal velocities must be equal, requiring 
$\varGamma _{i_{image}}=\varGamma _i$. Due to the 
$1/|\boldsymbol {x}|^3$ scaling of the ratio of the AG and QG velocities, the AG contribution from image vortices decays rapidly as vortex depth increases. The asymptotic horizon provides a minimum depth for QG
$^{+1}$ point vortices with a solid surface.
While image vortices keep the vertical velocity at the boundary zero, they do, in general, impact the vertical velocity below the boundary. Two physical vortices corresponds to a system with 
$N=4$, which will, in general, have vertical vortex advection. However, due to the symmetry of the images, the vertical advection of two physical vortices with a solid horizontal boundary is zero. Each vortex, vortex 1, say, has a vertical velocity from three pair-vortex interactions, 
$2$ and 
$2_{image}$, 
$2$ and 
$1_{image}$ and 
$1_{image}$ and 
$2_{image}$. Inspection of the dynamical equations shows that the advecting vertical velocity for each of these pairs is zero. Thus a solid horizontal boundary does not induce vertical vortex motion in a system of two physical vortices.
6. Interesting QG$^{+1}$ point vortex point-vortex numerical solutions
In both two dimensions and QG, there is a large literature of point-vortex and passive tracer solutions. A detailed study of QG
$^{+1}$ trajectories is beyond the scope of this work. Here, we choose a small number of configurations and present numerical simulations that highlight some of the new features of QG
$^{+1}$ point-vortex dynamics. All trajectories presented here are outside the asymptotic horizon. All QG
$^{+1}$ simulations are performed with 
$\textit {Ro} = 0.2$.
6.1. Passive particle motion induced by a same-sign vortex pair
As discussed above, QG
$^{+1}$ same-sign vortex pairs have no vertical vortex advection and co-rotate in horizontal planes. Passive tracer trajectories in the field of a QG
$^{+1}$ same-sign pair can experience vertical advection if the vortices are tilted, while if the vortices are vertically aligned, the vertical advection is zero. Point vortices, being non-diffusive, technically stir the fluid; mixing requires diffusion. Point-vortex stirring does, however, give a good indication of the mixing that would ensue in a vortex gas with small diffusion. A detailed study of stirring would require investigating the complex 3-D nature of passive particle trajectories. Here, we focus on the simpler vertical excursion of passive particles, the difference between the maximum and minimum height of a particle along its trajectory. This simpler measure already shows interesting behaviour.
We consider a tilted same-sign vortex pair with equal circulations 
$\varGamma _1 = \varGamma _2 = 4{\rm \pi}$ at locations 
$\boldsymbol {x} = \pm (1.5, 0, 1.5)$ and with 
$\textit {Ro} = 0.2$. The period 
$\tau$ of QG co-rotation is 
$\tau _{QG} \approx 240$. For this configuration, the AG correction slows the vortices down and the QG
$^{+1}$ co-rotation period is 
$\tau _{QG^{+1}}\approx 242$.
Figures 2–4 show trajectories from three different passive particle initial conditions chosen to display a variety of behaviours. The horizontal motion (panel a) looks quasiperiodic rather than chaotic. In a frame co-rotating with the vortices the motion becomes simpler and appears periodic (panel b). The initial condition of the passive particle in figure 4 is particularly interesting as it displays large slow vertical excursions. For this third initial condition, the upward (downward) moving branches of the trajectory occur with 
$y_{\textit {Co-Rotating}} < 0$ (
$y_{\textit {Co-Rotating}} > 0$), and the height extrema occur as the trajectory crosses the 
$y_{\textit {Co-Rotating}}=0$-axis.
Figure 2. Vortex and passive particle trajectories for a pair of same-sign vortices starting at 
$\pm (1.5,0,1.5)$ with 
$\varGamma _1 = \varGamma _2 = 4{\rm \pi}$, a passive particle starting at 
$(-3.0, 1.0, -1.8)$, 
$\textit {Ro} = 0.2$ and 
$0 \le t \le 1000$. (a) Horizontal projection of vortex trajectories (black line), passive trajectory (blue line), initial vortex locations (black circles) and initial passive location (blue circle). (b) Projection of the trajectories and initial positions onto horizontal coordinates co-rotating with the vortices. (c) Time series of the vortex and passive heights, 
$z(t)$.
Figure 3. Same as figure 2 for a passive particle starting at 
$(4.0, 0.1, 1.0)$.
Figure 4. Same as figure 2 for a passive particle starting at 
$(0.2, 2.7, 0.2)$ but with the time span in panel (a) being 
$0 \le t \le 10\,000$ and the time span in panels (b) and (c) being 
$0 \le t \le 100\,000$.
We explore these large vertical excursions by looking at the maximum vertical excursion of passive particles starting on the midplane of the vortices 
$z=0$ (figure 5). The symmetry of the figure reflects the symmetry of initial vortex positions and the symmetry of QG
$^{+1}$ dynamics. One sees that there is an 
$O(1)$-sized region of midplane passive initial conditions with 
$O(1)$ vertical excursions larger than the vertical vortex separation.
Figure 5. Vertical excursion, 
$\text {Max}(z(t)) - \text {Min}(z(t))$ of passive particles as a function of their initial position 
$(x,y,0)$. The black circles indicate the initial horizontal position of the vortices at 
$\pm (1.5, 0, 1.5)$, with 
$\varGamma _1 = \varGamma _2 = 4{\rm \pi}$ and 
$\textit {Ro} = 0.2$. The black diamond is the location of the line of passive initial conditions shown in figure 6.
We next look at the vertical excursion of passives with initial conditions along a vertical line, 
$\boldsymbol {x}(0) = (0, 3.82, z_0)$ going through q point near the maximum of figure 5 and denoted by a black diamond in that figure. The heights (figure 6) show that the region of large vertical excursions extends vertically across the entire vortex separation. This suggests that there are 3-D 
$O(1)$-sized lobes that are stirred in the vertical. Trajectories starting on this line appear to fall into two categories: those with 
$O(1)$ vertical excursions on long time scales, and those with very small vertical excursions and no long time scales.
Figure 6. Height vs time of passive particles in the field of a same-sign vortex pair. The vortices are initially at 
$\pm (1.5,0,1.5)$, have 
$\varGamma _1 = \varGamma _2=4 {\rm \pi}$, and 
$\textit {Ro} = 0.2$. The initial conditions of the passive particles are equally spaced in 
$z$ along a vertical line with 
$(x,y) = (0,3.82)$, indicated by the black diamond in figure 5.
6.2. Three same-sign vortices
Here we investigate the motion of a system composed of three same-sign vortices. We make no attempt to fully explore the large parameter space and instead focus on a few configurations which have interesting behaviour. The three vortices have identical circulations, 
$\varGamma _i = 4{\rm \pi}$ and are equally distributed in the vertical coordinate. The upper and lower vortices’ initial positions lie along a line in the 
$x-z$ plane, at positions 
$\pm (1.5, 0, 1.5)$. The middle vortex initially lies on the midplane, offset in the 
$y$-direction, at 
$(0, y_0,0)$.
We first consider an initial condition where all three vortices lie along a line, i.e. 
$y_0 = 0$. For this configuration, both QG and QG
$^{+1}$ point vortex point-vortex dynamics are relatively simple due to the symmetry of the initial condition. The induced velocities of the upper and lower vortices are horizontal and these vortices co-rotate. The middle vortex is stationary. The QG period of rotation for the upper and lower vortices is 
$\tau _{QG}=48.0$ and for QG
$^{+1}$ vortices with 
$\textit {Ro} = 0.2$, 
$\tau _{QG^{+1}} = 52.3$.
Next consider an initial condition with the middle vortex slightly offset from the origin 
$y_0=0.1$. The QG trajectories appear quasiperiodic (figure 7). The horizontal projection of the QG
$^{+1}$ trajectories also appear quasiperiodic but vortices now also display small vertical oscillations with periods similar to the horizontal rotation time scale (figure 8).
Figure 7. Trajectories of a system of three QG point vortices, 
$\textit {Ro} = 0.0$, initially at 
$\pm (1.5, 0, 1.5)$ and 
$(0,y_0,0)$ with 
$y_0=0.1$. The circulations are equal, 
$\varGamma _i = 4{\rm \pi}$. (a) Horizontal projection of the trajectory of the upper vortex for 
$0 \le t \le 1500$. The black circle indicates the initial position. (b) Same as panel (a) for the middle vortex. (c) Same as panel (a) for the lower vortex.
Figure 8. Trajectories of a system of three QG
$^{+1}$ point vortices with 
$\textit {Ro} = 0.2$, initially at 
$\pm (1.5, 0.1, 1.5)$ and 
$(0,y_0,0)$ with 
$y_0 = 0.1$. The circulations are equal, 
$\varGamma _i = 4{\rm \pi}$. (a) Horizontal projection of the trajectory of the upper vortex for 
$0 \le t \le 1500$. The black circle indicates the initial position. (b) Same as panel (a) for the middle vortex. The horizontal projection of the lower vortex trajectory is qualitatively similar to panel (a). (c) Height vs time of the three vortices.
Now consider a system where the middle vortex initially has a significantly larger horizontal offset, 
$y_0 = 3$ (figure 9). The QG dynamics still looks quasiperiodic. The QG
$^{+1}$ trajectories (figure 10) for this initial condition display similar horizontal trajectories as in the QG case for the time shown, a few horizontal rotation times. The notable new feature is the large 
$O(1)$ vertical oscillation which causes the vortices to exchange their vertical positions. The vertical motion appears periodic with a period 
$O(100)$ times longer than the horizontal rotation time scale. As the vortex heights slowly change, the shorter time scale horizontal motions change their character slowly but remain qualitatively similar.
Figure 9. Same as figure 7 (three QG point vortices) but with 
$y_0=3$.
Figure 10. Same as figure 8 (three QG
$^{+1}$ point vortices) but with 
$y_0 = 3$.
The vertical oscillations have an interesting structure as a function of the initial condition offset 
$y_0$ (figure 11). There is a region of 
$y_0$ with large vertical oscillations similar to those seen above for 
$y_0=3$. Note that this region is similar in location to the region found in the previous section where passive particles in the field of two vortices have large vertical excursions, but here the boundary of the region appears sharper. The amplitude of the vertical oscillation of the three vortices are qualitatively similar but with small quantitative differences. The periods of the oscillation of the three vortices are identical at each 
$y_0$.
Figure 11. Amplitude (a) and period (b) of the vertical oscillation of the vortex initially at 
$z=0$ in a system of three QG
$^{+1}$ point vortices initially at 
$\pm (1.5, 0, 1.5)$ and 
$(0, y_0, 0)$ as a function of the offset 
$y_0$. The vortices have 
$\varGamma _i = 4{\rm \pi}$ and 
$\textit {Ro} = 0.2$. In panel (a), the solid line is the amplitude of the initially upper and lower vortices and the dashed line is the amplitude of the initially middle vortex. The periods of the three vortices are identical.
6.3. chaotic motion
The above trajectories look quasiperiodic. Their Lyapunov exponents appear to converge to zero as the integration time grows. They are most likely not chaotic. However, like 2-D and QG point vortices, QG
$^{+1}$ point vortices do have chaotic motion. Four-vortex configurations can be found with chaotic trajectories with non-zero Lyapunov exponents (figure 12). The horizontal trajectories of the chaotic QG
$^{+1}$ vortices are qualitatively similar to that of their corresponding chaotic QG vortices. For this chaotic trajectory, the QG
$^{+1}$ vortices have 
$O(1)$ irregular vertical motion on long time scales.
Figure 12. Trajectories of a system of four vortices with 
$\varGamma _i = 4{\rm \pi}$ with positions chosen at random. (a) Horizontal motion of one vortex under QG dynamics, 
$\textit {Ro} = 0$, for 
$0 \le t \le 1500$. (b) Horizontal motion of the same vortex under QG
$^{+1}$ dynamics starting from the same initial condition with 
$\textit {Ro} = 0.2$. (c) Vertical motion of this QG
$^{+1}$ vortex.
7. Discussion
We have developed a new class of point-vortex dynamics, 3-D AG balanced point vortices. The solutions are found as point-vortex idealizations of vortex-gas solutions of a specific set of asymptotic balance equations, the QG
$^{+1}$ equations. The inclusion of AG dynamics in rapidly rotating, stably stratified flows brings both qualitative and quantitative differences. Horizontal velocities are perturbed by small, order-Rossby-number, corrections to QG horizontal flow. Vertical velocities are also small but as they are perturbations from zero, they produce qualitatively new phenomena.
We have seen that some configurations of QG
$^{+1}$ point vortices, and of passive particles in the field of QG
$^{+1}$ point vortices, have large 
$O(1)$ vertical transport. As the vertical velocity is small, large vertical transport requires correspondingly long times. We have seen that the regions of initial conditions with large vertical transport can be 
$O(1)$ in size and these regions can have relatively sharp boundaries. Vertical transport and mixing are important processes in rotating stratified turbulent flow and we expect that QG
$^{+1}$ point vortex point-vortex dynamics will prove useful in understanding AG vertical transport and mixing across many geophysical and astrophysical settings.
It is well known that 2-D and QG point-vortex dynamics have Hamiltonian structures which provide insight into the dynamical behaviour of the vortices. This fact, together with the behaviour of the trajectories seen here, leads us to think it likely that QG
$^{+1}$ point vortex point-vortex dynamics is Hamiltonian. There is a large literature on Hamiltonian fluid dynamics, e.g. Morrison (Reference Morrison1998). Based on the theoretical understanding of Hamiltonian fluid dynamics, we expect that the imposition of balance and the asymptotic expansion in QG
$^{+1}$ will preserve the Hamiltonian structure present in the primitive equations. We also note that the balance formulation of (Holm Reference Holm1996) has an explicit Hamiltonian structure. We leave the development of any potential Hamiltonian structure in QG
$^{+1}$ point vortices for future work.
Acknowledgements
The author would like to acknowledge fascinating and useful discussions with J.-L. Thiffeault, P. J. Morrison and C. R. Doering.
Funding
This work was performed in part at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. This work utilized resources from the University of Colorado Boulder Research Computing Group, which is supported by the National Science Foundation (awards ACI-1532235 and ACI-1532236), the University of Colorado Boulder and Colorado State University.
Declaration of interests
The author reports no conflict of interest.
Data availability statement
Mathematica code to calculate the far-field equations in § 4 is available at https://github.com/JeffreyWeiss/QGplus1PointVortices.
Open access
