1. Introduction
We introduce and analyse two new waves satisfying Euler's equations in Lagrangian coordinates – one on the 
$\beta$-plane and another in the 
$\gamma$ approximation. These flows are planetary vortex waves and exhibit behaviour that is analogous to, but distinct from, motions on the plane, due to variable rotation. The variable rotation induces a Lagrangian mean flow for these waves. For particular configurations in the 
$\gamma$ approximation, these waves have a jet-like mean flow and describe polar jet streams. Exact solutions in the Lagrangian frame allow us to examine highly nonlinear properties of these waves and connect them to dynamically relevant quantities, such as the vorticity.
Motion on a rotating surface with non-zero curvature is distinct from its planar counterpart due to variations in the planetary vorticity that support wave motions. These waves represent oscillations in the amount of planetary and relative vorticity about lines of constant latitude. The Coriolis parameter describing the planetary vorticity is considered in a hierarchy of approximations by Taylor expanding it about a point of latitude: the lowest order of approximation is the 
$f$-plane, the first-order approximation is the 
$\beta$-plane, and the second-order approximation is denoted the 
$\gamma$ approximation. Interestingly, the 
$\beta$-plane at the equator corresponds to an embedding space with a metric that has zero curvature. That is, curvature effects are necessary for a variable Coriolis parameter, but have no direct impact on the geometry of the embedding space. The 
$\gamma$ approximation describes motion near the poles of a rotating sphere, where 
$\beta$ vanishes. Solutions to these equations exist on a plane with a variable rotation rate that again supports wave motion.
A remarkable set of exact solutions to Euler's equations on the plane was found by Abrashkin & Yakubovich (Reference Abrashkin and Yakubovich1984). This class of solutions includes Gerstner's trochoidal wave, the Kirchhoff ellipse, and a class of solutions known as Ptolemaic, or polygonal, vortices, which are a rotationally symmetric analogue of Gerstner waves (Guimbard & Leblanc Reference Guimbard and Leblanc2006). For 
$z=x+\mathrm {i} y$, where 
$(x,y)$ are Cartesian coordinates, and 
$\mathrm {i}^2 =-1$, the Abrashkin & Yakubovich (Reference Abrashkin and Yakubovich1984) solutions are of the form 
$z=f(s){\rm e}^{{\mathrm {i} \omega _1 t}}+g(\bar {s}){\rm e}^{{\mathrm {i} \omega _2 t}}$, where 
$s=a+\mathrm {i} b$ for 
$(a,b)$ Lagrangian particle labels, 
$\omega _{1}$ and 
$\omega _{2}$ are constants, and the overline denotes complex conjugate. The functions 
$f,g$ are harmonic, but this is the only restriction on their form. When viewed in a certain reference frame, both Gerstner and Ptolemaic waves may be written as permanent progressive solutions of the form 
$z=z(a-\omega t,b)$.
A warning – exact solutions in one frame generally cannot be written in closed form in a different frame. An example of this is the Gerstner wave (Lamb Reference Lamb1932), which has a relatively simple form in the Lagrangian frame but cannot be written in closed form in the Eulerian frame, as this reduces to solving Kepler's equation: one would have to find 
$a$ as a function of 
$x$ in closed form in the equation 
$x=a- \sin a$, which is not possible. In two dimensions, vorticity is conserved along fluid particles so there is a close connection between the Lagrangian mean flow and the vorticity. This mass flux Doppler shifts the wave frequency, and for certain vorticity distributions the particle trajectories take relatively simple forms.
On the 
$\beta$-plane at the equator and in the 
$\gamma$ approximation, we find a direct analogue of Gerstner's wave, but now the frequencies must depend on the particle label 
$b$ in order to satisfy the vorticity equation. For finite-amplitude waves, a nonlinear correction to the wave phase speed arises, and this can be interpreted as being due to a Lagrangian mean velocity, which leads to a departure in behaviour from their planar counterparts. These waves exist only over a region of the embedding space, as the vorticity contours always form cusps at some critical value of the parameters. In the plane, this is handled by taking the vorticity to be constant outside a critical bounding vorticity contour. For example, Gerstner's wave takes the flow to be irrotational outside a bounding vorticity contour that is specified to be the free surface. On surfaces approximating the rotating sphere, we take the ambient flow immediately outside the waves to have constant potential vorticity, then take the flow to be zonal in the far field. This yields solutions that are well defined over the entire embedding space.
Taking our solutions as initial conditions in the fully spherical rotating Navier–Stokes equations, using the scheme developed in Salmon & Pizzo (Reference Salmon and Pizzo2023), we examine the temporal evolution of these solutions. Unlike on the embedding spaces, we now connect the solutions in the 
$\gamma$ approximation to those on the sphere by choosing the zonal flow in the outer region to be one that exactly solves the equations of motion on the rotating sphere. The numerical experiments are a strong test of these solutions and yield solutions that range from permanent progressive waves to waves that rapidly overturn and break, generating turbulent motion.
The plan of the paper is as follows. In § 2 and Appendix A, the equations of motion are presented. In § 3, new exact waves are discussed. These solutions are used as initial conditions for the fully spherical rotating equations of motion in § 4. The connection between these flows and jet streams is discussed in § 5, after which the conclusions of the paper are given.
2. The 
$\beta$ and $\gamma$ approximations
We examine two-dimensional incompressible flow with coordinates 
$(x,y)$ obeying conservation of mass (a derivation of these equations from the full equations on a rotating sphere is presented in Appendix A):
\begin{equation} \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} = 0, \end{equation}
for 
$(u,v)$ the velocity, and conservation of potential vorticity
\begin{equation} \frac{\partial q}{\partial t}+u\,\frac{\partial q}{\partial x}+v\, \frac{\partial q}{\partial y}=\frac{{\rm D} q}{{\rm D} t}=0, \end{equation}where the potential vorticity is given by
\begin{equation} q=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}+F(x,y), \end{equation}
and 
${\rm D}/{\rm D} t =\partial _t +\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }$.
We consider the dynamics in two locations. First, near the equator 
$F = \beta y$, with the Coriolis constant being zero at this location, corresponds to a 
$\beta$-plane approximation. Similarly, we analyse flow near the poles, where 
$F = f_0-\gamma (x^2+y^2)$ corresponds to the 
$\gamma$ approximation. Following Phillips (Reference Phillips1973), latitude-dependent scale factors do not arise in these classical expansions in the 
$\beta$-plane because we are at the equator, and in the 
$\gamma$ approximation because we additionally assume that 
$\gamma \gg 1$. For a more thorough discussion, see Dellar (Reference Dellar2011) and the derivation presented in § A.4. Throughout this paper, we repeatedly refer to ‘the’ 
$\beta$-plane approximation (and ‘the’ 
$\gamma$ approximation), with the understanding that this is one of many such approximations that bears this name. Here, 
$\beta,\gamma,f_0$ are constants, and we take the deformation length scale to be infinite, so there is no vorticity associated with vertical stretching.
The main focus of this paper is on solutions in the 
$\gamma$ approximation. Both of the approximations contain variable rotation and make a simplifying assumption about the curvature of the embedding surface. These approximations require non-zero surface curvature to possess variable rotation rates, but the curvature of the geometry alone is not present in both approximations.
Consider fluid particles labelled by coordinates 
$\boldsymbol {a}=(a,b)$ and rewrite the Eulerian governing equations in the Lagrangian reference frame 
$(x(a,b,\tau ),y(a,b,\tau ))$, where 
$\tau$ is time and we emphasize that 
$\partial _{\tau }$ with 
$(a,b)$ fixed is the time derivative following a fluid particle so that 
$\partial _{\tau }={\rm D}/{\rm D} t$.
We can assign labelling coordinates so that 
${\rm d}(\text {mass})={\rm d} a \,{\rm d} b$. Defining the reciprocal of the non-dimensionalized mass density of the fluid 
$\rho$ as 
$J=J(a,b)$, we have
\begin{equation} \frac{1}{\rho}\equiv J(a,b) =\frac{\partial (x,y)}{\partial(a,b)}, \end{equation}where
\begin{equation} \frac{\partial (x,y)}{\partial(a,b)}\equiv x_ay_b-x_by_a. \end{equation}(For a more thorough discussion of this map, see Salmon Reference Salmon1998.) The labels follow the fluid, from which it follows that
\begin{align} \frac{\partial J}{\partial \tau} &= \frac{\partial(x_{\tau},y)}{\partial (a,b)}+ \frac{\partial(x,y_{\tau})}{\partial(a,b)} = \left(\frac{\partial(x_{\tau},y)}{\partial(x,y)}+ \frac{\partial(x,y_{\tau})}{\partial(x,y)}\right)\frac{\partial (x,y)}{\partial(a,b)}\nonumber\\ & =\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right) \frac{\partial (x,y)}{\partial(a,b)}=0, \end{align}
where the term in parentheses vanishes due to the incompressibility of the fluid, i.e. (3.1a,b). Therefore, 
$J(a,b)$ is a time-independent function. Additionally, we require that 
$J(a,b)$ does not change sign so that the mass density 
$\rho$ stays finite.
Equation (2.2) requires the potential vorticity to be conserved on fluid particles. The relative vorticity in the Lagrangian frame is
\begin{equation} v_x-u_y=\frac{\partial (v,y)}{\partial(x,y)}+\frac{\partial (u,x)}{\partial(x,y)}= \frac{\partial (a,b)}{\partial(x,y)}\left(\frac{\partial (x_{\tau},x)}{\partial(a,b)}+ \frac{\partial (y_{\tau},y)}{\partial(a,b)}\right), \end{equation}so (2.3) becomes
\begin{equation} q(a,b)\,J(a,b) = \frac{\partial(x_{\tau},x)}{\partial(a,b)}+ \frac{\partial(y_{\tau},y)}{\partial(a,b)}+J(a,b)\,F(x,y). \end{equation}
Our task is to find 
$(x,y)$ that satisfy (2.4) and (2.8) subject to the constraints that 
$J(a,b)$ and 
$q(a,b)$ are time-independent, and 
$J(a,b)$ does not change sign.
Although the fluid acceleration is simpler in the Lagrangian frame, the above analysis shows that the nonlinear mapping between frames can turn linear operators, such as a partial derivative, into nonlinear operators. This leads to coupled partial differential equations that are different in form from their Eulerian counterparts, therefore we must employ different strategies to find classes of exact solutions. Although the Lagrangian reference frame is considered here, there is the possibility of finding exact solutions in a frame that is neither purely Eulerian nor Lagrangian, but a mixture of the two (see the formulation given by Virasoro Reference Virasoro1981).
To better understand flow on the 
$\beta$-plane and in the 
$\gamma$ approximation in the Lagrangian reference frame, we present zonal flow solutions to these equations. These flows prove useful for patching together solutions over the entire embedding space, as is discussed in § 3.3. In the 
$\beta$-plane, zonal flow takes the form
\begin{equation} x= a +U(b)\,\tau,\quad y = b. \end{equation}
Here, 
$J(a,b)=1$ and 
$q(a,b) =U'+\beta b$, where 
$U$ is the unspecified mean flow, and a prime represents a derivative with respect to 
$b$. Contours of constant potential vorticity are horizontal lines.
In the 
$\gamma$ approximation, zonal flow takes the form
\begin{equation} x = A_0\sin (ka-(\omega+\varOmega(b))\tau )\,{\rm e}^{kb},\quad y = A_0\cos (ka-(\omega+\varOmega(b))\tau )\,{\rm e}^{kb}, \end{equation}
where 
$A_0$, 
$\omega$ and 
$k$ are constants. These particle locations imply
\begin{equation} J =A_0^2k^2 \,{\rm e}^{2kb},\quad q =f_0-\frac{\gamma}{k^2}\, J+2(\omega+\varOmega)+\frac{\varOmega'}{k}. \end{equation}The velocity is entirely in the azimuthal direction, and contours of constant potential vorticity take the form of circles. Zonal flow solutions for the rotating sphere are presented in Appendix A.
We can also write down Rossby wave solutions in the Lagrangian frame for the 
$\beta$-plane approximation:
\begin{equation} x=a,\quad y = b+\epsilon \cos (ka-\omega \tau), \end{equation}
where 
$\epsilon, k,\omega$ are yet to be specified constants. These expansions imply 
$J(a,b)=1$, while
\begin{equation} q(a,b) = \beta (b+\epsilon \cos (ka-\omega \tau))+\epsilon \omega k \cos (k a-\omega \tau). \end{equation}
For 
$q(a,b)$ to be 
$\tau$-independent, we require
\begin{equation} \omega ={-}\frac{\beta}{k}, \end{equation}in agreement with the classical result for the Rossby wave dispersion relationship in the Eulerian frame.
In the 
$\gamma$ approximation, the situation is complicated by the geometry of the domain. In the Eulerian frame, Leblond (Reference Leblond1964) and Nof (Reference Nof1990) show that the vorticity contours have radial dependence in the form of Bessel functions, while the dispersion relationship is a function of roots of these Bessel functions. As far as we can tell, there is no simple Lagrangian analogue here, an illustration that there are solutions in the Eulerian frame that cannot be simply written in the Lagrangian frame. The situation is even more complicated for Rossby–Haurwitz waves on the sphere.
3. Exact solutions
We now present a description of waves on the 
$\beta$-plane and in the 
$\gamma$ approximation. They cannot cover the entire embedding space as the function 
$J(a,b)$ would then change sign and the vorticity contours would form cusps. Instead, we write down the particle locations of these waves over some limited region of the embedding space. We then describe ambient conditions outside these areas, where we specify the flow to have constant potential vorticity over a limited region, before being smoothly connected with a region of zonal flow. We consider each of these regions of the flow in turn.
3.1. Waves on the $\beta$-plane
Our strategy is to take known planar solutions from Abrashkin & Yakubovich (Reference Abrashkin and Yakubovich1984) and to allow the frequency to have label dependence. We first look for waves in the 
$\beta$-plane to build up intuition on how to tackle the more complicated 
$\gamma$ approximation. The analogue of the Gerstner wave on the 
$\beta$-plane is given by
\begin{equation} x= a+ U(b)\,\tau - \epsilon \sin \theta \,{\rm e}^{kb},\quad y=b+\epsilon \cos \theta \,{\rm e}^{kb}, \end{equation}where
\begin{equation} \theta = ka - ( \omega-k\,U(b) ) \tau, \end{equation}
and the Lagrangian mean flow 
$U$ is specified by
\begin{equation} U(b)=\frac{\beta}{2k^2}\left(\epsilon^2k^2\,\frac{{\rm e}^{2kb}}{2}-k b\right). \end{equation}
Here, 
$\epsilon$ is the amplitude of the wave, 
$k,\omega$ are the (constant) wavenumber and frequency, and 
$b\le 0$. We note that the mean velocity goes to 
$\infty$ as 
$b\to - \infty$, so that these flows, analogous to constant shear flow in an infinite fluid in the plane, are not physical over the entire embedding space.
These waves imply
\begin{equation} J = 1-\epsilon^2k^2 \,{\rm e}^{2kb} \end{equation}and
\begin{equation} qJ = 2\omega \epsilon^2k^2 \,{\rm e}^{2kb}+\beta \left(\frac{1}{2k}-\epsilon^4 k^3 \,{\rm e}^{4kb}+b\right). \end{equation}
$J$ does not change sign as long as 
$\epsilon k< 1$, which also corresponds to lines of constant potential vorticity remaining single-valued and differentiable. In the limit that 
$\beta =0$, we return to the Gerstner wave solution in the plane. The particle trajectories for these waves are no longer perfect circles as they are in the plane, and instead are trochoids with a distinct shape from the trochoidal contours of constant vorticity. Note that in the plane, gravity acts as a restoring force for the waves that oscillate vertically and attenuate in amplitude away from the free surface. On the 
$\beta$-plane, planetary vorticity is the restoring force, with the waves oscillating about lines of fixed latitude and decaying far away from the equator.
For 
$b\ll 0$, we have 
$J\to 1$, 
$U\to -\beta b/2k$ and 
$q \to \beta (1/2k +b)$ so that 
$U$ and 
$q$ are linear in 
$b$ and easily related to one another far below the limiting contour (
$b=0$). The curvature 
$U''=\beta \epsilon ^2 k^2 \,{\rm e}^{2kb}$ is largest near the surface then rapidly attenuates with depth as figure 1(b) shows. The geometry of lines of constant Lagrangian mean flow and vorticity are shown in figure 1 for 
$k=2$ and 
$\epsilon =1/2$.
Figure 1. (a) Geometry of lines of constant vorticity and Lagrangian mean velocity for the exact Gerstner wave solution on the 
$\beta$-plane. Here, 
$k=2$, 
$\epsilon =1/2$, 
$b\le 0$ and 
$\omega =0$. (b) The first two plots show the Lagrangian mean velocity for these waves as a function of 
$b$ and its curvature, with primes representing derivatives with respect to 
$b$. The last two plots show the vorticity and its curvature. Both curvature terms rapidly go to zero for 
$b\ll 0$.
The dispersion relationship (3.3) is distinct from Rossby waves as it has 
$b$ dependence. For 
$\epsilon$ small, 
$U\approx -\beta b/2k$, so that these waves have an inverse wavenumber dependence, like the Rossby wave, and also propagate to the left when 
$b<0$. We can then interpret these waves as being incident on a steady shear flow. Additionally, solutions can be made to cover the entire domain, decaying exponentially away from the contour 
$b=0$, by taking 
$x\to x$ and 
$y\to - y$, and choosing the parameter 
$\epsilon$ and an interval for 
$b$ so that the waves do not intersect.
3.2. Ptolemaic vortex waves in the $\gamma$ approximation
The main focus of this paper is on solutions in the 
$\gamma$ approximation near the poles, corresponding to 
$F=f_0 - \gamma (x^2+y^2)$, where in our non-dimensionalized coordinates the radius of the sphere that we are approximating is 1 (see Appendix A). Our strategy is to once again start from a planar solution and take the rotation rate to depend on a particle label. Planetary Ptolemaic vortex waves, with the nomenclature based on analogous solutions in the plane (Abrashkin & Yakubovich Reference Abrashkin and Yakubovich1984; Abrashkin, Zenkovich & Yakubovich Reference Abrashkin, Zenkovich and Yakubovich1996; Guimbard & Leblanc Reference Guimbard and Leblanc2006) and the fact that the particle trajectories can be epitrochoids (i.e. the curve generated by tracing a point along a radial line of a circle rolling along another circle), are described by
$$\begin{gather} x/A_0={\rm e}^{kb}\sin \theta -\epsilon \,{\rm e}^{nk b}\sin n\theta, \end{gather}$$
$$\begin{gather}y/A_0={\rm e}^{kb}\cos\theta + \epsilon \,{\rm e}^{nkb}\cos n\theta, \end{gather}$$where
\begin{equation} \theta = ka-(\omega+\varOmega(b)) \tau, \end{equation}
for 
$k,\omega, A_0$ constants. We define 
$\tilde {x}=x/A_0$ and 
$\tilde {y}=y/A_0$, and drop the tildes from here onwards.
If
\begin{equation} \varOmega(b)= \frac{\gamma }{2n}\,({\rm e}^{2kb}-\epsilon^2 n{\rm e}^{2knb}), \end{equation}then
\begin{equation} J = k^2({\rm e}^{2kb}-\epsilon^2 n^2 {\rm e}^{2nkb}) \end{equation}and
\begin{align} qJ &= \frac{k^2
}{n}\,(\gamma ({-}2+n)\,{\rm
e}^{4kb}+\gamma(1+2n)\epsilon^4 n^3 {\rm e}^{4knb} \nonumber\\
&\quad -(f_0+2\omega)n{\rm e}^{2knb}-(f_0-2n\omega)n^3\epsilon^2
{\rm e}^{2knb}).
\end{align}
The term 
$\varOmega$ should be interpreted as follows. Define 
$\varphi$ implicitly as
\begin{equation} \tan \varphi = \frac{y}{x}; \end{equation}then
\begin{equation} \langle \varphi_{\tau}\rangle=\frac{1-n}{2}\,(\omega+\varOmega), \end{equation}
for 
$\langle \ \rangle \equiv (2{\rm \pi} )^{-1}\int _0^{2{\rm \pi} } \,{\rm d} a$, so that 
$\varOmega$ is related to the mean azimuthal angular rate of change. The zonal velocity corresponding to this, at this level of approximation, is given by 
$r\varphi _{\tau }$, where 
$r^2=x^2+y^2$. Therefore, even if 
$\varOmega$ is relatively weak, for planetary flows 
$r$ can be large when put in dimensional form, implying that these mean flows need not be weak. Depending on the signs of 
$n$ and 
$b$, this velocity does not necessarily diverge over the region in the embedding space over which the waves are defined, unlike the scenario in the 
$\beta$-plane.
The mean zonal velocity is
$$\begin{gather} \mathcal{U}\equiv
\langle r\varphi_{\tau}\rangle = \frac{1}{2{\rm \pi}}
\int_0^{2{\rm \pi}} \frac{{\rm
e}^{kb}({\rm e}^{{2k(1\!-\!n)b}}-n\epsilon^2 +(1\!-\!n)\epsilon
{\rm e}^{{k(1\!-\!n)b}} \cos ((1\!+\!n)\theta))(\omega +
\varOmega)}{\sqrt{{\rm e}^{{2k(1-n)b}}+ \epsilon^2 +2\epsilon
{\rm e}^{{k(1\!-\!n)b}} \cos ((1\!+\!n)\theta) }}\,{\rm d} a,
\end{gather}$$which may be integrated to give
\begin{align} \mathcal{U}
&=\frac{{\rm
e}^{knb}}{\rm \pi}\,\bigl((1-n)({\rm e}^{{(1-n)kb}}+\epsilon)\,E(m) +(1+n)({\rm e}^{{(1-n)kb}}-\epsilon)\,K(m)\bigr)
(\omega+\varOmega),
\end{align}
where 
$K,E$ are the complete elliptical integrals of the first and second kind, respectively, and
\begin{equation} m = \frac{4 \epsilon
{\rm e}^{{(1-n)kb}}}{({\rm e}^{{(1-n)kb}}+\epsilon)^2}.
\end{equation}
The exact solutions allow us to write down properties of these waves, such as their mean zonal velocity, even in highly nonlinear scenarios when 
$\epsilon$ is not much less than 1. When 
$n>0$, we take the label 
$b$ to be negative so that the waves are described for 
$b\in (-\infty,0)$, and when 
$n<0$, we take 
$b\in (0,\infty )$.
The parameters 
$n$ and 
$\epsilon$ must be chosen to keep 
$J$ same signed and the contours of 
$q$ single-valued. These contours are of the form (at 
$b=t=0$)
$$\begin{gather} x=\sin \theta_0 -\epsilon \sin n\theta_0, \end{gather}$$
$$\begin{gather}y=\cos\theta_0 + \epsilon \cos n\theta_0, \end{gather}$$
for 
$\theta _0 = ka$, and have a cusp when 
$x_a=y_a=0$, or
\begin{equation} k\cos \theta_0 -\epsilon nk \cos n\theta_0={-}k\sin \theta_0 -\epsilon nk \sin n\theta_0=0. \end{equation}
A solution to this is 
$\theta _0=0$ and 
$\epsilon n=-1$, so we must have
\begin{equation} |\epsilon n|< 1. \end{equation} Figure 2 shows two Ptolemaic vortex waves. Figures 2(a,b) show the geometry of lines of constant zonal Lagrangian mean flow 
$\mathcal {U}$ and potential vorticity 
$q$. The first wave has 
$n=-4$ and 
$b\ge 0$. Figures 2(b,d) display 
$n=2$ with 
$b\le 0$. The blue line segments represent the location where 
$q$ and 
$\mathcal {U}$ are computed in the bottom row. Both waves have 
$\omega =0$ and potential vorticity and zonal mean Lagrangian flows that change monotonically.
Figure 2. Two Ptolemaic vortex waves. (a,b) The geometry of lines of constant potential vorticity and Lagrangian zonal mean flow. The blue line segments show the transect where 
$\mathcal {U}$ (in black, with values given by the left ordinate) and 
$q$ (in red, with values given by the right ordinate) are plotted in (c,d). (a,c) A wave with 
$n=-4$, 
$A_0=0.2$, 
$\epsilon =0.2$, 
$\gamma = 1$, 
$\omega =0$ and 
$b\ge 0$. The zonal Lagrangian mean flow and potential vorticity, as a function of 
$r=\sqrt {x^2+y^2}$, both increase away from the boundary of the wave. (b,d) A wave with 
$n=2$, 
$A_0=0.2$, 
$\epsilon =0.4$, 
$\gamma =1$, 
$\omega =0$ and 
$b\le 0$. The potential vorticity and Lagrangian zonal mean flow decrease in magnitude towards the boundary of the wave.
In figure 3, we take 
$n=5$, 
$b\le 0$, and choose 
$\omega$ to be non-zero. The mean flow and potential vorticity are shown in figure 3(b). Both contain a critical point. The mean flow is peaked and represents a jet-like flow. The mean flow for various values of 
$\omega$ is shown in figure 3(c). Increasing values of 
$\omega$ have critical points at increasing values of 
$r$.
Figure 3. (a) Ptolemaic wave with 
$n=5$, 
$\epsilon =0.1$, 
$A_0=0.2$, 
$\gamma =1$ and 
$b\le 0$. Geometry of contours of constant Lagrangian mean flow and potential vorticity. The blue line is where the mean flow and potential vorticity are computed in (b,c). (b) Here, 
$\omega =0.004$ is non-zero, yielding a localized maximum in the zonal Lagrangian mean flow and potential vorticity, corresponding to jet stream motion near the centre of this vortex wave. (c) The Lagrangian mean flow for 
$\omega$ varying between 
${\pm }0.016$ in equal increments, showing that the location of the maximum Lagrangian mean velocity varies as a function of 
$\omega$.
The nonlinear term in 
$\varOmega$ arises due to the variable Coriolis parameter. The amplitude squared dependence of the frequency corresponds to the existence of a four wave resonance (cf. the triad resonance of Rossby waves) and has stability implications for these waves (see Pizzo et al. (Reference Pizzo, Lenain, Rømcke, Ellingsen and Smeltzer2023) for analogous considerations of surface gravity waves). Linear stability analysis for these waves is pursued elsewhere – here, we generate solutions that may be tested numerically – as a much stronger check on the persistence of these waves.
3.3. Ambient conditions and far-field behaviour
The waves presented above cannot exist everywhere when 
$\epsilon$ is non-zero, as the function 
$J$ would change sign and the potential vorticity contours would become multi-valued. Therefore they exist over some range of label space 
$(a,b)$. This is analogous to the scenario described by Abrashkin & Yakubovich (Reference Abrashkin and Yakubovich1984) and discussed in more detail in Salmon (Reference Salmon2020). We extend the discussion presented in Salmon (Reference Salmon2020) to generate solutions that are valid over the entire 
$\beta$-plane and 
$\gamma$-plane by prescribing ambient conditions outside the waves.
In the plane, the flow may be taken to have uniform vorticity outside the bounding curve given by 
$b=\text {constant}$. A simple analogue is to assume that our planetary waves are surrounded by regions of uniform potential vorticity, which may be subsequently connected to far-field regions that have zonal flow. Physically, regions of constant potential vorticity can occur from wave breaking that leads to potential vorticity mixing and homogenization where the vorticity contours have most curvature, that is, near their bounding contours.
In both scenarios, we must specify regions where the potential vorticity 
$\nabla ^2 \psi +f_0-\gamma (x^2+y^2)$ is constant, where 
$\psi$ is the stream function associated with the relative vorticity of the flow in this yet-unspecified region evaluated near the equator (the 
$\beta$ approximation considered in this paper) or near the poles (the 
$\gamma$ approximation considered in this paper). The potential vorticity can be discontinuous across this boundary, but the velocity must remain continuous. This sets the boundary condition on 
$\psi$: the tangent velocity of the particle locations provided in §§ 3.1 and 3.2 should be equal to 
$\psi _{\tilde {n}}$, for 
$\tilde {n}$ the normal to the contour. Then 
$\psi$ is determined at each point in time, and specifies a solution everywhere valid on the 
$\beta$- and 
$\gamma$-planes. Of central importance here is the fact that we need not specify the actual trajectories in this region, which can be complicated because of the geometry, as knowledge of 
$q$ is enough to specify the problem once the continuity of the velocity field has been ensured.
The region over which the potential vorticity is uniform must be finite for the relative vorticity to stay bounded. Subsequently, when 
$n>0$, we can connect the constant potential vorticity region to a region with the zonal flow solutions presented in § 2. The choice of 
$\omega$ and 
$\varOmega$ for the zonal flow solutions again allows the velocity field to remain continuous. Figure 4 shows a sketch of these scenarios, depending on the sign of 
$n$ in the 
$\gamma$ approximation. For 
$n<0$, the polar cap is taken to have constant potential vorticity, and outside this, the flow takes the form of a Ptolemaic wave. When 
$n>0$, the inner region is the Ptolemaic wave, while immediately outside this wave the flow has constant potential vorticity, and in the far field the flow is zonal. Figures 4(a,b) are on the 
$\gamma$-plane, while figure 4(c) shows the configuration on the southern hemisphere. For this scenario, we take the far field to be zonal flow solutions on the fully rotating sphere (see (A27)).
Figure 4. Black indicates contours of potential vorticity and Lagrangian mean flow, for (a,b) the 
$\gamma$ approximation, and (c) the stereographic coordinates on the southern hemisphere of the sphere. (a) The case for 
$n=-2$. The red dashed line partitions the two regions: inside the polar cap, the solution has uniform potential vorticity while beyond that the potential vorticity is given by (3.11). As the radius of the contours gets large, the Ptolemaic solutions tend towards zonal flow. (b) A Ptolemaic solution with 
$n=5$. There are three regions here: the innermost region has potential vorticity given by (3.11); the region between the two dashed red lines has uniform potential vorticity; the outer region is zonal flow in the 
$\gamma$ approximation. (c) The same as (b), but now we take the flow in the outer region to describe zonal flow on the rotating sphere, as given by (A27). The entire southern hemisphere is shown in stereographic coordinates 
$(\xi,\eta )$.
4. Numerics
Salmon & Pizzo (Reference Salmon and Pizzo2023) presented a numerical model of two-dimensional flow on the sphere using stereographic coordinates and a generalization of Arakawa's method developed by Salmon & Talley (Reference Salmon and Talley1989). This solver can be run inviscidly or viscously, and provides an opportunity to examine the solutions presented above. The numerical formulation is Eulerian in nature, so we must carefully provide the initial potential vorticity conditions in this frame. In addition to mass and potential vorticity conservation, the Gauss constraint requires the potential vorticity to integrate to zero over the sphere.
By construction, Arakawa's method conserves discrete analogues of energy, potential vorticity and enstrophy when the viscosity vanishes. The equations of motion are solved in two sets of stereographic coordinates in the unit disc, one corresponding to the northern hemisphere, and the other to the southern hemisphere. At the equator, we require the solutions to be continuous. These discs are covered in quadrilateral elements, and the nodal values of the stream function and potential vorticity serve as the dependent variables of the model.
Viscosity on non-Euclidean surfaces has been a source of confusion as various reasonable sounding constraints lead to different forms of the final form of the viscosity. Following the recommendation of Gilbert, Riedinger & Thuburn (Reference Gilbert, Riedinger and Thuburn2014), we use a viscosity that conserves angular momentum, which takes a particularly simple form in stereographic coordinates. There is an error in the form of the viscosity presented in Salmon & Pizzo (Reference Salmon and Pizzo2023), which is corrected in Salmon & Pizzo (Reference Salmon and Pizzo2024).
We integrate the waves presented in § 3. As discussed in the previous section, for 
$n<0$, we require a constant potential vorticity region at the polar cap surrounded by the Ptolemaic waves. For 
$n>0$, we prescribe a region of constant potential vorticity outside the wave region. Then a third region of potential vorticity, using the known exact zonal flow solutions on the sphere (presented in Appendix A), is chosen so that the velocity remains continuous and the vorticity tends towards zero as the equator is approached.
In all scenarios, we mirror the potential vorticity distribution in the northern hemisphere with that in the southern hemisphere to automatically satisfy the Gauss constraint. We choose 
$\omega$ so that the potential vorticity exactly vanishes at the equator. This allows us to avoid introducing a vortex sheet at the equator. The initial conditions are chosen to illustrate the range of behaviours that these waves exhibit, but we do not perform a thorough examination of parameter space. We start by illustrating a solution where 
$n<0$ that has stable vorticity contours for the relatively small value of 
$|n\epsilon |$ considered. Then we consider a solution where 
$n>0$ but 
$|n\epsilon |$ is taken to be large, so that the waves rapidly overturn and break. Finally, we examine vorticity contours that again have permanent progressive form but in this instance take 
$n>0$.
First, in figure 5, we show the evolution of waves with 
$\epsilon = 0.075$, 
$A_0=0.2$, 
$n=-4$ and 
$\gamma = 10$. Figure 5(a) shows the evolution of the contour at 
$b=0$, in black, and the theoretical prediction in red. The two curves are in agreement. Note that these curves are shown in the stereographic plane for the fully spherical geometry. The evolution of the potential vorticity in the southern hemisphere is shown in figure 5(b).
Figure 5. Contours of constant potential vorticity for a Ptolemaic wave with 
$\epsilon =0.075$, 
$a_0=0.2$, 
$n=-4$ and 
$\gamma = 10$ integrated over a period of revolution. The region in the polar cap has constant potential vorticity. (a) The flow of one contour in black, versus the theoretical prediction given by the dashed red line. (b) The evolution in the southern hemisphere. Note that 
$(\xi,\eta )$ are the stereographic coordinates as defined in Appendix A.
In figure 6, we show the evolution of a wave with 
$\epsilon =0.4$, 
$A_0=0.2$, 
$n=2$ and 
$\gamma = 10$. The evolution of the potential vorticity is shown. These waves are unstable and break almost immediately. The potential vorticity contours become multi-valued, and the flow near the polar vortex becomes turbulent and potential vorticity from higher latitudes is advected to lower latitudes.
Figure 6. Contours of potential vorticity for Ptolemaic waves with 
$\epsilon =0.4$, 
$A_0=0.2$, 
$n=2$ and 
$\gamma = 10$. The region outside the waves starts with uniform potential vorticity, while the flow in the far field is zonal. These waves rapidly go unstable and overturn and break, generating turbulent flow and advecting high potential vorticity into the region outside the waves, showing that mass is being transported by these breaking events.
Finally, in figure 7, we consider a wave with 
$\epsilon =0.03$, 
$A_0=0.2$, 
$n=$, 
$\gamma = 10$, 
$\omega = 0.004$, following the configuration shown in figure 3. The geometry of the solutions is reminiscent of the polar jet stream on the north pole of Saturn, as is discussed further in § 5. The flow retains its structure over the integration time, which exceeds a wave period. Viewed in an inertial frame, the waves are rotating.
Figure 7. Contours of potential vorticity for Ptolemaic waves with 
$\epsilon =0.03$, 
$A_0=0.2$, 
$n=5$ and 
$\gamma = 10$. The flow outside the waves has uniform potential vorticity, before becoming zonal in the far field. The Lagrangian mean flow corresponds to the distribution drawn in figure 3, with a jet-like structure. The configuration is approximately stable over the integration time. The geometry of the solutions is reminiscent of the polar jet stream on the north pole of Saturn, as is discussed further in § 5.
5. Discussion
The generation of the vortex waves considered in this paper, and their degree of symmetry, arises as an accident of the initial conditions. Their geometry, stability and potential vorticity distributions in the numerical simulations invite comparison between these flows and polar jet streams on Earth and other planets, such as Saturn. The meandering of the jet stream is still an active area of research (Nakamura & Huang Reference Nakamura and Huang2018), particularly under a changing climate (Hoskins & Woollings Reference Hoskins and Woollings2015). The flows studied here are two-dimensional and unforced, so any comparison with the real atmosphere is admittedly qualitative.
The Voyager spacecraft first observed the hexagonal flow pattern on Saturn's north pole (Godfrey Reference Godfrey1988). Later, the Cassini mission (Ingersoll Reference Ingersoll2020) captured the hexagon in more detail (Sayanagi et al. Reference Sayanagi, Baines, Dyudina, Fletcher, Sánchez-Lavega and West2018). These observations show that the hexagon pattern is quasi-stationary in a frame rotating with the atmosphere, even though the jet velocity is approximately 125 m s
$^{-1}$. The pattern has been observed for more than three decades, implying a long-time stability. Additionally, the relative vorticity gradient is extremely strong across the jet, i.e. it is a relative vorticity front (see the criteria discussed in § 12.7.3 of Sayanagi et al. Reference Sayanagi, Baines, Dyudina, Fletcher, Sánchez-Lavega and West2018). The hexagonal solutions presented in figure 7 are qualitatively consistent with these criteria, but a quantitative comparison is left for future work.
Exact solutions have been proposed on the sphere in an Eulerian frame (Crowdy & Cloke Reference Crowdy and Cloke2003; Crowdy Reference Crowdy2004) that require singular vorticity distributions in the form of point vortices. The 
$\gamma$ approximation was employed by Siegelman, Young & Ingersoll (Reference Siegelman, Young and Ingersoll2022) to explain the vortex crystal structure found on the pole of Jupiter. The waves in this paper are chosen over parameter regimes so that they avoid singular behaviour in the form of potential vorticity contour cusps or multi-valued solutions.
We have presented two new waves in Euler's equations on the 
$\beta$-plane at the equator and in the 
$\gamma$ approximation. They are motivated by analogous motions on the plane, but differ crucially in that they are describing wave motion where the restoring force is planetary vorticity. Solutions are constructed with prescribed potential vorticity distributions outside the wave fields. Several examples are considered, but our analysis was not a thorough examination of the broad range of solutions that can be constructed. Additionally, only two-dimensional flows have been considered here – analogous exact planar solutions can be unstable to three-dimensional high-frequency perturbations (Leblanc Reference Leblanc2004; Guimbard & Leblanc Reference Guimbard and Leblanc2006).
Acknowledgements
We thank the anonymous referees for their helpful comments, which have led to a significantly improved paper. N.P. thanks Professors J. Wurtele and R. Littlejohn for hosting him in the Department of Physics at UC Berkeley, where part of this work was completed.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of governing equations
We now derive the equations presented in § 2.
A.1. General equations on the sphere in stereographic coordinates
Consider two-dimensional incompressible flow on the unit sphere, where
\begin{equation} X^2+Y^2+Z^2 = 1 \end{equation}
for 
$(X,Y,Z)$ the Cartesian coordinates. We use the stereographic coordinates 
$(\xi,\eta )$ to describe these motions. To derive the governing equations on a surface 
$(X,Y,Z)=(X(\xi,\eta ),Y(\xi,\eta ),Z(\xi,\eta ))$, begin by taking the square of a line element
\begin{equation} {\rm d} s^2 ={\rm d} X_i \,{\rm d} X_i =\frac{\partial X_i}{\partial \xi_j}\, \frac{\partial X_i}{\partial \xi_k}\,{\rm d} \xi_j \,{\rm d} \xi_k=g_{jk}\,{\rm d} \xi_j \,{\rm d} \xi_k, \end{equation}
where 
$(\xi _1,\xi _2)=(\xi,\eta )$, 
$g_{jk}$ is the metric tensor, and Einstein summation is assumed throughout. When the coordinates 
$(\xi,\eta )$ are orthogonal, the metric tensor is diagonal and we have
\begin{equation} {\rm d} s^2 = h_1^2 \,{\rm d} \xi^2 + h_2^2 \,{\rm d} \eta^2, \end{equation}where
\begin{equation} g_{11}\equiv h_1^2 = \left(\frac{\partial X}{\partial \xi}\right)^2+ \left(\frac{\partial Y}{\partial \xi}\right)^2+\left(\frac{\partial Z}{\partial \xi}\right)^2 \end{equation}and
\begin{equation} g_{22}\equiv h_2^2 = \left(\frac{\partial X}{\partial \eta}\right)^2+ \left(\frac{\partial Y}{\partial \eta}\right)^2+\left(\frac{\partial Z}{\partial \eta}\right)^2 \end{equation}
are the tensor components. The kinetic energy of a particle on this surface can be found by dividing (A3) by 
$2({\rm d} t)^2$, so that
\begin{equation} \frac{1}{2}\left(\frac{{\rm d} s}{{\rm d} t}\right)^2 \equiv \frac{1}{2}\,(U^2+V^2) = \frac{1}{2}\left(h_1^2\left(\frac{{\rm d} \xi}{{\rm d} t}\right)^2+h_2^2 \left(\frac{{\rm d} \eta}{{\rm d} t} \right )^2 \right), \end{equation}where the velocities are given by
\begin{equation} U = h_1\,\frac{{\rm d} \xi}{{\rm d} t},\quad V= h_2\,\frac{{\rm d} \eta}{{\rm d} t}. \end{equation}The stereographic coordinates are (see e.g. Needham Reference Needham1997)
\begin{equation} \xi =\frac{X}{1-Z},\quad \eta = \frac{Y}{1-Z}, \end{equation}with inverse mapping
\begin{equation} X=\frac{2\xi}{1+\xi^2+\eta^2},\quad Y=\frac{2\eta}{1+\xi^2+\eta^2},\quad Z = \frac{-1+\xi^2+\eta^2}{1+\xi^2+\eta^2}, \end{equation}and the metric components of this conformal map are
\begin{equation} h_1=h_2 \equiv h = \frac{2}{1+\xi^2+\eta^2}. \end{equation}Next, a unit of area of the fluid is
\begin{equation} {\rm d} A = h^2 \,{\rm d} \xi \,{\rm d} \eta, \end{equation}
so the kinetic energy 
$T$ is
\begin{align} T &= \int \frac{1}{2}\,(U^2+V^2) h^2 \,{\rm d} \xi \,{\rm d} \eta\nonumber\\ &= \int \frac{1}{2}\,(U^2+V^2) J \,{\rm d} a \,{\rm d} b \end{align}
in the Lagrangian reference frame, where 
$(a,b)$ are particle labels, and we define
\begin{equation} J \equiv h^2\,\frac{\partial (\xi,\eta)}{\partial (a,b)}. \end{equation} The time independence of 
$J$ can be established directly by assuming that the area of the fluid does not change. As the fluid has constant density, this implies conservation of mass. Equation (A11) implies
\begin{equation} h'^2 \,{\rm d} \xi' \,{\rm d} \eta' = h^2\,{\rm d} \xi \,{\rm d} \eta, \end{equation}where
\begin{equation} \xi' = \xi(a,b) +\delta \tau\,\xi_{\tau}(a,b),\quad \eta' = \eta(a,b) +\delta \tau\,\eta_{\tau}(a,b) \end{equation}
and 
$h' =h(\xi ',\eta ')$. Equation (A14) is equivalent to
\begin{equation} h'^2\,\frac{\partial (\xi', \eta')}{\partial (a,b)} = h^2\,\frac{\partial (\xi, \eta)}{\partial (a,b)}, \end{equation}
so to first order in 
$\delta \tau$, (A16) implies
\begin{equation} \frac{\partial (h^2\xi_{\tau},\eta )}{\partial (a,b)}- \frac{\partial (h^2\eta_{\tau},\xi )}{\partial (a,b)} = 0, \end{equation}from which it follows that
\begin{equation} \frac{\partial J}{\partial \tau} = 0, \end{equation}
with 
$J$ defined in (A13). For the mapping to be invertible, 
$J$ must not change sign.
To find the vorticity, begin with the circulation 
$\varGamma$ of the fluid:
\begin{equation} \varGamma = \oint \frac{{\rm d} X_i}{{\rm d} t}\,{\rm d} X_ i, \end{equation}where the contour is moving with the fluid on our surface. Using the definition of the fluid velocity and a Green's identity, we have
\begin{align} \varGamma &= \int hU \,{\rm d} \xi +hV \,{\rm d} \eta\nonumber\\ &=\iint \left(\frac{\partial (hV)}{\partial\xi}- \frac{\partial (hU)}{\partial\eta}\right) {\rm d} \xi \,{\rm d} \eta\nonumber\\ &=\iint h^2 q \,{\rm d} \xi\,{\rm d} \eta, \end{align}where
\begin{equation} q =\frac{1}{h^2}\left(\frac{\partial}{\partial \xi} (h^2\eta_{t})-\frac{\partial}{\partial \eta}(h^2\xi_{t})\right). \end{equation}In Lagrangian coordinates, this becomes
\begin{equation} q = \frac{1}{h^2}\left(\frac{\partial(h^2\xi_{\tau},\xi)}{\partial(\xi,\eta)}+ \frac{\partial(h^2\eta_{\tau},\eta)}{\partial(\xi,\eta)} \right) = \frac{1}{h^2}\,\frac{\partial(a,b)}{\partial(\xi,\eta)} \left(\frac{\partial(h^2\xi_{\tau},\xi)}{\partial(a,b)}+ \frac{\partial(h^2\eta_{\tau},\eta)}{\partial(a,b)}\right), \end{equation}which can be rewritten as
\begin{equation} qJ = \frac{\partial(h^2\xi_{\tau},\xi)}{\partial(a,b)}+ \frac{\partial(h^2\eta_{\tau},\eta)}{\partial(a,b)}, \end{equation}
where 
$q=q(a,b)$, which follows from Kelvin's circulation theorem and implies that the vorticity is conserved along fluid particles. When 
$h=1$, the mass and vorticity constraints reduce to those found on the plane.
We rewrite the system in a reference frame rotating at angular velocity 
$\varOmega$ about the 
$z$-axis. This rotation does not change the mass conservation equation, but the potential vorticity becomes (see (34) in Salmon & Pizzo Reference Salmon and Pizzo2023)
\begin{equation} q\to q+2\varOmega z = q+ 2\varOmega\,\frac{-1+\xi^2+\eta^2}{1+\xi^2+\eta^2}, \end{equation}
where in Lagrangian coordinates, 
$q$ is defined by (A23).
Zonal flow may also be written down for this system, and takes the form
\begin{equation} \xi = A_0\sin (ka-(\omega+\varOmega(b))\tau )\,{\rm e}^{kb},\quad \eta = A_0\cos (ka-(\omega+\varOmega(b))\tau )\,{\rm e}^{kb}, \end{equation}where
$$\begin{gather} J= k^2 \,\text{sech}^2 (kb), \end{gather}$$
$$\begin{gather}qJ=2\varsigma \,\text{tanh} (kb) +k \,\text{sech}^2 (kb)\,({-}2k \,\text{tanh} (kb)\, (\omega+\varOmega)+\varOmega'). \end{gather}$$
Just as in the 
$\beta$ and 
$\gamma$ approximations, we may choose 
$\omega$ and 
$\varOmega (b)$ such that the velocity has a prescribed value at the boundary, and the vorticity contours take the form of zonal flow. We have utterly failed at finding more interesting exact solutions on the sphere in Lagrangian coordinates.
A.2. Approximate governing equations
We may perform the asymptotics at the order of the Lagrangian, so that the conservation laws are readily available, or to the final system of equations presented in the previous subsection. We choose the latter for clarity of presentation. There are two assumptions used to simplify the system (see the related discussion in Phillips (Reference Phillips1973), and also the informative overview in Dellar (Reference Dellar2011)). First, we assume that disturbances are small deviations from some rest position. That is,
\begin{equation} \xi = \xi_0+\epsilon\,\tilde{\xi} (a,b,t),\quad \eta = \eta_0+\epsilon\,\tilde{\eta} (a,b,t) \end{equation}
for 
$\epsilon \ll 1$ and 
$(\xi _0,\eta _0)$ constants.
Following Phillips (Reference Phillips1973), our second condition is that 
$\varOmega \gg O(U/2R)$, where 
$U$ is a characteristic velocity scale, and 
$R$ is the planetary radius, here taken to be 1. We impose this constraint on 
$\varOmega$ by writing
\begin{equation} \varOmega = \frac{\omega}{\epsilon}, \end{equation}which yields the classical form of these approximations without additional curvature terms (Dellar Reference Dellar2011).
A.3. The $\beta$-plane
Take 
$\xi _0$ and 
$\eta _0$ to be non-zero so that (A13) and (A23), with 
$q$ defined by (A24), become
\begin{equation} J =\epsilon^2 h_0^2[\tilde{\xi},\tilde{\eta}] +O(\epsilon^3) \end{equation}and
$$\begin{gather} qJ = \epsilon^2 h_0^2 ([\tilde{\xi}_{\tau},\tilde{\xi}]+[\tilde{\eta}_{\tau},\tilde{\eta}])\nonumber\\ \quad {}+8\epsilon h_0^4[\tilde{\xi},\tilde{\eta}] \omega t(1-(\eta_0^2+\xi_0^2)^2+4\epsilon ({-}2+\xi_0^2+\eta_0^2) (\eta_0\widetilde{\eta_1}+\xi_0 \tilde{\xi}))+O(\epsilon^3), \end{gather}$$where
\begin{equation} h_0 = \frac{2}{1+\xi_0^2+\eta_0^2} \end{equation}and
\begin{equation} [\xi,\eta]\equiv \frac{\partial(\xi,\eta)}{\partial(a,b)}. \end{equation}
To illustrate the structure of the equations, first let 
$\xi _0=0$ and 
$\eta _0=1$, corresponding to flow near the equator, which yields the equations
\begin{equation} J=\epsilon^2[\tilde{\xi},\tilde{\eta}] \end{equation}and
\begin{equation} qJ = \epsilon^2 ([\tilde{\xi}_{\tau},\tilde{\xi}]+[\tilde{\eta}_{\tau},\tilde{\eta}])+2 \omega \tilde{\eta} J. \end{equation}
Let 
$\epsilon^2 \hat{J}=J$, and drop the hats for clarity of presentation. Defining (where we recall that the radius of the sphere is 1, which must be taken into account to establish the equivalence of the units of this relationship)
\begin{equation} 2\omega \equiv \beta, \end{equation}the governing equations become
\begin{equation} J = [\xi,\eta] \end{equation}and
\begin{equation} qJ = [\xi_{\tau},\xi]+[\eta_{\tau},\eta] +\beta J \eta, \end{equation}
which describe our 
$\beta -$plane approximation near the equator. When 
$\beta =0$, we return the equations governing motion in the plane (Abrashkin & Yakubovich Reference Abrashkin and Yakubovich1984; Salmon Reference Salmon2020).
Next, let 
$\xi _0=0$ and keep 
$\eta _0$ arbitrary. The equations of motion become
$$\begin{gather} J=h_0^2 \epsilon^2 [\tilde{\xi},\tilde{\eta}], \end{gather}$$
$$\begin{gather}qJ = \epsilon^2h_0^2([\tilde{\xi}_{\tau},\tilde{\xi}]+ [\tilde{\eta}_{\tau},\tilde{\eta}])+2\varOmega J (\mu +\epsilon \lambda \tilde{\eta}), \end{gather}$$where
\begin{equation} \mu = h_0^2(1-h_0),\quad \lambda ={-}4h_0^2\eta_0({-}2+\eta_0^2), \end{equation}and
\begin{equation} h_0=\frac{2}{1+\eta_0^2}. \end{equation}
Defining 
$\epsilon ^2 \hat{J} =J$ and again letting 
$2\omega = \beta$, the governing equations become (dropping the hats again)
$$\begin{gather} J = [\xi,\eta], \end{gather}$$
$$\begin{gather}qJ = [\xi_{\tau},\xi]+[\eta_{\tau},\eta] +\mu\beta J + \beta\lambda \eta J. \end{gather}$$
Redefining 
$q$ as 
$q-\mu \beta$ removes the first term from the right-hand side of the vorticity equation, while 
$\lambda$ can be absorbed into the definition of 
$\beta$, leaving a system of equations that is mathematically equivalent to (A37) and (A38).
A.4. The $\gamma$ approximation
To find our approximate equations of motion near the (south) pole, let 
$\xi _0=\eta _0=0$ and 
$\varOmega =\omega /\epsilon ^2$, so that
$$\begin{gather} J = 4\epsilon^2 [\tilde{\xi},\tilde{\eta}]+O(\epsilon^3), \end{gather}$$
$$\begin{gather}qJ = 8\omega [\tilde{\xi},\tilde{\eta}]+ 4\epsilon^2 ([\tilde{\xi}_{\tau},\tilde{\xi}]+[\tilde{\eta}_{\tau},\tilde{\eta}]) -8 \omega J (\tilde{\xi}^2+\tilde{\eta}^2)+O(\epsilon^3). \end{gather}$$
Let 
$\epsilon ^2 \hat {J}= J$, 
$\hat {\xi }=2\tilde {\xi }$, 
$\hat {\eta }=2\tilde {\eta }$, and define
\begin{equation} 2\omega = \gamma \end{equation}so that our governing equations become, dropping the hats,
$$\begin{gather} J = [\xi,\eta], \end{gather}$$
$$\begin{gather}qJ = f_0 J+ [\xi_{\tau},\xi]+[\eta_{\tau},\eta] -\gamma J (\xi^2+\eta^2), \end{gather}$$
as we concluded from mapping the Eulerian formulation presented in § 2. Note that in the text, we take 
$f_0=\gamma /\epsilon ^2$, which maybe re-absorbed into the definition of 
$q$. This is in agreement with the Eulerian form presented by Nof (Reference Nof1990) and Siegelman et al. (Reference Siegelman, Young and Ingersoll2022).
Open access
