# Charge distribution and spin textures in magic-angle twisted bilayer
graphene

A.O. Sboychakov  Institute for Theoretical and Applied Electrodynamics,
Russian Academy of Sciences, Moscow, 125412 Russia  A.V. Rozhkov  Institute
for Theoretical and Applied Electrodynamics, Russian Academy of Sciences,
Moscow, 125412 Russia  A.L. Rakhmanov  Institute for Theoretical and Applied
Electrodynamics, Russian Academy of Sciences, Moscow, 125412 Russia

###### Abstract

We examine the coexisting spin and charge density waves as a possible ground
state of the magic-angle twisted bilayer graphene. When interactions are not
included, the spectrum of the material has 4 (8 if spin is taken into account)
almost flat almost degenerate bands. Interactions break down the degeneracy
forming an order parameter which is usually assumed to be a spin density wave
with a preset spin structure. Here we take into account a possible charge
density wave contribution to the order parameter, that is, inhomogeneous
distribution of the charge density within a twisted graphene supercell. We
also calculate self-consistently the spin structure of the order parameter. We
find that the density wave order is stable in the whole doping range from $-4$
to $+4$ extra electrons per supercell. The spin texture changes from collinear
at zero doping to almost coplanar at finite doping. The density wave order
shows nematic distortion when we dope the system. We demonstrate that the
local spin magnetization is much stronger than the charge density variation,
unless the doping exceeds $3$ extra electrons or holes per supercell.

###### pacs:

73.22.Pr, 73.22.Gk, 73.21.Ac

After the discovery of Mott insulating states NatureMott2018 ;
MottSCNature2019  and superconductivity MottSCNature2019 ; NatureSC2018  in
twisted bilayer graphene (tBLG), this material is intensively studied. In tBLG
one graphene layer is rotated with respect to another one by a twist angle
$\theta$. As a result, the system has a superstructure for certain
commensurate twist angles ourBLGreview2016 . The low-energy electronic
properties substantially depend on $\theta$. At the so-called first magic
angle $\theta_{c}\cong 1^{\circ}$ the energy spectrum is characterized by 4 (8
if spin is included) almost flat bands close to the Fermi level. This flatness
makes the system very susceptible to electron-electron interactions.
Interactions lead to the spontaneous violation of certain symmetries. The
nature of non-superconducting many-body states in tBLG is not yet definitively
known. Different candidates were proposed for the observed many-body
insulating states in the tBLG Philips2018 ; PhysRevB.98.081102 ;
ChiralSDW_SC2018 ; AFMMottSC2019 ; FMorderPRL2019 ; OurtBLGPRB2019 ;
NematicPRB2020 ; PhysRevLett.128.247402 ; PhysRevLett.128.156401 ;
JETPLetters.112.651 . The spin density wave (SDW) is among them
ChiralSDW_SC2018 ; AFMMottSC2019 ; OurtBLGPRB2019 ; NematicPRB2020 .

In our previous papers OurtBLGPRB2019 ; NematicPRB2020  we assumed the planar
SDW order to be the ground state of the system in the doping range $-4<x<4$
electrons per supercell. This choice was not arbitrary. It is known that in
the tBLG at small twist angle the electrons at the Fermi level are located in
the AA regions of the superlattice cell ourBLGreview2016 ; dSPRB . The ground
state of the AA stacked bilayer graphene should be antiferromagnetic (AFM)
AAPRL ; AAPRB2013 . Based on the assumption of the SDW order, our simulations
reproduced qualitatively well the dependence of the conductivity on the doping
observed experimentally. We also showed NematicPRB2020  that doping the system
away from the charge neutrality point leads to the formation of the nematic
state, which is also detected by the STM MottNematicNature2019 ;
KerelskyNematicNature2019 ; STMNature2019 ; WongNematic2020 .

In the present paper we extend the variety of possible order parameters,
considering the non-coplanar SDW order on a backdrop of the inhomogeneous
charge distribution. We show that the spin texture changes from collinear to
almost coplanar with doping. It is demonstrated that the symmetry of the SDW
order is reduced with doping. This leads to the formation of the nematic spin
state. At the same time, the rotational symmetry of the charge density is
almost unaffected by doping. The nematic state can manifests itself in the
non-symmetric spatial distribution of the local density of states, which was
observed experimentally MottNematicNature2019 ; KerelskyNematicNature2019 ;
STMNature2019 ; WongNematic2020 . Finally, we show that local spin density is
much stronger than the charge density variation unless the doping exceeds 3
extra electrons or holes per supercell.

Model Hamiltonian. The periodic superstructure exists in the tBLG when
$\theta$ satisfies the condition

$\displaystyle\cos\theta=(3m_{0}^{2}+3m_{0}r+r^{2}/2)/(3m_{0}^{2}+3m_{0}r+r^{2}),$
(1)

with $m_{0}$ and $r$ being mutually coprime positive integers ourBLGreview2016
. The superlattice cell has a form of a right rhombus. When $r=1$ the
superlattice cell coincides with the moiré cell. In our study we consider only
such superstructures. When the twist angle is small, the supercell can be
considered as consisting of regions with AA, AB, and BA stacking
ourBLGreview2016 ; NanoLettTB .

We use the following model Hamiltonian of tBLG:

$\displaystyle H$ $\displaystyle=$
$\displaystyle\\!\\!\\!\sum_{{\mathbf{nm}ij\atop\alpha\beta\sigma}}\\!t(\mathbf{r}_{\mathbf{n}}^{i\alpha};\mathbf{r}_{\mathbf{m}}^{j\beta})d^{{\dagger}}_{\mathbf{n}i\alpha\sigma}d^{\phantom{{\dagger}}}_{\mathbf{m}j\beta\sigma}+U\\!\sum_{{\mathbf{n}i\alpha\sigma}}\\!n_{\mathbf{n}i\alpha\uparrow}n_{\mathbf{n}i\alpha\downarrow}+$
(2)
$\displaystyle\frac{1}{2}\\!\mathop{{\sum}^{\prime}}_{{\mathbf{nm}ij\atop\alpha\beta\sigma\sigma^{\prime}}}\\!V(\mathbf{r}_{\mathbf{n}}^{i\alpha}-\mathbf{r}_{\mathbf{m}}^{j\beta})n_{\mathbf{n}i\alpha\sigma}n_{\mathbf{m}j\beta\sigma^{\prime}}\,.$

Here $d^{{\dagger}}_{\mathbf{n}i\alpha\sigma}$
($d^{\phantom{{\dagger}}}_{\mathbf{n}i\alpha\sigma}$) are the creation
(annihilation) operators of the electron with spin $\sigma$ ($=\uparrow$,
$\downarrow$) at the unit cell $\mathbf{n}$ in the layer $i$ ($=1,2$) in the
sublattice $\alpha$ ($={\cal A,B}$), while
$n_{\mathbf{n}i\alpha\sigma}=d^{{\dagger}}_{\mathbf{n}i\alpha\sigma}d^{\phantom{{\dagger}}}_{\mathbf{n}i\alpha\sigma}$.
The first term in Eq. (2) is the single-particle tight-binding Hamiltonian
with $t(\mathbf{r}_{\mathbf{n}}^{i\alpha};\mathbf{r}_{\mathbf{m}}^{j\alpha})$
being the amplitude of the electron hopping from site in the position
$\mathbf{r}_{\mathbf{m}}^{j\beta}$ to the site
$\mathbf{r}_{\mathbf{n}}^{i\alpha}$. The second term in Eq. (2) describes the
on-site (Hubbard) interaction of electrons, while the last term corresponds to
the inter-site Coulomb interaction (the prime means that the elements with
$\mathbf{r}_{\mathbf{n}}^{i\alpha}=\mathbf{r}_{\mathbf{m}}^{j\beta}$ should be
excluded).

Now we have to choose a parametrization of the hopping amplitudes. We keep
only nearest-neighbor terms for the intralayer hopping with $t=-2.57$ eV. The
inter-plane hopping amplitudes are parameterized by the following Slater-
Koster formula for $p_{z}$ electrons

$t(\mathbf{r};\mathbf{r}^{\prime})=\frac{\left[(\mathbf{r}-\mathbf{r}^{\prime})\mathbf{e}_{z}\right]^{2}}{|\mathbf{r}-\mathbf{r}^{\prime}|^{2}}V_{\sigma}(\mathbf{r}-\mathbf{r}^{\prime})\,,$
(3)

where $\mathbf{e}_{z}$ is the unit vector transverse to the layers,

$V_{\sigma}(\mathbf{r})=t_{0}e^{-(|\mathbf{r}|-d)/r_{0}}F_{c}(|\mathbf{r}|)\,,\;\;F_{c}(r)=\frac{1}{1+e^{(r-r_{c})/l_{c}}}\,,$
(4)

and the cutoff function $F_{c}(r)$ is introduced to nullify the hopping
amplitudes at distances larger than $r_{c}$. We use $r_{c}=4.92$ Å,
$l_{c}=0.2$ Å. The parameter $t_{0}$ defines the largest interlayer hopping
amplitude. We choose $t_{0}=0.37$ eV (this value was used to describe the AB
bilayer graphene ourBLGreview2016 ). The parameter $r_{0}$ describes how fast
the hopping amplitudes decay inside the region $r<r_{c}$. We choose
$r_{0}=0.34$ Å.

At the first magic angle $\theta=\theta_{c}$ the low-energy band structure
consists of 4 almost flat almost degenerate bands separated by energy gaps
from lower and higher dispersive bands. The bandwidth $W$ of these flat bands
has a minimum at $\theta=\theta_{c}$. For model parameters used here we
numerically determine that $W=1.8$ meV, the gaps between flat and dispersive
bands are $\approx 2.5$ meV, and $\theta_{c}=1.08^{\circ}$ (which is close to
the experimentally observed and corresponds to the superstructure with
$m_{0}=30$ and $r=1$).

Spin and Charge Densities. When the system has flat bands crossing the Fermi
level, the interactions become very important. The interactions break the
symmetry of the single particle Hamiltonian inducing some order parameter. We
start with the SDW-like ordering, with multicomponent order parameter. First,
we introduce on-site order parameters

$\displaystyle\Delta_{\mathbf{n}i\alpha}$ $\displaystyle=$ $\displaystyle
U\langle
d^{{\dagger}}_{\mathbf{n}i\alpha\uparrow}d^{\phantom{{\dagger}}}_{\mathbf{n}i\alpha\downarrow}\rangle\,,$
$\displaystyle\Delta^{z}_{\mathbf{n}i\alpha}$ $\displaystyle=$ $\displaystyle
U\left[\langle
d^{{\dagger}}_{\mathbf{n}i\alpha\uparrow}d^{\phantom{{\dagger}}}_{\mathbf{n}i\alpha\uparrow}\rangle-\langle
d^{{\dagger}}_{\mathbf{n}i\alpha\downarrow}d^{\phantom{{\dagger}}}_{\mathbf{n}i\alpha\downarrow}\rangle\right]\\!/2.$
(5)

The quantity $\Delta_{\mathbf{n}i\alpha}$ is complex, while
$\Delta^{z}_{\mathbf{n}i\alpha}$ is real. They describe on-site magnetization

$\mathbf{S}_{\mathbf{n}i\alpha}=\left[\operatorname{Re}(\Delta_{\mathbf{n}i\alpha}),\,\operatorname{Im}(\Delta_{\mathbf{n}i\alpha}),\,\Delta^{z}_{\mathbf{n}i\alpha}\right]/U.$
(6)

Parameters $\Delta_{\mathbf{n}i\alpha}$ and $\Delta^{z}_{\mathbf{n}i\alpha}$
are controlled by the Hubbard interaction. We take $U=2t$. This value is
somewhat smaller than the critical value for a single-layer graphene
transition into a mean-field AFM state MF_Uc_sorella1992 , $U_{c}=2.23t$.
Thus, our Hubbard interaction is rather strong, but not too strong to open a
gap in single layer graphene.

In a graphene layer, each atom in one sublattice has three nearest neighbors
belonging to another sublattice. For this reason we consider three types of
in-plane nearest-neighbor order parameters of the SDW type,
$A^{(\ell)}_{\mathbf{n}i\sigma}$ and $A^{z(\ell)}_{\mathbf{n}i}$
($\ell=1,\,2,\,3$), corresponding to three different links connecting the
nearest-neighbor sites. These order parameters are defined as follows

$\displaystyle A^{(\ell)}_{\mathbf{n}i\sigma}$ $\displaystyle=$ $\displaystyle
V_{\rm nn}\langle d^{{\dagger}}_{\mathbf{n}+\mathbf{n}_{\ell}i{\cal
A}\sigma}d^{\phantom{{\dagger}}}_{\mathbf{n}i{\cal B}\bar{\sigma}}\rangle\,,$
(7) $\displaystyle A^{z(\ell)}_{\mathbf{n}i}$ $\displaystyle=$
$\displaystyle\frac{V_{\rm nn}}{2}\left(\langle
d^{{\dagger}}_{\mathbf{n}+\mathbf{n}_{\ell}i{\cal
A}\uparrow}d^{\phantom{{\dagger}}}_{\mathbf{n}i{\cal
B}\uparrow}\rangle-\langle d^{{\dagger}}_{\mathbf{n}+\mathbf{n}_{\ell}i{\cal
A}\downarrow}d^{\phantom{{\dagger}}}_{\mathbf{n}i{\cal
B}\downarrow}\rangle\right),$

where $\mathbf{n}_{1}=(0,\,0)$, $\mathbf{n}_{2}=(1,\,0)$,
$\mathbf{n}_{3}=(0,\,1)$, $\bar{\sigma}=-\sigma$, and $V_{\rm
nn}=V(|\bm{\delta}|)$ is the in-plane nearest-neighbor Coulomb repulsion
energy ($\bm{\delta}$ is the vector connecting ${\cal A}$ and ${\cal B}$ sites
in the graphene unit cell). We take $V_{\text{nn}}/U=0.59$, in agreement with
Ref. Wehling . The order parameters (7) define spins on links connecting in-
plane nearest neighbor sites according to

$\displaystyle\\!\\!\mathbf{S}^{(\ell)}_{\mathbf{n}i}=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}\\!\bm{\sigma}_{\sigma\sigma^{\prime}}\langle
d^{{\dagger}}_{\mathbf{n}+\mathbf{n}_{\ell}i{\cal
A}\sigma}d^{\phantom{{\dagger}}}_{\mathbf{n}i{\cal
B}\sigma^{\prime}}\rangle+{\rm c.c.}$
$\displaystyle\\!\\!=\frac{1}{V_{\text{nn}}}\\!\left(\frac{1}{2}\operatorname{Re}[A^{(\ell)}_{\mathbf{n}i\uparrow}+A^{(\ell)}_{\mathbf{n}i\downarrow}],\,\frac{1}{2}\operatorname{Im}[A^{(\ell)}_{\mathbf{n}i\uparrow}-A^{(\ell)}_{\mathbf{n}i\downarrow}],A^{z(\ell)}_{\mathbf{n}i}\right),$

where $\bm{\sigma}$ is a vector of the Pauli matrices. In contrast to Ref.
NematicPRB2020 , the spins $\mathbf{S}_{\mathbf{n}i\alpha}$ and
$\mathbf{S}^{(\ell)}_{\mathbf{n}i}$ are allowed to have the $z$ components. In
Ref. NematicPRB2020  we also considered inter-layer nearest-neighbor SDW order
parameters. The calculations showed, however, that these components are by
order of magnitude smaller than $A^{(\ell)}_{\mathbf{n}i\sigma}$ (which are
smaller than $\Delta_{\mathbf{n}i\alpha}$). In present paper we neglect such
order parameters.

Besides the SDW order parameters, we consider here the charge-density-wave-
like contributions. First, we take into account that the charges are not
uniformly distributed inside the superlattice cell and introduce the quantity

$\Delta^{c}_{\mathbf{n}i\alpha}=\frac{U}{2}\left[\sum_{\sigma}\langle
d^{{\dagger}}_{\mathbf{n}i\alpha\sigma}d^{\phantom{{\dagger}}}_{\mathbf{n}i\alpha\sigma}\rangle-1\right].$
(9)

This value can be considered as the on-site potential due to charge
inhomogeneity. Note that $\Delta^{c}_{\mathbf{n}i\alpha}$ is finite even in
the absence of any symmetry breaking since the sites inside a supercell are
non-identical. Thus, it cannot be considered as an order parameter. In our
simulations $\Delta^{c}_{\mathbf{n}i\alpha}$ is normalized according to

$\frac{2}{U{\cal
N}_{sc}}\sum_{\mathbf{n}i\alpha}\Delta^{c}_{\mathbf{n}i\alpha}=x\,,$ (10)

where ${\cal N}_{sc}$ is the number of supercells in the system, and $x$ is
the doping level, that is, the number of extra electrons or holes per one
supercell. Besides $\Delta^{c}_{\mathbf{n}i\alpha}$, we introduce the in-plane
nearest-neighbor (inter-site) potentials

$A^{c(\ell)}_{\mathbf{n}i}=\frac{V_{\rm nn}}{2}\sum_{\sigma}\langle
d^{{\dagger}}_{\mathbf{n}+\mathbf{n}_{\ell}i{\cal
A}\sigma}d^{\phantom{{\dagger}}}_{\mathbf{n}i{\cal B}\sigma}\rangle\,.$ (11)

For all these charge and spin distributions we assume the same periodicity as
the supercell periodicity.

Figure 1: The spatial distribution of the order parameters ${\cal
D}_{\mathbf{n}1\alpha}$ [(a), (e)] and ${\cal A}^{(\ell)}_{\mathbf{n}1}$ [ (b)
– (d) and (f) – (h)] (for definition, see the text) calculated at $x=0$ [(a) –
(d)] and $x=-2$ [(e) – (h)]. Doping reduces the symmetry of the order
parameters from $C_{6}$ down to $C_{2}$.

Mean-field iteration scheme. To calculate the order parameters and potentials,
we use a mean-field approach. It is based on the decoupling of the quadratic
terms in the Hamiltonian (2):

$\displaystyle\\!\\!\\!\\!\\!\\!n^{\phantom{{\dagger}}}_{\mathbf{n}i\alpha\sigma}n^{\phantom{{\dagger}}}_{\mathbf{m}j\beta\sigma^{\prime}}\\!\rightarrow\\!n^{\phantom{{\dagger}}}_{\mathbf{n}i\alpha\sigma}\langle
n^{\phantom{{\dagger}}}_{\mathbf{m}j\beta\sigma^{\prime}}\rangle+n^{\phantom{{\dagger}}}_{\mathbf{m}j\beta\sigma^{\prime}}\langle
n^{\phantom{{\dagger}}}_{\mathbf{n}i\alpha\sigma}\rangle-$
$\displaystyle\\!\\!\\!\\!\\!\\!\langle
n^{\phantom{{\dagger}}}_{\mathbf{n}i\alpha\sigma}\rangle\langle
n^{\phantom{{\dagger}}}_{\mathbf{m}j\beta\sigma^{\prime}}\rangle-d^{{\dagger}}_{\mathbf{n}i\alpha\sigma}d^{\phantom{{\dagger}}}_{\mathbf{m}j\beta\sigma^{\prime}}\langle
d^{{\dagger}}_{\mathbf{m}j\beta\sigma^{\prime}}d^{\phantom{{\dagger}}}_{\mathbf{n}i\alpha\sigma}\rangle-$
$\displaystyle\\!\\!\\!\\!\\!\\!d^{{\dagger}}_{\mathbf{m}j\beta\sigma^{\prime}}d^{\phantom{{\dagger}}}_{\mathbf{n}i\alpha\sigma}\langle
d^{{\dagger}}_{\mathbf{n}i\alpha\sigma}d^{\phantom{{\dagger}}}_{\mathbf{m}j\beta\sigma^{\prime}}\rangle+$
$\displaystyle\\!\\!\\!\\!\\!\\!\langle
d^{{\dagger}}_{\mathbf{n}i\alpha\sigma}d^{\phantom{{\dagger}}}_{\mathbf{m}j\beta\sigma^{\prime}}\rangle\langle
d^{{\dagger}}_{\mathbf{m}j\beta\sigma^{\prime}}d^{\phantom{{\dagger}}}_{\mathbf{n}i\alpha\sigma}\rangle.$

As a result, we obtain a mean-field Hamiltonian $H^{\text{MF}}$. Numerical
algorithm for finding the order parameters and the potentials is similar to
that used in Ref. NematicPRB2020 . First, we write down the electronic
operators in the momentum representation PankratovPRB2013

$\displaystyle
d^{\phantom{{\dagger}}}_{\mathbf{pG}i\alpha\sigma}=\frac{1}{\sqrt{{\cal{N}}}}\sum_{\mathbf{n}}\exp{[-i(\mathbf{p}+\mathbf{G})\mathbf{r}_{\mathbf{n}}^{i}]}d_{\mathbf{n}i\alpha\sigma}\,,$
(13)

where ${\cal N}$ is the number of the graphene unit cells in one layer of the
sample, $\mathbf{r}_{\mathbf{n}}^{i}$ is the position of the $\mathbf{n}$-th
unit cell of the $i$-th layer, the momentum $\mathbf{p}$ lies in the first
Brillouin zone of the superlattice, and $\mathbf{G}$ are the reciprocal
vectors of the superlattice lying in the first Brillouin zone. The number of
vectors $\mathbf{G}$ is equal to $N_{\rm sc}={\cal N}/{\cal N}_{sc}$ for each
graphene layer.

The mean-field Hamiltonian $\hat{H}_{\mathbf{p}}$ is $N_{R}\times N_{R}$
matrix, where $N_{R}=8N_{\rm sc}$ (factor $8$ is due to the spin, layer, and
sublattice indices). The rank of this matrix is too large to perform numerical
integration over the quasimomentum $\mathbf{p}$ for a realistic time. Indeed,
for the first magic angle $\theta_{c}=1.08^{\circ}$ we have $N_{R}=22328$. For
this reason, we use simplifications. The main contribution to the order
parameters comes from the low-energy states. Consequently, the contributions
from other states can be approximated. In the limit of uncoupled ($t_{0}=0$)
graphene layers and zero order parameters, the matrix $\hat{H}_{\mathbf{p}}$
is block-diagonal with the $2\times 2$ matrices on its diagonal

$-t\left(\begin{array}[]{cc}0&f_{\mathbf{p}+\mathbf{G}}^{i}\\\
f_{\mathbf{p}+\mathbf{G}}^{i*}&0\end{array}\right),$ (14)

where
$f^{i}_{\mathbf{p}}=1+e^{-i\mathbf{pa}^{i}_{1}}+e^{-i\mathbf{pa}^{i}_{2}}$,
$\mathbf{a}^{i}_{1,2}$ are the unit vectors of the $i$-th layer. The
eigenenergies of such a matrix are $\pm t|f_{\mathbf{p}+\mathbf{G}}^{i}|$. The
interlayer hopping amplitudes and order parameters are much smaller than $t$.
As long as we are interested in low-energy features, we can use the truncated
matrix $\hat{H}^{\prime}_{\mathbf{p}}$ excluding the rows and columns in
$\hat{H}_{\mathbf{p}}$ containing elements with
$t|f_{\mathbf{G+p}}^{i}|>E_{R}$, where $E_{R}$ is the cutoff energy. The rank
of the truncated matrix is $N^{\prime}_{R}<N_{R}$. The eigenenergies
$E^{(S)}_{\mathbf{p}}$ of $\hat{H}^{\prime}_{\mathbf{p}}$, which lie close to
$\pm E_{R}$, are calculated with significant errors. To fix this problem, we
take into account only bands with $|E^{(S)}_{\mathbf{0}}|<E_{0}$, where
$E_{0}<E_{R}$. The number of such bands is $N_{0}<N^{\prime}_{R}$. We use
$E_{0}=0.2t$, $E_{R}=0.4t$ ($N_{0}=480$, $N^{\prime}_{R}=720$). Calculations
with smaller and larger $E_{R}$ and $E_{0}$ show that the results are almost
independent of these quantities.

The contribution to the total energy from the discarded states
$E^{(S)}_{\mathbf{p}}<-E_{0}$ must be accounted for separately. Since $E_{0}$
is much larger than order parameters, this can be done perturbatively. The
leading corrections to the total energy are quadratic in SDW order parameters
$\Delta_{\mathbf{n}i\alpha}$, $\Delta^{z}_{\mathbf{n}i\alpha}$,
$A^{(\ell)}_{\mathbf{n}i\sigma}$, and $A^{z(\ell)}_{\mathbf{n}i}$. We assume
that the proportionality coefficients are identical for all order parameters
and are equal to

$-\frac{1}{V_{c}(E_{0})}=-\frac{1}{2}\int_{E_{0}}^{3t}\\!\\!\\!\\!dE\,\frac{\rho_{0}(E)}{E}\,,$
(15)

where $\rho_{0}(E)$ is the single-layer graphene density of states. Such a
correction can be taken into account by the replacement in the mean-field
Hamiltonian

$\displaystyle H^{\text{MF}}$ $\displaystyle\rightarrow$ $\displaystyle
H^{\prime\text{MF}}-\sum_{\mathbf{n}i\alpha}\frac{|\Delta_{\mathbf{n}i\alpha}|^{2}}{V_{c}(E_{0})}-\sum_{\mathbf{n}i\alpha}\frac{|\Delta^{z}_{\mathbf{n}i\alpha}|^{2}}{V_{c}(E_{0})}$
(16)
$\displaystyle-\sum_{\mathbf{n}i{\ell}\sigma}\frac{|A^{(\ell)}_{\mathbf{n}i\sigma}|^{2}}{V_{c}(E_{0})}-\sum_{\mathbf{n}i{z\ell}\sigma}\frac{|A^{z(\ell)}_{\mathbf{n}i}|^{2}}{V_{c}(E_{0})},$

where $H^{\prime\text{MF}}$ is the Hamiltonian in the truncated basis. The
contributions to $\Delta^{c}$ and $A^{c(\ell)}$ from the bands with
$E^{(S)}_{\mathbf{p}}<-E_{0}$ are almost independent of the interactions and
ordering. In the limit of uncoupled layers and zero order parameters, the
contribution to $\Delta^{c}_{\mathbf{n}i\alpha}$ from each excluded band
equals $U/(4N_{sc})$. As a result, we obtain the self-consistency condition
for $\Delta^{c}_{\mathbf{n}i\alpha}$

$\displaystyle\frac{2}{U}\Delta^{c}_{\mathbf{n}i\alpha}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{sc}}\\!\mathop{{\sum}^{\prime}}_{S}\\!\\!\sum_{\mathbf{GG}^{\prime}\sigma}\\!\\!\int\\!\frac{d^{2}\mathbf{p}}{v_{RBZ}}\Phi^{(S)*}_{\mathbf{pG}i\alpha\sigma}\Phi^{(S)}_{\mathbf{pG}^{\prime}i\alpha\sigma}\times$
(17) $\displaystyle
e^{-i(\mathbf{G}-\mathbf{G}^{\prime})\mathbf{r}_{\mathbf{n}}^{i}}\Theta\left(\mu-E^{(S)}_{\mathbf{p}}\right)-N_{0}/8N_{sc},$

where $\Phi^{(S)}_{\mathbf{pG}i\alpha\sigma}$ are eigenfunctions of
$\hat{H}^{\prime}_{\mathbf{p}}$, $\Theta(x)$ is the step function, $v_{RBZ}$
is the area of the supercell’s Brillouin zone, and $\mu$ is the chemical
potential. The summation in Eq. (17) is performed over $N_{0}$ bands and
truncated basis, and the equation for chemical potential is

$\displaystyle\mathop{{\sum}^{\prime}}_{S}\\!\\!\int\\!\frac{d^{2}\mathbf{p}}{v_{RBZ}}\Theta\left(\mu-E^{(S)}_{\mathbf{p}}\right)-\frac{N_{0}}{2}=x\,.$
(18)

Figure 2: The on-site spins $\mathbf{S}_{\mathbf{n}1\alpha}$ in sublattice
${\cal A}$ (a) and ${\cal B}$ (b) at $x=-2$. Spins lie in the $xz$ plane. For
the purpose of the visual representation, the rotation of the spins from $xz$
plane to $xy$ plane is performed. Only AA region of the supercell is shown.

Figure 3:  The textures of spins on the link,
$\mathbf{S}^{(\ell)}_{\mathbf{n}1}$, calculated at $x=-2$. Panels (a) – (c)
show the $xy$ projections of the spins for different values of link parameter
$\ell=1,2,3$, while panels (d) – (f) depicts the $xz$ projections (after
rotation to the $xy$ plane). Only AA region of the supercell is shown.

The contribution to $A^{c(\ell)}_{\mathbf{n}i}$ from excluded bands is non-
zero (moreover, one can check that it is non-zero even for single-layer
graphene). In the limit of uncoupled layers this contribution is independent
of $\mathbf{n}$ and $i$. Thus, one can write
$A^{c(\ell)}_{\mathbf{n}i}=\bar{A}+\delta\\!A^{c(\ell)}_{\mathbf{n}i}$, where
$\bar{A}$ can be estimated as

$\bar{A}=\frac{V_{\text{nn}}}{2}\int\frac{d^{2}\mathbf{k}}{v_{BZ}}e^{-i\varphi_{\mathbf{k}}^{i}}\,\Theta\left(t|f_{\mathbf{k}}^{i}|-E_{0}\right),$
(19)

where $\varphi_{\mathbf{k}}^{i}=\arg(f_{\mathbf{k}}^{i})$, and $\bar{A}$ is
independent of $i$. For parameters chosen, we have
$\bar{A}/V_{\text{nn}}=0.26$. The parameter $\bar{A}$ renormalizes the in-
plane nearest-neighbor hopping amplitude $t$. We assume that this
renormalization is already absorbed into the value $t=-2.57$ eV, and thus can
be ignored. As for $\delta\\!A^{c(\ell)}_{\mathbf{n}i}$, we add quadratic term
to the effective Hamiltonian, similar to that in Eq. (16).

To calculate the order parameters and the potentials we minimize the total
energy ${\cal E}$ iteratively using the method of successive approximations.
In each iteration, we calculate $E^{(S)}_{\mathbf{p}}$ and
$\Phi^{(S)}_{\mathbf{pG}i\alpha\sigma}$ of the matrix
$\hat{H}^{\prime}_{\mathbf{p}}$, and the gradients $\partial\cal
E/\partial\lambda=\langle\partial H^{\text{MF}}/\partial\lambda\rangle$, where
$\lambda=\Delta_{\mathbf{n}i\alpha}$, $\Delta^{z}_{\mathbf{n}i\alpha}$,
$A^{(\ell)}_{\mathbf{n}i\sigma}$, $A^{z(\ell)}_{\mathbf{n}i}$, or
$\delta\\!A^{c(\ell)}_{\mathbf{n}i}$. These gradients are used to calculate
new values of the aforementioned quantities. New value of
$\Delta^{c}_{\mathbf{n}i\alpha}$ is found using Eq. (17) with $\mu$ from Eq.
(18).

Results. Our numerical calculations reveal that the SDW order exists in the
system at any doping level. However, the spatial distributions of the SDW
magnetization are qualitatively different for different $x$.

Let us start from $x=0$. Figures 1(a) – (d) depicts the spatial distributions
of ${\cal D}_{\mathbf{n}i\alpha}=U|\mathbf{S}_{\mathbf{n}i\alpha}|$ and ${\cal
A}^{(\ell)}_{\mathbf{n}i}=V_{\text{nn}}|\mathbf{S}^{(\ell)}_{\mathbf{n}i}|$
inside the superlattice cell for $i=1$. These quantities describe the absolute
values of the on-site spins and spins on the links, correspondingly. We see
that the order parameters are non-zero only in the AA region. At zero doping,
the order parameters spatial distributions have rather symmetric form. The
area of non-zero ${\cal D}_{\mathbf{n}1\alpha}$ has a shape of a ‘rounded’
hexagon. This hexagon is invariant under rotation on $60^{\circ}$ around the
center of the AA region $\mathbf{R}_{0}=(\mathbf{R}_{1}+\mathbf{R}_{2})/3$,
with $\mathbf{R}_{1,2}$ being the elementary vectors of the superlattice.
Figures 1(b) – (d) demonstrate that spin magnetization on the links remain
finite within the areas shaped like a dumbbell. The rotation on $180^{\circ}$
and also on $120^{\circ}$ around $\mathbf{R}_{0}$ preserves ${\cal
A}^{(\ell)}_{\mathbf{n}1}$.

When $x=0$ all spins, $\mathbf{S}_{\mathbf{n}i\alpha}$ and
$\mathbf{S}^{(\ell)}_{\mathbf{n}i}$, are collinear. For layer $1$, all on-site
spins $\mathbf{S}_{\mathbf{n}1{\cal A}}$ of the sublattice ${\cal A}$ point
along $x$ axis, while in the sublattice ${\cal B}$ they have the opposite
direction. The same is true for layer $2$. Therefore, the on-site spin texture
has the collinear AFM arrangement. The spins on the links
$\mathbf{S}^{(\ell)}_{\mathbf{n}i}$ also form a kind of AFM structure: in one
part of the dumbbell $\mathbf{S}^{(\ell)}_{\mathbf{n}i}$ directed along $x$
axis, while in another part of the dumbbell they have the opposite direction.

Figure 4:  The spatial distribution of the on-site potential
$\Delta^{c}_{\mathbf{n}1{\cal A}}$ [(a), (e)], $\Delta^{c}_{\mathbf{n}1{\cal
B}}$ [(b), (f)], $\Delta^{c}_{\mathbf{n}2{\cal A}}$ [(c), (g)], and
$\Delta^{c}_{\mathbf{n}2{\cal B}}$ [(d), (h)], calculated at $x=0$ [(a) – (d)]
and at $x=-2$ [(e) – (h)].

The order parameters distributions change with doping. In Fig. 1(e) – (h) we
plot the spatial distribution of ${\cal D}_{\mathbf{n}1\alpha}$ and ${\cal
A}^{(\ell)}_{\mathbf{n}1}$ at hole doping $x=-2$. We see that now no hexagonal
symmetry exists. The profile of ${\cal D}_{\mathbf{n}1\alpha}$ is stretched
unidirectionally, and $120^{\circ}$ rotation is no longer a symmetry of ${\cal
A}^{(\ell)}_{\mathbf{n}1}$. However, all these order parameters are symmetric
under rotation on $180^{\circ}$. In doped systems the spin textures are no
longer collinear. The on-site spins turn out to be coplanar. We believe that
this is not an artifact since our numerical procedure accounts for non-
coplanar textures. In our simulations the on-site spins
$\mathbf{S}_{\mathbf{n}i\alpha}$ lie in the $xz$ plane (see Fig. 2), they form
a helical AFM structure. The spins on links,
$\mathbf{S}^{(\ell)}_{\mathbf{n}i}$, form a non-coplanar structures. However,
almost all of these spins lie in the $xz$ plane, and only small fraction of
them have $y$ components. Figures 3(a) – (c) show the directions of spins
$\mathbf{S}^{(\ell)}_{\mathbf{n}i}$ projected to the $xy$ plane, while Figs.
3(d) – (f) show the directions of spins projected to the $xz$ plane (with
subsequent rotation to the $xy$ plane). We see that almost all spins lie in
the $xz$ plane. The spins which violate the coplanarity lie along three lines
passing through the center of the AA region. These lines are visible in Fig.
3(a) – (c).

Figure 5:  The spatial distribution of the inter-site potentials $\delta
A^{c(\ell)}_{\mathbf{n}1}$ for three different $\ell$ , calculated at x=0 [(a)
– (c)] and at $x=-2$ [(d) – (f)].

In Fig. 4 we plot the spatial distributions of the on-site potential
$\Delta^{c}_{\mathbf{n}i\alpha}$ calculated at zero doping [panels (a) – (d)]
and at $x=-2$ [panels (e) – (h)]. The distributions are shown for each layer
and sublattice separately. At $x=0$ there is an excess of electrons in the AA
region and lack of electrons in the AB (corners of the supercell) and BA (area
centered at $\mathbf{R}_{BA}=2(\mathbf{R}_{1}+\mathbf{R}_{2})/3$) regions of
the superlattice cell. At $x=-2$, the area, where excess of electrons is
observed, expands; extra electrons also appear in the AB (for
$\Delta^{c}_{\mathbf{n}1{\cal B}}$ and $\Delta^{c}_{\mathbf{n}2{\cal A}}$) and
BA (for $\Delta^{c}_{\mathbf{n}1{\cal A}}$ and $\Delta^{c}_{\mathbf{n}2{\cal
B}}$) regions. At zero doping the charge distributions remain invariant under
rotation on $120^{\circ}$ around the point $\mathbf{R}_{0}$. The potentials
$\Delta^{c}_{\mathbf{n}i{\cal A}}$ and $\Delta^{c}_{\mathbf{n}i{\cal B}}$ are
transformed into each other under the rotation on $60^{\circ}$ around the
point $\mathbf{R}_{0}$. The same is true approximately at $x=-2$. Thus, the
nematicity does not manifest itself through the charge distribution.

Figure 5 shows the spatial distributions of the inter-site potentials $\delta
A^{c(\ell)}_{\mathbf{n}1}$, calculated at $x=0$ and $x=-2$. The profiles of
$\delta A^{c(\ell)}_{\mathbf{n}1}$ are approximately the same for $x=0$ and
$x=-2$, and only absolute values of $\delta A^{c(\ell)}_{\mathbf{n}1}$ change.
These spatial distributions are approximately symmetric under rotation on
$180^{\circ}$ around $\mathbf{R}_{0}$. The inter-site potentials are
transformed into each other under the rotation on $120^{\circ}$ around the
point $\mathbf{R}_{0}$.

Thus, the doping reduces the symmetry of the SDW order from $C_{6}$, which is
the symmetry of the lattice, to $C_{2}$ indicating the appearance of the
nematic state NematicPRB2020 . At the same time the symmetry of the charge-
related quantities $\Delta^{c}$ and $A^{c(\ell)}$ is virtually unaffected by
doping. These parameters are integrated quantities that are influenced by
contributions of high-energy states, the latter being insensitive to the low-
energy symmetry breaking. At the same time, the electron nematic state can
manifest itself as a symmetry reduction of the local density of states, which
can be detected in an STM experiment. This is confirmed by our simulations.
Figure 6 shows the spatial distribution of the local density of states
calculated close to half filling (at $x=-2$ the density of states is almost
zero). The spatial profile is stretched indicating the appearance of the
nematic state. Such a feature was observed in the experiments
MottNematicNature2019 ; KerelskyNematicNature2019 ; STMNature2019 ;
WongNematic2020 .

Finally, we compare the values of the spin and charge density variations. At
zero doping, the on-site spin order parameter ${\cal D}_{\mathbf{n}i\alpha}$
reaches the maximum about $34$ meV. At the same time, the on-site charge
density variation $\Delta^{c}_{\mathbf{n}i\alpha}$ is about $6$ times smaller.
The SDW order parameters decrease down to zero, when the doping changes from
$x=0$ to $x=\pm 4$. At the same time, the values of $\Delta^{c}$ do not change
much with $x$. However, even at $x=\pm 2$ the maximum value of ${\cal
D}_{\mathbf{n}i\alpha}$ is about $3$ times larger than that of
$|\Delta^{c}_{\mathbf{n}i\alpha}|$, and only at $|x|\approx 3$, the on-site
SDW order parameter becomes comparable to $|\Delta^{c}_{\mathbf{n}i\alpha}|$.

Figure 6: The local density of states calculated close to half filling
($x=-1.75$). The spatial profile is stretched indicating the appearance of the
electron nematic state.

Conclusions. We study numerically the coexisting spin and charge density waves
in the magic-angle tBLG in the doping range $x=\pm 4$ extra electrons per
supercell using the mean-field approach. The single-electron spectrum of the
material has 8 almost flat almost degenerate bands. The electron-electron
interaction breaks down the symmetries of the single-particle state forming a
set of order parameters. We calculate self-consistently the charge
distribution in the supercell and spin structure of the SDW order parameters.
We found that the SDW order is stable in the whole doping range. The spin
texture of the SDW order parameters depends crucially on $x$ changing from
collinear at $x=0$ to almost coplanar at finite doping. The ground state of
the doped system has the nematic symmetry.

We are grateful to the Joint Supercomputer Center of the Russian Academy of
Sciences (JSCC RAS) for the computational resources provided. The data
analysis and analytical calculations were funded by the Russian Science
Foundation (project No. 22-22-00464 https://rscf.ru/en/project/22-22-00464/).

## References

  * (1) Y. Cao, V. Fatemi, A. Demir et al., Nature 556, 80 (2018).
  * (2) X. Lu, P. Stepanov, W. Yang et al., Nature 574, 653 (2019).
  * (3) Y. Cao, V. Fatemi, S. Fang et al., Nature 556, 43 (2018).
  * (4) A. Rozhkov, A. Sboychakov, A. Rakhmanov, and F. Nori, Phys. Rep. 648, 1 (2016).
  * (5) B. Padhi, C. Setty, and P. W. Phillips, Nano letters 18, 6175 (2018).
  * (6) M. Ochi, M. Koshino, and K. Kuroki, Phys. Rev. B 98, 081102 (2018).
  * (7) C.-C. Liu, L.-D. Zhang, W.-Q. Chen, and F. Yang, Phys. Rev. Lett. 121,217001 (2018).
  * (8) T. Huang, L. Zhang, and T. Ma, Science Bulletin 64, 310 (2019).
  * (9) K. Seo, V. N. Kotov, and B. Uchoa, Phys. Rev. Lett. 122, 246402 (2019).
  * (10) A. O. Sboychakov, A. V. Rozhkov, A. L. Rakhmanov, and F. Nori, Phys. Rev. B 100, 045111 (2019).
  * (11) A. O. Sboychakov, A. V. Rozhkov, A. L. Rakhmanov, and F. Nori, Phys. Rev. B 102, 155142 (2020).
  * (12) S. Zhang, X. Lu, and J. Liu, Phys. Rev. Lett. 128, 247402 (2022).
  * (13) G. Wagner, Y. H. Kwan, N. Bultinck et al., Phys. Rev. Lett. 128, 156401 (2022).
  * (14) A. O. Sboychakov, A. V. Rozhkov, K. I. Kugel, and A. L. Rakhmanov, JETP Letters 112, 651 (2020).
  * (15) J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. B 86, 155449 (2012).
  * (16) A. L. Rakhmanov, A. V. Rozhkov, A. O. Sboychakov, and F. Nori, Phys. Rev. Lett. 109, 206801 (2012).
  * (17) A. O. Sboychakov, A. V. Rozhkov, A. L. Rakhmanov, and F. Nori, Phys. Rev. B 88, 045409 (2013).
  * (18) Y. Choi, J. Kemmer, Y. Peng et al., Nature Physics 15, 1174 (2019).
  * (19) A. Kerelsky, L. J. McGilly, D. M. Kennes et al., Nature 572, 95 (2019).
  * (20) Y. Xie, B. Lian, B. Jäck, et al., Nature 572, 101 (2019).
  * (21) D. Wong, K. P.Nuckolls, M. Ohet et al., Nature 582, 198 (2020).
  * (22) G. Trambly de Laissardière, D. Mayou, and L. Magaud, Nano Letters 10, 804 (2010).
  * (23) S. Sorella and E. Tosatti, Europhysics Letters 19, 699 (1992).
  * (24) T. O. Wehling, E. Şaşıoğglu, C. Friedrich et al., Phys. Rev. Lett. 106, 236805 (2011).
  * (25) S. Shallcross, S. Sharma, and O. Pankratov, Phys. Rev. B 87, 245403 (2013).


# Point vortex dynamics in three-dimensional ageostrophic balanced flows

Jeffrey B. Weiss1<EMAIL_ADDRESS>1Department of Atmospheric and
Oceanic Sciences, University of Colorado, Boulder, CO 80302 USA

###### Abstract

Geophysical turbulent flows, characterized by rapid rotation, quantified by
small Rossby number, and stable stratification, often self-organize into a
collection of coherent vortices, referred to as a vortex gas. The lowest order
asymptotic expansion in Rossby number is quasigeostrophy which has purely
horizontal velocities and cyclone-anticylone antisymmetry. Ageostrophic
effects are important components of many geophysical flows and, as such, these
phenomena are not well-modeled by quasigeostrophy. The next order correction
in Rossby number, which includes ageostrohpic affects, is so-called balanced
dynamics. Balanced dynamics includes ageostrophic vertical velocity and breaks
the geostrophic cyclone-anticylone antisymmetry. Point vortex solutions are
well known in 2d and quasigeostrophic dynamics and are useful for studying the
vortex gas regime of geophysical turbulence. Here we find point vortex
solutions in fully three dimensional continuously stratified QG+1 dynamics, a
particular formulation of balanced dynamics. Simulations of QG+1 point
vortices show several interesting features not captured by quasigeostrophic
point vortices including significant vertical transport on long timescales.
The ageostrophic component of QG+1 point vortex dynamics renders them useful
in modeling flows where quasigeostrophy filters out important physical
processes.

##  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, 1984; McWilliams
et al., 1994; McWilliams & Weiss, 1994; Dritschel et al., 1999; Carton, 2001;
Viúdez, 2016). A population of physically well-separated coherent vortices, a
“vortex gas”, can be approximated using a hierarchy of idealized models (Kida,
1981; Melander et al., 1986; Meacham, 1992; Meacham et al., 1998; Miyazaki et
al., 2001; Dritschel et al., 2004). The simplest member of the hierarchy is
the point vortex model (Helmholtz, 1867; Kirchhoff, 1876; Aref, 2007) in which
the vorticity is considered to be concentrated at a collection of isolated
points in space. Point vortex models are useful idealizatons for geophysical
and astrophysical turbulence and are used to model atmospheres and oceans on
the Earth, planets and exoplanets, and astrophysical disks (Aref & Stremler,
2001; Carton, 2001; Lucarini et al., 2014; Hirt et al., 2018; G. Abrahamyan,
2020).

The rotation in planetary fluids is characterized by a non-dimensional Rossby
number, $\mbox{{Ro}}=U/fL$, the ratio of the planetary rotation timescale
$1/f$ to the advective timescale $U/L$, where $U$ is the advective velocity
scale, $L$ is a characteristic length scale, $f=2\Omega\sin(\theta_{lat})$ is
the Coriolis parameter, $\Omega$ is the planetary rotation rate, and
$\theta_{lat}$ is the latitude. A flow whose advective timescale is much
longer than its rotation timescale has $\mbox{{Ro}}\ll 1$. Large scale
midlatitude flows in the Earth’s atmosphere typically have $\mbox{{Ro}}\sim
0.1$, flows in the Earth’s oceans typically have $\mbox{{Ro}}\sim 0.01$, while
flows on planets and exoplanets often have similarly small Rossby numbers
(Showman et al., 2010). 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-know 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(\mbox{{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 (1980);
Allen (1993); Holm (1996); Muraki et al. (1999), and references therein).
These models have $O(\mbox{{Ro}}^{2})$ differences and thus are, formally,
equivalently correct through $O(\mbox{{Ro}})$. Different approaches to balance
models choose different quantities as a distinguished variable. For example,
(Allen, 1993) chooses pressure. Here we investigate the balance model of
Muraki et al. (1999) (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 nonzero $O(\mbox{{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 (2015); McKiver & Dritschel (2016); Mahdinia et
al. (2017); Reinaud & Dritschel (2018); McKiver (2020); Sokolovskiy et al.
(2020), and references therein).

Point vortex models have a long history in 2d fluid dynamics (see Aref (2007)
for a review). They are also well developed in layered and continuously
stratified QG dynamics (Gryanik, 1983, 1991; Gryanik et al., 2000, 2006;
Reznik & Kizner, 2007). Due to the constraints of QG dynamics, QG point vortex
models are limited to modeling 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 2D and QG point vortices, but go
beyond them by inheriting from balance models $O(\mbox{{Ro}})$ vertical
velocities and broken cyclone-anticyclone antisymmetry. In section 2, we
review the aspects of QG+1 theory relevant for the development of QG+1 point
vortex dynamics. In section 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, section 4.2 shows that the self-advection is zero
in the point vortex limit, and section 4.3 summarizes the QG+1 point vortex
equations of motion. Section 5 describes some of the properties of QG+1 point
vortex dynamics obtained by analyzing 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, three-dimensional
hydrostatic fluid on an $f$-plane (i.e. constant Coriolis parameter $f$)
described by the primitive equations. We restrict ourselves to flows with a
constant Brunt-Väisälä frequency $N$. The dynamics is nondimensionalized as in
MSR1999. The primitive equations are asyptotically expanded in terms of the
small Rossby number Ro. The leading order QG fields will be denoted by the
subscript “0”.111We note that, technically, QG appears at $O(\mbox{{Ro}})$ in
the full asymptotic expansion. The $O(1)$ term is uninteresting and is
typically removed, following which the QG quantities are scaled by Ro and are
rendered $O(\mbox{{Ro}}^{0})$. At next order, ageostrophic (AG) affects enter.
AG fields will be denoted by the subscript “1”.

The QG+1 model is a specific instance of an iterated balance model (Allen,
1993). 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 (1993) chose pressure as the
distinguished variable. QG+1 chooses potential vorticity (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 $\psi_{0}$ is obtained by
solving a Poisson equation with $q_{0}$ as the source,
$\nabla^{2}\psi_{0}=q_{0}$, where
$\nabla=(\partial_{x},\partial_{y},\partial_{z})$ is the three-dimensional
gradient operator. The QG potential vorticity is thus the sum of the relative
vorticity $\left(\partial_{xx}+\partial_{yy}\right)\psi_{0}$ and stretching
vorticity $\partial_{zz}\psi_{0}$. Conservation of QG potential vorticity
under QG flow is given by

$\frac{D_{0}q_{0}}{Dt}=\partial_{t}q_{0}+\left(\boldsymbol{u}_{0}\cdot\nabla\right)q_{0}=0,$
(1)

where $D_{0}/Dt$ indicates the material derivative of fluid parcels advected
with the QG velocity.

MSR1999 develop the QG+1 approximation by including next order terms in Ro,
making use of the structure of the primitive and QG equations. Two additional
scalar potentials are introduced, $F$ and $G$, from which, with $\psi$, a
vector potential $\boldsymbol{A}$ is constructed which determines the
velocity,

$\boldsymbol{u}=\nabla\times\boldsymbol{A},\qquad\boldsymbol{A}=\begin{pmatrix}-G\\\
\hphantom{-}F\\\ -\psi\\\ \end{pmatrix},$ (2)

where the signs in $F$ and $G$ are chosen to match those in MSR1999. The
velocity and vector potential are expanded asymptotically in Rossby number

$\displaystyle\boldsymbol{u}$ $\displaystyle=\boldsymbol{u}_{0}+\mbox{{Ro}}\
\boldsymbol{u}_{1}+O(\mbox{{Ro}}^{2})$ (3) $\displaystyle\boldsymbol{A}$
$\displaystyle=\boldsymbol{A}_{0}+\mbox{{Ro}}\
\boldsymbol{A}_{1}+O(\mbox{{Ro}}^{2})$ (4)

where the vector potentials are

$\boldsymbol{A}_{0}=\begin{pmatrix}\hphantom{-}0\\\ \hphantom{-}0\\\
-\psi_{0}\\\ \end{pmatrix},\quad\boldsymbol{A}_{1}=\begin{pmatrix}-G_{1}\\\
\hphantom{-}F_{1}\\\ -\psi_{1}\\\ \end{pmatrix},$ (5)

i.e., $F_{0}=G_{0}=0$. The velocities and vector potentials are related by
$\boldsymbol{u}_{n}=\nabla\times\boldsymbol{A}_{n}$ at each order $n$. The
curl relation between the velocity and the vector potential guarantees that
the fluid is incompressible. The potential temperature is given by the
negative divergence of the vector potential,
$\theta=-\nabla\cdot\boldsymbol{A}=\theta_{0}+\mbox{{Ro}}\ \theta_{1}$.

The dynamics is governed by potential vorticity conservation, where the full
potential vorticity $q$ is advected by the approximate velocity

$\frac{Dq}{Dt}=\frac{\partial q}{\partial
t}+\left((\boldsymbol{u}_{0}+\mbox{{Ro}}\
\boldsymbol{u}_{1}+O(\mbox{{Ro}}^{2}))\cdot\nabla\right)q=0$ (6)

The potential vorticity can be partitioned into QG and AG contributions,
$q=q_{0}+\mbox{{Ro}}\ q_{1}+O(\mbox{{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}$.

The QG and AG vector potentials are solutions to Poisson equations,
$\nabla^{2}\boldsymbol{A}_{n}=\boldsymbol{S}_{n}$. The QG source is
$\boldsymbol{S}_{0}=(0,0,-q_{0})$. MSR1999 showed that the QG+1 source is an
operator acting on the QG streamfunction $\psi_{0}$:

$\displaystyle\boldsymbol{S}_{1}$
$\displaystyle=\boldsymbol{S}_{1}(\psi_{0},\psi_{0})=\begin{pmatrix}-2J\left(\partial_{z}\psi_{0},\partial_{y}\psi_{0}\right),\\\
\hphantom{-}2J\left(\partial_{z}\psi_{0},\partial_{x}\psi_{0}\right),\\\
-q_{1}+\left(\nabla^{2}\psi_{0}\right)\partial_{zz}\psi_{0}-\left|\nabla\left(\partial_{z}\psi_{0}\right)\right|^{2},\end{pmatrix}.$
(7)

where $J$ is the usual horizontal Jacobian operator,
$J(f,g)=(\partial_{x}f)(\partial_{y}g)-(\partial_{y}f)(\partial_{x}g)$.
$\boldsymbol{S}_{1}$ thus takes the form of a bilinear non-symmetric vector
operator,

$\displaystyle\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),$
(8)
$\displaystyle\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}),$
(9) $\displaystyle\boldsymbol{S}_{1}(f,g)\neq\boldsymbol{S}_{1}(g,f),$ (10)

where $f$’s and $g$’s are scalar functions of $\boldsymbol{x}$ and $c$’s are
constants. Note that since the AG vector depends quadratically on the QG
streamfunction, it manifestly breaks the cyclone-anticylone (anti)symmetry of
QG:

$\displaystyle\boldsymbol{A}_{0}\xrightarrow[q_{0}\to-
q_{0}]{}-\boldsymbol{A}_{0},$ (11)
$\displaystyle\boldsymbol{A}_{1}\xrightarrow[q_{0}\to-
q_{0}]{}+\boldsymbol{A}_{1},$ (12)

and so $\boldsymbol{A}=\boldsymbol{A}_{0}+\mbox{{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 three-dimensional infinite
domain and require that the velocity and potential temperature go to zero at
infinity.

The QG+1 equations completed by solvability constraints, which, in the
infinite domain consider here, are

$\displaystyle\iiint q_{0}\ d^{3}\boldsymbol{x}$
$\displaystyle=\lim_{z\to\infty}\iint\left.\theta_{0}\right|^{+z}_{-z}dxdy,$
(13) $\displaystyle\iiint q_{1}\ d^{3}\boldsymbol{x}$
$\displaystyle=\lim_{z\to\infty}\iint\left[\left.\theta_{1}+\partial_{z}\psi_{0}\left(\partial_{xx}\psi_{0}+\partial_{yy}\psi_{0}\right)\right]\right|^{+z}_{-z}dxdy,$
(14)

where equation 13 guarantees the invertibility of the QG potential vorticity-
streamfunction relation, and Eq 14 guarantees the invertibility of the
$z$-component of equation 7.

##  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 $\Gamma_{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 shapes,
and their circulations are assumed to be constant. The PV field is then

$q(\boldsymbol{x})=\sum_{i=1}^{N}\Gamma_{i}\hat{q}(\boldsymbol{x}-\boldsymbol{x}_{i}(t)),$
(15)

where $\hat{q}$ is normalized to have unit circulation,
$\int\hat{q}\,d^{3}\boldsymbol{x}=1$. Potential vorticity conservation,
equation (6), then requires the vortex locations move with the local velocity,

$\frac{d\boldsymbol{x}_{i}(t)}{dt}=\boldsymbol{u}(\boldsymbol{x}_{i})=\boldsymbol{u}_{0}(\boldsymbol{x}_{i})+\mbox{{Ro}}\
\boldsymbol{u}_{1}(\boldsymbol{x}_{i})+O(\mbox{{Ro}}^{2}),$ (16)

where we have used

$\frac{\partial\hat{q}(\boldsymbol{x}-\boldsymbol{x}_{i}(t))}{\partial
t}=-\left.\nabla\hat{q}(\boldsymbol{x})\right|_{\boldsymbol{x}=\boldsymbol{x}-\boldsymbol{x}_{i}}\cdot\frac{d\boldsymbol{x}_{i}(t)}{dt}$
(17)

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

$\displaystyle F_{0}(\boldsymbol{x},t)$
$\displaystyle=\sum_{i=1}^{N}\Gamma_{i}\hat{F}_{0}(\boldsymbol{x}-\boldsymbol{x}_{i}),$
(18) $\displaystyle F_{1}(\boldsymbol{x},t)$
$\displaystyle=\sum_{i,j=1}^{N}\Gamma_{i}\Gamma_{j}\hat{F}_{1}(\boldsymbol{x}-\boldsymbol{x}_{i},\boldsymbol{x}-\boldsymbol{x}_{j}).$
(19)

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
2d and QG point vortices. In 2d 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 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}$

$\displaystyle\hat{F}_{1}^{s}(\boldsymbol{x}-\boldsymbol{x}_{i})$
$\displaystyle=\hat{F}_{1}(\boldsymbol{x}-\boldsymbol{x}_{i},\boldsymbol{x}-\boldsymbol{x}_{i}),$
(20)
$\displaystyle\hat{F}_{1}^{p}(\boldsymbol{x}-\boldsymbol{x}_{i},\boldsymbol{x}-\boldsymbol{x}_{j})$
$\displaystyle=\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})$
(21)

resulting in

$F_{1}(\boldsymbol{x},t)=\sum_{i=1}^{N}\Gamma_{i}^{2}\hat{F}_{1}^{s}(\boldsymbol{x}-\boldsymbol{x}_{i})+\sum_{\begin{subarray}{c}i=1\\\
j=i+1\end{subarray}}^{N}\Gamma_{i}\Gamma_{j}\hat{F}_{1}^{p}(\boldsymbol{x}-\boldsymbol{x}_{i},\boldsymbol{x}-\boldsymbol{x}_{j})$
(22)

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:

$\displaystyle\boldsymbol{\hat{S}}_{0}(\boldsymbol{x})=\begin{pmatrix}0\\\
0\\\
-\nabla^{2}\hat{\psi}_{0}(\boldsymbol{x})\end{pmatrix}=\begin{pmatrix}0\\\
0\\\ -\hat{q}(\boldsymbol{x})\end{pmatrix}$ (23)
$\displaystyle\boldsymbol{\hat{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}))$
(24)
$\displaystyle\boldsymbol{\hat{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$ (25)
$\displaystyle\nabla^{2}\boldsymbol{\hat{A}}_{n}^{\\{s,p\\}}=\hat{\boldsymbol{S}}_{n}^{\\{s,p\\}},$
(26)
$\displaystyle\hat{\boldsymbol{u}}_{n}^{\\{s,p\\}}=\nabla\times\boldsymbol{\hat{A}}_{n}^{\\{s,p\\}},$
(27)

where the operator $\boldsymbol{S}_{1}$ is given by equation (7).

Vortex-gas advection, equation (16), becomes

$\frac{d\boldsymbol{x}_{i}(t)}{dt}=\sum_{j=1}^{N}\Gamma_{j}\hat{\boldsymbol{u}}_{0}(\boldsymbol{x}_{i}-\boldsymbol{x}_{j})\\\
+\mbox{{Ro}}\left(\sum_{j=1}^{N}\Gamma_{j}^{2}\hat{\boldsymbol{u}}_{1}^{s}(\boldsymbol{x}_{i}-\boldsymbol{x}_{j})+\sum_{\begin{subarray}{c}j=1\\\
k=j+1\end{subarray}}^{N}\Gamma_{j}\Gamma_{k}\hat{\boldsymbol{u}}_{1}^{p}(\boldsymbol{x}_{i}-\boldsymbol{x}_{j},\boldsymbol{x}_{i}-\boldsymbol{x}_{k})\right).$
(28)

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 2d 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 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
potential vorticity. We will show below that this choice is consistent with
the solvability constraint equation (14). 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\pi|\boldsymbol{x}|$, and the
QG vector potential shape is $\boldsymbol{\hat{A}}_{0}=(0,0,-\hat{\psi}_{0})$.
The single-vortex AG Poisson source $\boldsymbol{\hat{S}}_{1}^{s}$ from a
vortex at the origin can be directly calculated, leading to

$\boldsymbol{\hat{S}}_{1}^{s}=\frac{1}{8\pi^{2}|\boldsymbol{x}|^{8}}\begin{pmatrix}3xz\\\
3yz\\\ -\left(x^{2}+y^{2}+4z^{2}\right)/2\\\ \end{pmatrix}$ (29)

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 $\boldsymbol{\hat{S}}_{1}^{s}$ is proportional to
$|\boldsymbol{x}|^{8}$, we expect the denominator of
$\boldsymbol{\hat{A}}_{1}^{s}$ to be proportional to $|\boldsymbol{x}|^{6}$.
Dimensionally, $\boldsymbol{\hat{S}}_{1}^{s}\sim 1/L^{6}$, so
$\boldsymbol{\hat{A}}_{1}^{s}\sim 1/L^{4}$. This means the numerator of
$\boldsymbol{\hat{A}}_{1}^{s}$ must be quadratic in $(x,y,z)$. We thus set the
numerator of each component of $\boldsymbol{\hat{A}}_{1}^{s}$ to a general
quadratic $ax^{2}+by^{2}+cz^{2}+dxy+exz+fyz$, and take the Laplacian of this
general $\boldsymbol{\hat{A}}_{s}^{2}$. Requiring that
$\boldsymbol{\hat{A}}_{1}^{s}$ solves the Poisson equation determines the
values of the constants in the quadratic. The result is

$\boldsymbol{\hat{A}}_{1}^{s}(\boldsymbol{x})=\frac{z}{16\pi^{2}|\boldsymbol{x}|^{6}}\begin{pmatrix}x\\\
y\\\ -z/2\\\ \end{pmatrix}$ (30)

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 (30) can each be written as the product of two functions

$\boldsymbol{\hat{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},$ (31)

where $\boldsymbol{h}=\boldsymbol{x}/4\pi|\boldsymbol{x}|^{3}$ and it’s
subscript represents the corresponding vector component. We guess that
$\boldsymbol{\hat{A}}_{1}^{p}$ takes a similar symmetrized form

$\boldsymbol{\hat{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}.$
(32)

Direct calculation of $\nabla^{2}\boldsymbol{\hat{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 (13) can be
integrated analytically and verified. We have not found an analytic form for
the integral on the RHS of equation (14). However, the scaling behavior of the
fields lets us deduce that the integral is zero. Since
$\boldsymbol{\hat{A}}_{1}$ decays as $|\boldsymbol{x}|^{-4}$,
$\hat{\theta}_{1}$ decays as $|\boldsymbol{x}|^{-5}$. $\hat{\psi}_{0}$ decays
as $|\boldsymbol{x}|^{-1}$, so the second term inside the integral also decays
as $|\boldsymbol{x}|^{-5}$. The RHS 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 2d 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 dynamics,
there are three self-advection velocities: $\hat{\boldsymbol{u}}_{0}(0)$ the
QG velocity at a vortex center induced by its own vorticity distribution;
$\hat{\boldsymbol{u}}_{1}^{s}(0)$, the AG velocity at a vortex center 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 center 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
3d infinite-domain Green’s function, written as

$G(\boldsymbol{x},\boldsymbol{x}^{\prime})=\frac{-1}{4\pi|\boldsymbol{x}-\boldsymbol{x}^{\prime}|},$
(33)

and the Green’s function solution of the Poisson equation,

$\displaystyle\nabla^{2}\boldsymbol{A}_{n}(\boldsymbol{x})=\boldsymbol{S}_{n}(\boldsymbol{x})\,\Longleftrightarrow\,\boldsymbol{A}_{n}(\boldsymbol{x})=\int
G(\boldsymbol{x},\boldsymbol{x}^{\prime})\boldsymbol{S}_{n}(\boldsymbol{x}^{\prime})dV^{\prime}.$
(34)

The velocity is then

$\boldsymbol{u}_{n}(\boldsymbol{x})=\nabla\times\boldsymbol{A}_{n}(\boldsymbol{x})=\int\nabla_{\boldsymbol{x}}G(\boldsymbol{x},\boldsymbol{x}^{\prime})\times\boldsymbol{S}_{n}(\boldsymbol{x}^{\prime})d^{3}\boldsymbol{x}^{\prime},$
(35)

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

$\hat{\boldsymbol{u}}_{0}(0)=\lim_{r_{0}\to
0}\int[\nabla_{\boldsymbol{x}}G(\boldsymbol{x},\boldsymbol{x}^{\prime})]_{\boldsymbol{x}=0}\times\boldsymbol{\hat{S}}_{0}(\boldsymbol{x}^{\prime})d^{3}\boldsymbol{x}^{\prime}.$
(36)

The symmetries in the vortex shape and the Green’s function render the
argument of the integral in equation (36) 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 $\boldsymbol{\hat{S}}_{1}^{s}(\boldsymbol{x}^{\prime})$. The
derivatives in equation (7), allow one to determine the symmetry of the
components of $\boldsymbol{\hat{S}}_{1}^{s}(\boldsymbol{x}^{\prime})$ 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}\neq 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

$\displaystyle\hat{q}^{r_{0}}$
$\displaystyle=\frac{e^{-r^{2}/2r_{0}^{2}}}{(2\pi r_{0}^{2})^{3/2}},$ (37)
$\displaystyle\hat{\psi}_{0}^{r_{0}}$
$\displaystyle=G(r)\hbox{Erf}(r/(\sqrt{2}r_{0})),$ (38)

where $r=|\boldsymbol{x}|$,
$\hbox{Erf}(s)=\frac{2}{\sqrt{\pi}}\int_{0}^{s}e^{-t^{2}}\,dt$ 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 center, 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 $\boldsymbol{\hat{S}}_{1}^{p}$:

$\boldsymbol{\hat{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}))$
(39)

Expanding $\boldsymbol{\hat{S}}_{1}^{p}$ in a Taylor series near the origin,
$r\ll r_{0}$, gives $\boldsymbol{\hat{S}}_{1}^{p}\sim O(r^{0})$. This then
implies that $\boldsymbol{\hat{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 2d 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:

$\displaystyle\frac{d\boldsymbol{x}_{i}}{dt}$
$\displaystyle=\boldsymbol{u}(\boldsymbol{x}_{i}),$ (40)
$\displaystyle\boldsymbol{u}(\boldsymbol{x}_{i})$
$\displaystyle=\sum_{j=1}^{N}{}^{{}^{\prime}}\left(\Gamma_{j}\hat{\boldsymbol{u}}_{0}(\boldsymbol{x}_{i}-\boldsymbol{x}_{j})+\mbox{{Ro}}\
\Gamma_{j}^{2}\hat{\boldsymbol{u}}_{1}^{s}(\boldsymbol{x}_{i}-\boldsymbol{x}_{j})\right)$
$\displaystyle\qquad+\mbox{{Ro}}\sum_{\begin{subarray}{c}j=1\\\
k=j+1\end{subarray}}^{N}{}^{{}^{\prime}}\Gamma_{j}\Gamma_{k}\hat{\boldsymbol{u}}_{1}^{p}(\boldsymbol{x}_{i}-\boldsymbol{x}_{j},\boldsymbol{x}_{i}-\boldsymbol{x}_{k})$
(41)

where $\sum{}^{{}^{\prime}}$ denotes skipping terms with $j=i$ and $k=i$. The
single-vortex and pair-vortex velocities are

$\displaystyle\hat{\boldsymbol{u}}_{0}(\boldsymbol{x})$
$\displaystyle=\frac{1}{4\pi|\boldsymbol{x}|^{3}}\begin{pmatrix}-y\\\ x\\\
\,0\end{pmatrix},$ (42)
$\displaystyle\hat{\boldsymbol{u}}_{1}^{s}(\boldsymbol{x})$
$\displaystyle=\frac{x^{2}+y^{2}-8z^{2}}{16\pi^{2}|\boldsymbol{x}|^{8}}\begin{pmatrix}-y\\\
x\\\ \,0\end{pmatrix},$ (43)
$\displaystyle\hat{\boldsymbol{u}}_{1}^{p}(\boldsymbol{x}_{1},\boldsymbol{x}_{2})$
$\displaystyle=\frac{\tilde{\boldsymbol{u}}}{16\pi^{2}|\boldsymbol{x}_{1}|^{5}|\boldsymbol{x}_{2}|^{5}},$
(44)

where

$\tilde{\boldsymbol{u}}=\begin{pmatrix}3|\boldsymbol{x}_{1}|^{2}(y_{1}z_{2}^{2}+2y_{2}z_{1}z_{2})+3|\boldsymbol{x}_{2}|^{2}(y_{2}z_{1}^{2}+2y_{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}+2x_{2}z_{1}z_{2})-3|\boldsymbol{x}_{2}|^{2}(x_{2}z_{1}^{2}+2x_{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}$
(45)

We note that, as is the case for 2d 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 $\Gamma/4\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 $\Gamma$ 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 $\Gamma$ is its
circulation. For a collection of point vortices with similar circulations,
$\Gamma_{i}\sim O(\Gamma)$, there is thus a minimum distance between QG+1
vortices and between passives and vortices, an“asymptotic horizon”, if you
will, $r_{a}\sim(\Gamma/4\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 $\Gamma$’s,
the asymptotic horizon becomes more complex.

##  5 Properties of QG+1 point vortex dynamics

###  5.1 vortex pairs

The qualitative behavior of a system composed of only two vortices, $N=2$ can
be understood from the dynamical equation, Eqs. (40)-(45). We begin by noting
that for two vortices, it is not possible to have $i$, $j$, $k$ which satisfy
$i\neq j$, $i\neq k$, and $k>j$. As a result, the second sum in equation (41)
has no terms and the pair-vortex velocity is zero.

Both 2d 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.
$\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 2d 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 $\Gamma^{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 relative magnitudes of the AG and QG vortex speeds,
$|\boldsymbol{u}_{1}^{s}|/|\boldsymbol{u}_{0}|$ scales as
$\mbox{{Ro}}\Gamma/{|\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 2d and QG vortices with $\Gamma_{1}=-\Gamma_{2}$, called hetons
in the oceanographic context, propagate along straight lines. Due to the
dependence of the AG velocity on $\Gamma^{2}$, the two vortices no longer have
the same horizontal velocity and they travel in circles, similar to opposite-
sign unequal magnitude vortices in 2d and QG, with a radius of curvature that
scales as $1/\mbox{{Ro}}$ (figure 1).

Figure 1: QG (dashed) and QG+1 (solid) heton trajectories viewed from above
for initial vortex positions $\pm(1,0,1)$ with $\Gamma=\mp 4\pi$ and
$\mbox{{Ro}}=0.2$. The height and separations of the vortices comprising QG
and a QG+1 heton are constant. The main figure shows the horizontal projection
of the trajectories for $0\leq t\leq 500$ and the inset shows the projection
for the QG+1 heton for $0\leq t\leq 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 nonzero
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 $\Gamma_{i_{image}}=\Gamma_{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 numerical solutions

In both 2d 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 $\mbox{{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 three-dimensional 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 behavior.

Figure 2: Vortex and passive particle trajectories for a pair of same-sign
vortices starting at $\pm(1.5,0,1.5)$ with $\Gamma_{1}=\Gamma_{2}=4\pi$, a
passive particle starting at $(-3.0,1.0,-1.8)$, $\mbox{{Ro}}=0.2$, and $0\leq
t\leq 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\leq t\leq 10,000$ and the time span in panels b) and
c) being $0\leq t\leq 100,000$.

We consider a tilted same-sign vortex pair with equal circulations
$\Gamma_{1}=\Gamma_{2}=4\pi$ at locations $\boldsymbol{x}=\pm(1.5,0,1.5)$ and
with $\mbox{{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$.

Figs. 2-4 show trajectories from three different passive particle initial
conditions chosen to display a variety of behaviors. 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_{Co-Rotating}<0$ ($y_{Co-Rotating}>0$), and the height extrema
occur as the trajectory crosses the $y_{Co-Rotating}=0$-axis.

Figure 5: Vertical excursion, $\hbox{Max}(z(t))-\hbox{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 $\Gamma_{1}=\Gamma_{2}=4\pi$ and $\mbox{{Ro}}=0.2$.. The black diamond is
the location of the line of passive initial conditions shown in figure 6
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
$\Gamma_{1}=\Gamma_{2}=4\pi$, and $\mbox{{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.

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.

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 three-dimensional $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 timescales, and
those with very small vertical excursions and no long timescales.

###  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 behavior. The
three vortices have identical circulations, $\Gamma_{i}=4\pi$ and are equally
distributed in the vertical coordinate. The upper and lower vortices’ initial
position 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
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 $\mbox{{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 disply small vertical oscillations with periods similar
to the horizontal rotation timescale. (figure 8).

Figure 7: Trajectories of a system of three QG point vortices,
$\mbox{{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, $\Gamma_{i}=4\pi$. a) Horizontal
projection of the trajectory of the upper vortex for $0\leq t\leq 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 $\mbox{{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, $\Gamma_{i}=4\pi$. a) Horizontal projection of the
trajectory of the upper vortex for $0\leq t\leq 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. Figure 9: Same as
figure 7 but with $y_{0}=3$. Figure 10: Same as figure 8 but with $y_{0}=3$.

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 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 timescale. As the vortex heights slowly
change, the shorter timescale horizontal motions change their character slowly
but remain qualitatively similar.

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 (panel a) and period (panel 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 $\Gamma_{i}=4\pi$ and $\mbox{{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

Figure 12: Trajectories of a system of 4 vortices with $\Gamma_{i}=4\pi$ with
positions chosen at random. a) Horizontal motion of one vortex under QG
dynamics, $\mbox{{Ro}}=0$, for $0\leq t\leq 1500$. b) Horizontal motion of the
same vortex under QG+1 dynamics starting from the same initial condition with
$\mbox{{Ro}}=0.2$. c) Vertical motion of this QG+1 vortex.

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 2d and QG point vortices, QG+1 point vortices do have
chaotic motion. Four-vortex configurations can be found with chaotic
trajectories with nonzero 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 timescales.

##  7 Discussion

We have developed a new class of point vortex dynamics, three dimensional
ageostrophic 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 ageostrophic 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 dynamics will prove useful in understanding
ageostrophic vertical transport and mixing across many geophysical and
astrophysical settings.

It is well known that 2d and QG point vortex dynamics have Hamiltonian
structures which provide insight into the dynamical behavior of the vortices.
This fact, together with the behavior of the trajectories seen here, leads us
to think it likely that QG+1 point vortex dynamics is Hamiltonian. There is a
large literature on Hamiltonian fluid dynamics, e.g., Morrison (1998). 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, 1996) 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 Jean-Luc Thiffeault, Philip J. Morrison, and Charles 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 authors report no conflict of interest.

[Data availability statement]Mathematica code to calculate the far-field
equations in section 4 is available at
https://github.com/JeffreyWeiss/QGplus1PointVortices.

[Author ORCID]J.B. Weiss, https://orcid.org/0000-0002-0706-861X

## References

  * Allen (1993) Allen, J S 1993 Iterated Geostrophic Intermediate Models. Journal of Physical Oceanography 23 (11), 2447–2461.
  * Aref (2007) Aref, Hassan 2007 Point vortex dynamics: A classical mathematics playground. Journal of Mathematical Physics 48 (6), 065401.
  * Aref & Stremler (2001) Aref, Hassan & Stremler, Mark A 2001 Point Vortex Models and the Dynamics of Strong Vortices in the Atmosphere and Oceans. In Fluid Mechanics and the Environment: Dynamical Approaches (ed. John L Lumley), pp. 1–17. Berlin, Heidelberg: Springer Berlin Heidelberg.
  * Carton (2001) Carton, Xavier 2001 Hydrodynamical Modeling Of Oceanic Vortices. Surveys in Geophysics 22 (3), 179–263.
  * Dritschel et al. (2004) Dritschel, David G., Reinaud, Jean N. & McKiver, William J. 2004 The quasi-geostrophic ellipsoidal vortex model. Journal of Fluid Mechanics 505, 201–223.
  * Dritschel et al. (1999) Dritschel, David G, de la Torre Juárez, Manuel & Ambaum, Maarten H P 1999 The three-dimensional vortical nature of atmospheric and oceanic turbulent flows. Physics of Fluids 11 (6), 1512–1520.
  * G. Abrahamyan (2020) G. Abrahamyan, Martin 2020 Vortices in Rotating and Gravitating Gas Disk and in a Protoplanetary Disk. In Vortex Dynamics Theories and Applications (ed. Zambri Harun). IntechOpen.
  * Gryanik (1991) Gryanik, V.M. 1991 Dynamics of singular geostrophic vortices near critical points of currents in a N-layer model of the atmosphere (ocean) . Izvestiya atmospheric and oceanic physics 27, 517–526.
  * Gryanik (1983) Gryanik, V. M. 1983 Dynamics of localized vortex perturbations ”vortex charges” in a baroclinic fluid. Izvestiya atmospheric and oceanic physics 19, 347–352.
  * Gryanik et al. (2000) Gryanik, Vladimir M, Doronina, Tatiana N, Olbers, Dirk J & Warncke, Torsten H 2000 The theory of three-dimensional hetons and vortex-dominated spreading in localized turbulent convection in a fast rotating stratified fluid. Journal of Fluid Mechanics 423, 71–125.
  * Gryanik et al. (2006) Gryanik, V. M., Sokolovskiy, M. A. & Verron, J. 2006 Dynamics of heton-like vortices. Regular and Chaotic Dynamics 11 (3), 383–434.
  * Helmholtz (1867) Helmholtz, H 1867 LXIII. On Integrals of the hydrodynamical equations, which express vortex-motion. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 33 (226), 485–512.
  * Hirt et al. (2018) Hirt, Mirjam, Schielicke, Lisa, Müller, Annette & Névir, Peter 2018 Statistics and dynamics of blockings with a point vortex model. Tellus A: Dynamic Meteorology and Oceanography 70 (1), 1–20.
  * Holm (1996) Holm, Darryl D 1996 Hamiltonian balance equations. Physica D: Nonlinear Phenomena 98 (2), 379–414.
  * Kida (1981) Kida, Shigeo 1981 Motion of an Elliptic Vortex in a Uniform Shear Flow. Journal of the Physical Society of Japan 50 (10), 3517–3520.
  * Kirchhoff (1876) Kirchhoff, G. 1876 Vorlesungen über matematische Physik. Mechanik. Leipzig: Teubner.
  * Lucarini et al. (2014) Lucarini, Valerio, Blender, Richard, Herbert, Corentin, Ragone, Francesco, Pascale, Salvatore & Wouters, Jeroen 2014 Mathematical and physical ideas for climate science. Reviews of Geophysics 52 (4), 809–859.
  * Mahdinia et al. (2017) Mahdinia, Mani, Hassanzadeh, Pedram, Marcus, Philip S & Jiang, Chung-Hsiang 2017 Stability of three-dimensional Gaussian vortices in an unbounded, rotating, vertically stratified, Boussinesq flow: linear analysis. Journal of Fluid Mechanics 824, 97–134.
  * McKiver (2020) McKiver, William J 2020 Balanced ellipsoidal vortex at finite Rossby number. Geophysical & Astrophysical Fluid Dynamics 114 (4-5), 453–480.
  * McKiver & Dritschel (2016) McKiver, William J & Dritschel, David G 2016 Balanced solutions for an ellipsoidal vortex in a rotating stratified flow. Journal of Fluid Mechanics 802, 333–358.
  * McWilliams (1984) McWilliams, James C. 1984 The emergence of isolated coherent vortices in turbulent flow. Journal of Fluid Mechanics 146, 21–43.
  * McWilliams & Gent (1980) McWilliams, James C. & Gent, Peter R. 1980 Intermediate Models of Planetary Circulations in the Atmosphere and Ocean. Journal of the Atmospheric Sciences 37 (8), 1657–1678.
  * McWilliams & Weiss (1994) McWilliams, James C & Weiss, Jeffrey B 1994 Anisotropic geophysical vortices. Chaos: An Interdisciplinary Journal of Nonlinear Science 4 (2), 305–311.
  * McWilliams et al. (1994) McWilliams, James C, Weiss, Jeffrey B & Yavneh, Irad 1994 Anisotropy and Coherent Vortex Structures in Planetary Turbulence. Science 264 (5157), 410.
  * Meacham (1992) Meacham, S. P. 1992 Quasigeostrophic, ellipsoidal vortices in a stratified fluid. Dynamics of Atmospheres and Oceans 16 (3-4), 189–223.
  * Meacham et al. (1998) Meacham, S. P., Morrison, P. J. & Flierl, G. R. 1998 Hamiltonian moment reduction for describing vortices in shear. Physics of Fluids 9 (8), 2310.
  * Melander et al. (1986) Melander, M. V., Zabusky, N. J. & Styczek, A. S. 1986 A moment model for vortex interactions of the two-dimensional Euler equations. Part 1. Computational validation of a Hamiltonian elliptical representation. Journal of Fluid Mechanics 167 (-1), 95\.
  * Miyazaki et al. (2001) Miyazaki, Takeshi, Furuichi, Yu & Takahashi, Naoya 2001 Quasigeostrophic Ellipsoidal Vortex Model. Journal of the Physical Society of Japan 70 (7), 1942–1953.
  * Morrison (1998) Morrison, P J 1998 Hamiltonian description of the ideal fluid. Reviews of Modern Physics 70 (2), 467–521.
  * Muraki et al. (1999) Muraki, David J., Snyder, Chris & Rotunno, Richard 1999 The Next-Order Corrections to Quasigeostrophic Theory. Journal of the Atmospheric Sciences 56 (11), 1547–1560.
  * Reinaud & Dritschel (2018) Reinaud, Jean N & Dritschel, David G 2018 The merger of geophysical vortices at finite Rossby and Froude number. Journal of Fluid Mechanics 848, 388–410.
  * Reznik & Kizner (2007) Reznik, Gregory & Kizner, Ziv 2007 Two-layer quasi-geostrophic singular vortices embedded in a regular flow. Part 1. Invariants of motion and stability of vortex pairs. Journal of Fluid Mechanics 584, 185–202.
  * Showman et al. (2010) Showman, A. P., Cho, J. Y.-K. & Menou, K. 2010 Atmospheric circulation of exoplanets. In Exoplanets (ed. S. Seager), pp. 471–516. University of Arizona Press.
  * Sokolovskiy et al. (2020) Sokolovskiy, Mikhail A, Carton, Xavier J & Filyushkin, Boris N 2020 Mathematical Modeling of Vortex Interaction Using a Three-Layer Quasigeostrophic Model. Part 1: Point-Vortex Approach. Mathematics 8 (8).
  * Tsang & Dritschel (2015) Tsang, Yue Kin & Dritschel, David G. 2015 Ellipsoidal vortices in rotating stratified fluids: Beyond the quasi-geostrophic approximation. Journal of Fluid Mechanics 762, 196–231.
  * Viúdez (2016) Viúdez, A 2016 Isolated spheroidal geophysical vortices. Geophysical & Astrophysical Fluid Dynamics 110 (1), 50–77.