# Stratospheric planetary flows
from the perspective of the Euler equation on a rotating sphere

A. Constantin and P. Germain  Faculty of Mathematics, University of Vienna,
Oskar-Morgenstern-Platz 1, 1090 Vienna, Austria
<EMAIL_ADDRESS>Courant Institute of Mathematical Sciences, New
York University, 251 Mercer Street, New York, NY 10012, USA
<EMAIL_ADDRESS>
###### Abstract.

This article is devoted to stationary solutions of Euler’s equation on a
rotating sphere, and to their relevance to the dynamics of stratospheric flows
in the atmosphere of the outer planets of our solar system and in polar
regions of the Earth. For the Euler equation, under appropriate conditions,
rigidity results are established, ensuring that the solutions are either zonal
or rotated zonal solutions. A natural analogue of Arnold’s stability criterion
is proved. In both cases, the lowest mode Rossby-Haurwitz stationary solutions
(more precisely, those whose stream functions belong to the sum of the first
two eigenspaces of the Laplace-Beltrami operator) appear as limiting cases. We
study the stability properties of these critical stationary solutions. Results
on the local and global bifurcation of non-zonal stationary solutions from
classical Rossby-Haurwitz waves are also obtained. Finally, we show that
stationary solutions of the Euler equation on a rotating sphere are building
blocks for travelling-wave solutions of the 3D system that describes the
leading order dynamics of stratospheric planetary flows, capturing the
characteristic decrease of density and increase of temperature with height in
this region of the atmosphere.

Keywords: inviscid flow, rotating spherical coordinates.

AMS Subject Classifications (2020): 86A10, 76B47, 35Q35.

##  1\. Introduction

###  1.1. Euler’s equation on a rotating sphere

The Euler equation set on the 2-sphere $\mathbb{S}^{2}$, with standard metric,
in a frame rotating at speed $\omega\in\mathbb{R}$ about the polar axis, can
be written in terms of the stream function $\psi$ as

(${\mathcal{E}}_{\omega}$)
$\partial_{t}\,\Delta\psi+\frac{1}{\cos\theta}\,\left[-\partial_{\theta}\psi\,\partial_{\varphi}+\partial_{\varphi}\psi\,\partial_{\theta}\right]\left(\Delta\psi+2\omega\sin\theta\right)=0\,.$

(see for instance Section 13.4.1 in [31]). Here
$(\theta,\varphi)\in(-\frac{\pi}{2},\frac{\pi}{2})\times[0,2\pi)$ are the
latitude and longitude angles (see Figure 1), and $\Delta$ is the Laplace-
Beltrami operator, see Section 2 for a more thorough presentation of the
differential geometry of the sphere.

Figure 1. The rotating spherical coordinate system
$(r^{\prime},\varphi,\theta)$: $\theta\in[-\frac{\pi}{2},\,\frac{\pi}{2}]$ is
the angle of longitude, $\varphi\in[0,\pi]$ is the angle of latitude, and
$r^{\prime}=|OP|$ is the distance from the origin at planet’s center. The
North Pole is at $\theta=\frac{\pi}{2}$, the Equator is on $\theta=0$ and the
South Pole is at $\theta=-\frac{\pi}{2}$.

Conservation laws and energy estimates for equation (${\mathcal{E}}_{\omega}$)
are very similar to the more classical framework of the (two-dimensional)
Euclidean space, or of the torus. Therefore, one can use energy methods to
prove local well-posedness in $H^{s}$, $s>2$ (see for instance [40]), and then
an analogue of the Beale-Kato-Majda theorem [5] ensures global well-posedness
in $H^{s}$, $s>2$. For these matters we refer to Taylor [52], where this
program is explained in greater detail. Global existence being settled, we now
turn to more qualitative questions.

A recent line of research has focused on the behavior of $(E_{\omega})$ as
$\omega\to\infty$, proving convergence to zonal flows after time averaging,
see Cheng-Mahalov [14], Wirosoetisno [56], Taylor [52].

In the present paper, the focus will be on stationary solution and their
stability. Stationary solutions of (${\mathcal{E}}_{\omega}$) can be found by
solving the semilinear elliptic problem

(1.1) $\Delta\psi=F(\psi)-2\omega\sin\theta\,,$

where $F$ is a smooth function; throughout regions without critical points of
the stream function, any stationary solution comes about in this way (see
[16]). Two fundamental explicit classes of stationary solutions of (1.1) are
zonal flows and Rossby-Haurwitz planetary waves:

  * •

Zonal flows correspond to stream functions which only depend on the polar
angle $\theta$, $\psi=\psi(\theta)$. Their stability has been investigated in
Caprino-Marchioro [10] and Taylor [52]. Notice that one can use the invariance
of the equation through the action of $\mathbb{O}(3)$ to obtain non-zonal
solutions from zonal flows. In particular, this explains the apparently ad hoc
Ansatz made in [1].

  * •

Rossby-Haurwitz solutions of degree $k$ are given by the stream function

$\psi_{0}=\alpha\sin\theta+Y(\varphi,\theta),$

where $Y$ belongs to the $k$-th eigenspace $\mathbb{E}_{k}$ of the Laplace-
Beltrami operator and $\alpha\in\mathbb{R}$, solving (1.1) for

$F(\psi)=-k(k+1)\psi\quad\text{and}\quad\omega=\alpha\left(1-\frac{k(k+1)}{2}\right)\,.$

These are the classical Rossby-Haurwitz planetary waves, due to Craig [21],
who found the complete nonlinear solution corresponding to the solutions of
the linearized barotropic vorticity equation obtained by Rossby [47] and
Haurwitz [29] on the beta-plane and on the sphere, respectively (see the
discussion in [50]).

The degree $1$ modes are either zonal (of the form $\beta\sin\theta$) or
rotations of this zonal solution. These solutions can be thought-of as ground
states and will play a key role in this article, since they are distinguished
in many respects:

    * –

They are minimizers of the Dirichlet energy for fixed $L^{2}$ norm. In more
hydrodynamical terms, they minimize the enstrophy
$\int_{{\mathbb{S}}^{2}}|\Delta\psi|^{2}\,{\rm d}\sigma$ for fixed kinetic
energy $\int_{{\mathbb{S}}^{2}}|U|^{2}\,{\rm d}\sigma$, $U$ being the
associated velocity field and ${\rm d}\sigma=\cos\theta\,{\rm d}\theta{\rm
d}\varphi$ being the surface element on the sphere. This gives a first proof
of its stability.

    * –

The $L^{2}$-projection of any solution $\psi$ onto the subspace of ground
states is conserved by the flow of (${\mathcal{E}}_{\omega}$). This is
assertion $(iii)$ in Proposition 1, giving a second proof of the stability of
this ground state.

    * –

This solution is isochronal, in other words its Lagrangian flow is periodic.

The next-gravest modes have degree 2, and they comprise waves with a more
intricate latitude variation. These solutions are of considerable interest in
meteorology. For example, the wave obtained by setting $Y$ proportional to the
spherical harmonics $Y_{2}^{1}$ is commonly observed in the terrestrial
atmosphere, being known as the $5$-day wave since it travels westwards with a
period of about 5 days (see the field data in [30]). These waves are also
preponderant in the atmospheres of the outer planets of our solar system
(Jupiter, Saturn, Uranus, Neptune); see [22]. The instability of the Rossby-
Haurwitz waves is a key factor in the lack of predictability of the weather in
long-term forecasts [8]. Earlier attempts to study the stability of the degree
2 waves by numerical means are somewhat inconclusive, the results being partly
contradicting (see the discussion in [8]). Since the linear stability of the
zonal mode-two solutions was proved by Taylor [52], the issue of the nonlinear
stability or instability of the mode-two Rossby-Haurwitz waves, settled by
Theorem 7, is of outmost importance.

###  1.2. Stratospheric planetary flows

The dynamics of the planetary atmospheres in our solar system is intertwined
with its thermal structure, with temperature gradients driving specific
atmospheric flows which, in turn, modify the temperature field. To a large
extent, Earth’s weather is conditioned by the redistribution of the excess
solar insolation received by the tropical regions towards the poles, whereas
sunlight is not the main driver of the atmospheric motions of Jupiter, Saturn
and Neptune, all these planets radiating about twice as much energy as they
receive from the Sun (see [22, 39]). The atmosphere of these three planets is
primarily made up of a hydrogen-helium gas mixture and the dynamics is
dominated by zonal flows that feature a banded structure – flows of this type
are also common in terrestrial polar regions but the Earth’s atmospheric
circulation at midlatitudes is much more complicated (see [17]). Latitudinal
bands are also the main atmospheric features on Uranus but they are
considerably fewer and only visible in the infrared (see [39]): the lack of an
internal energy source results in less drastic changes of the atmospheric flow
pattern with latitude and an overall rather bland atmosphere – an additional
factor being that Uranus lies sideways, with its poles where its equator
should be, so that its icy core makes the temperatures globally uniformly low,
impeding the formation of localized flow patterns. Note that the strongest
atmospheric flows on any planet were measured on Neptune (reaching supersonic
speeds in equatorial regions). The wind patterns on Neptune and Uranus lack
Jupiter’s multiple zonal winds that flow alternatively in opposite directions:
there is a westward atmospheric flow at low latitudes and an eastward flow at
higher latitudes in each hemisphere (see Figure 2).

Figure 2. Variation of the mean zonal winds with latitude on the giant planets
of our solar system, measured relative to the planet’s rotation speed about
its polar axis (Credit: OpenStax CNX). The traces of methane (which absorbs
red light) in their upper atmosphere gives Uranus and Neptune a blue hue,
obscuring the visibility of specific flow patterns. These pictures show the
high altitude clouds just beneath the stratosphere (at the top of the
troposphere) – the only planetary atmosphere in our solar system transparent
enough to see through from space being that of the Earth.

The stratosphere of the Earth and of the outer planets of our solar system is
thermally stably stratified and presents a rapid decrease of density with
height. In large-scale atmospheric flow conditions, gas parcels move
adiabatically (i.e. without loosing or gaining heat), thus conserving
potential temperature (see [31]). Since stratospheric isentropic surfaces
(level sets of potential temperature) are practically of constant height (see
the discussions in [20, 36]), taking a spherical model for each specific
planet, we see that the motion of inviscid fluids on the surface of a rotating
sphere is relevant to the dynamics of the stratosphere. Consequently, the
study of the Euler equation on a rotating sphere offers insight into the
dynamics of the stratosphere of the outer planets in our solar system. We will
pursue this aspect in some detail in Section 6. Note that this issue is not
relevant for the inner planets of our solar system: Mercury has no atmosphere,
while the atmospheres of Venus and Mars lack a stratosphere. Without a stable
stratification (which is the hallmark of the stratosphere, where the energy
balance is primarily determined by absorption and emission of radiation), two-
dimensional flows on a rotating sphere fail to be pertinent for atmospheric
dynamics.

It is well-established that the Rossby-Haurwitz waves play an important role
in the large-scale dynamics of atmospheric flows. In particular, the zonal
spherical harmonics capture the pattern of the band structure in the upper
atmosphere of the outer planets, comprised of zonal flows whose direction
alternates from westward to eastward. But superimposed on these flows one
often observes non-zonal features (referred to as “eddies” in the atmospheric
sciences), like Jupiter’s Great Red Spot. First thought of as an exotic
occurrence, it is now appreciated that such long-lived vortices are frequently
encountered throughout the solar system. For example, several spots occur in
Saturn’s and Neptune’s atmosphere (like the Polar Hexagon and the Great Dark
Spot, respectively), while in Jupiter’s atmosphere there are upwards of a
dozen smaller vortices near the latitude 22∘S (where the Great Red Spot is
centred); see the discussion in [22]. The spherical harmonics by themselves
cannot cover this plethora of flows. It is therefore of interest to develop an
approach that can provide large families of Rossby-Haurwitz-like solutions.
This is precisely our aim when pursuing bifurcation from Rossby-Haurwitz
waves. Let us also point out that the scarcity of stable atmospheric flows
makes flows likely to be unstable also of great interest. In this context,
note that Jupiter’s Great Red Spot is confined by an eastward jet stream to
its south and a westward one to its north (which explains why it rotates
counterclockwise), and while a change of direction in zonal flows if often
indicative of their linear instability by means of a variant of Rayleigh’s
criterion (see [52]), only small changes were noted in the dynamics of the
Great Red Spot since 1831. This shows that even potentially unstable flow
patterns can be (relatively) long-lived. We would also like to emphasise the
importance of not performing the analysis within the flat geometry of the
$f$-plane or $\beta$-plane approximation. Even the Great Red Spot (large
enough to engulf Earth) is not an isolated vortex but rather a global system
involving involving one large anticyclone, and several smaller ones in the
same anticyclone zone, and a filamentary region in the adjacent equatorward
shear zone (see [22]).

###  1.3. Main results and organization of the article

####  1.3.1. Symmetry

Our first result addresses symmetries of the solutions of (1.1).

######  Theorem 1.

Consider $\psi$ a solution of (1.1) for some $\omega\in\mathbb{R}$. If
$F^{\prime}>-6$, then $\psi$ is a zonal flow, modulo a rotation in
$\mathbb{O}(3)$.

This theorem is an immediate consequence of Theorem 4 in Section 3. Note that
the number $-6$ is the second eigenvalue of the Laplace-Beltrami operator.
Such symmetry results are proved in Constantin-Drivas-Ginsberg [15] for
general Riemannian surfaces with Killing field, under the condition that
$F^{\prime}$ is larger than the smallest eigenvalue of the Laplace-Beltrami
operator, which on the sphere amounts to $F^{\prime}>-2$. The symmetric
structure of the sphere explains the improvement that we are able to obtain.

The sharpness of this result can be seen from the Rossby-Haurwitz solutions in
$\mathbb{E}_{1}+\mathbb{E}_{2}$ (the first and second eigenspaces of
$-\Delta$). Modulo $\mathbb{O}(3)$, they are of the type
$\psi_{0}=\alpha\sin\theta+Y$, where $\alpha\in\mathbb{R}$ and
$Y\in\mathbb{E}_{2}$, and correspond to $F^{\prime}=-6$; they are not zonal in
general. This is a first indication that Rossby-Haurwitz solutions in
$\mathbb{E}_{1}+\mathbb{E}_{2}$ play a key role in the qualitative study of
the 2D Euler problem on the sphere.

####  1.3.2. Stability of zonal solutions

Persistent zonal (east-west) flows are ubiquituous in planetary atmospheres.
While the Earth’s atmosphere typically exhibits meandering jets, the zonal
flows on the outer planets present a long-lived coherence that might be
indicative of stability. In Section 4, after a brief discussion of the
Rayleigh and Fjortoft necessary criteria for the linear instability of zonal
flows on a sphere, we present the proof of an Arnold-type nonlinear stability
result for zonal flows. While variations of this result can be found in the
research literature, with the possibly most comprehensive approach for general
rotating bodies due to Taylor [52], we present a simple proof for the case of
a rotating sphere. We also show that this result yields the stability of the
mean zonal flow patterns on Uranus and Neptune.

####  1.3.3. Stability of Rossby-Haurwitz solutions of degree 2

The classical theory of Arnold [3, 4] adresses the stability of stationary
solutions of the Euler equation on a planar domain $\Omega$. Let $\lambda$
denote the largest eigenvalue of the negative Dirichlet Laplacian $-\Delta$ on
$\Omega$. Then a solution of $\Delta\psi=F(\psi)$ is called Arnold stable of
type I if $-\lambda<F^{\prime}<0$, and Arnold stable of type II if
$F^{\prime}>0$. These solutions can be showed to be nonlinearly stable for the
2D Euler problem, through the construction of an energy functional.

Due to its rich symmetry structure, the relevant threshold for the sphere
becomes the second eigenvalue of $-\Delta$, instead of the first as in the
classical Arnold theory. Namely, the following result holds.

######  Theorem 2.

If $\psi$ is a solution of (1.1) with $-6<F^{\prime}<0$, then it is stable in
$H^{2}({\mathbb{S}}^{2})$.

This theorem corresponds to Theorem 5 in Section 4. It is natural to
conjecture that the transition to instability occurs as $\min\\{F^{\prime}\\}$
crosses the threshold of $-6$. This points to the importance of Rossby-
Haurwitz solutions in $\mathbb{E}_{1}+\mathbb{E}_{2}$, for which
$F^{\prime}=-6$. It was observed following Theorem 1 that these flows are in
some sense the ”first” genuine non-zonal stationary flows; we now see that
they appear at the transition to instability since we are able to analyze the
stability of Rossby-Haurwitz solutions of degree 2 precisely.

######  Theorem 3.

  * (i)

The set $\mathbb{E}_{1}+\mathbb{E}_{2}$ is stable in
$H^{2}({\mathbb{S}}^{2})$.

  * (ii)

Non-zonal Rossby-Haurwitz solutions of degree 2 are unstable in
$H^{2}({\mathbb{S}}^{2})$

  * (iii)

Zonal Rossby-Haurwitz solutions of degree 2 are stable in
$H^{2}({\mathbb{S}}^{2})$.

This theorem is a short version of Theorem 7 in Section 5. The notion of
stability considered here is always Lyapunov stability.

####  1.3.4. Bifurcation

Exact solutions of the vorticity equation (${\mathcal{E}}_{\omega}$) are very
useful for gaining insight into the dynamics, for validating models and for
evaluating numerical discretizations. Due to the scarcity of the available
explicit exact non-zonal solutions, we develop a bifurcation approach that
permits the construction of families of exact solutions to (1.1) for some
classes of parameter-dependent nonlinearities $F$ with specific structural
properties. To ensure that these solutions comprise non-zonal flows, we
implement a symmetry-breaking approach. We also derive a priori bounds of the
velocity and of the vorticity along the continuum of solutions constructed by
means of bifurcation (see Theorems 8-10).

####  1.3.5. Stratospheric flows

We develop a methodology for embedding the 2D flows studied hitherto into the
3D dynamics of the stratosphere. More precisely, solutions to
$\Delta\psi=F(\psi)$ on ${\mathbb{S}}^{2}$ correspond to geostationary
solutions $\widehat{\psi}=\psi(\varphi+\omega t,\theta)$ of the vorticity
equation (${\mathcal{E}}_{\omega}$) which are restrictions to the surface of
the sphere of solutions to the 3D system that describes the leading-order
dynamics of the stratosphere (see Theorem 11), where the density decreases and
the temperature increases with height above the tropopause. This feature
appears to be replicated by solutions to (1.1) for suitable forcings of the
gravity acting in the radial direction (typically encountered if the
geopotential surfaces are not spheres) but we do not pursue this direction in
the present paper.

###  1.4. Comparison with the torus

It is instructive to compare the results which have been stated above to what
is known for 2D-Euler on the torus, for which we adopt the parametrization
$(x,y)\in[0,2\pi]^{2}$ with periodic boundary conditions, the metric being the
standard one.

From a geometric and analytic viewpoint, the sphere and the torus share many
similarities: the eigenspaces of the Laplacian can be described explicitly and
isometries act transitively in both cases. However, the sphere is more
symmetric than the torus, since, for instance, all geodesics are equivalent up
to isometries; this ultimately results in a different picture for the
stability of stationary solutions of the Euler equation.

On the torus, there are two important classes of explicit stationary
solutions: shear flows, whose stream function only depends on $y$, and
eigenfunctions of the Laplacian; these are analogous to zonal flows and
Rossby-Haurwitz solutions, respectively.

The first eigenvalue of $-\Delta$ is $1$, with a distinguished eigenfunction
provided by $\sin y$, and also a shear flow (sometimes called Kolmogorov
flow). It is analogous to the stationary solution $\sin\theta$ on the sphere.
Just like for the sphere, solutions of $-\Delta\psi=F(\psi)$ with
$F^{\prime}>-1$ are constant. However, there are elements of the first
eigenspace of the Laplacian which are not shear, and there is furthermore a
rich and nontrivial family of stationary solutions bifurcating from $\sin y$,
see Coti-Zelati-Elgindi-Widmayer [18].

It is instructive to consider other aspect ratios: if we change the
parameterization to $[0,2\pi L]\times[0,2\pi]$, the Kolmogorov flow $\sin y$
is still a stationary solution, but new phenomena occur. First, the only
stationary solutions in a small neighborhood are shear flows [18]. Second, the
stability properties of the Kolmogorov flow depend on the aspect ratio $L$:

  * •

the flow is stable for $L>1$ for obvious energy reasons;

  * •

it can also be proved to be stable for $L=1$, see [4];

  * •

linear instability for $L<1$ has been the subject of a number of works,
starting with the sketch in Arnold and Meshalkin [2], followed by Belenkaya-
Friedlander-Yudovich [7] and Buttà-Negrini [9].

Finally, we mention that inviscid damping for the linearized problem around
the Kolmogorov flows was studied in Wei-Zhang-Zhao [54]; such questions are
certainly harder for the sphere, due to the more involved eigenfunction
decomposition for the Laplacian.

### Acknowledgements

The authors thank Vladimir Šverak for a very helpful communication on the
prescribed Gauss curvature equation on the sphere. While working on this
project, PG was visiting the University of Vienna, to which he is very
grateful. PG was supported by the NSF grant DMS-1501019, by the Simons
collaborative grant on weak turbulence, and by the Center for Stability,
Instability and Turbulence (NYUAD). AC was supported by the Austrian Science
Foundation (FWF) ”Wittgenstein-Preis Z 387-N”.

##  2\. The vorticity equation for inviscid flow on a rotating sphere

###  2.1. Differential geometry of the sphere

Since by the “hairy ball theorem” the $2$-sphere $\mathbb{S}^{2}$ does not
possess a continuously differentiable field of unit tangent vectors (see
[28]), it is not possible to cover $\mathbb{S}^{2}$ with one chart. However,
the standard longitude-latitude spherical coordinates
$(\varphi,\theta)\in(0,2\pi)\times(-\frac{\pi}{2},\frac{\pi}{2})$ provide us
(see Fig. 1) with a chart

$(\varphi,\theta)\in(0,2\pi)\times\big{(}-\tfrac{\pi}{2},\tfrac{\pi}{2}\big{)}\mapsto(\cos\varphi\cos\theta,\,\sin\varphi\cos\theta,\,\sin\theta)\in\mathbb{S}^{2}$

covering $\mathbb{S}^{2}$ with the half-circle $\varphi=0$ (the international
date line, including the poles) excised. A smooth atlas for $\mathbb{S}^{2}$
is obtained by coupling this with the chart

$(\tilde{\varphi},\tilde{\theta})\in(0,2\pi)\times\big{(}-\tfrac{\pi}{2},\tfrac{\pi}{2}\big{)}\mapsto(-\cos\tilde{\varphi}\cos\tilde{\theta},\,\sin\tilde{\theta},\,\sin\tilde{\varphi}\cos\tilde{\theta})\in\mathbb{S}^{2}$

covering $\mathbb{S}^{2}$ with the equatorial half-circle parametrized in
spherical coordinates by $\\{\theta=0\,,\,\varphi\in[\pi,2\pi]\\}$ excised:
the bijective transformation $(x,y,z)\mapsto(-x,z,y)$ between the above
parametrizations is equivalent to first rotating the Euclidean coordinate
system by $\frac{\pi}{2}$ about the $x$-axis and then by $\pi$ about the
$z$-axis.

Throughout this paper we will mostly rely only on spherical coordinates. The
double-valued ambiguity along the international date line of the chart
provided by the spherical coordinates can be resolved by assuming a periodic
dependence on the azimuthal angle $\varphi$. At several places in the
manuscript the use of spherical coordinates (more precisely, the fact that
longitude is not well-defined and latitude circles degenerate into a single
point at the poles) introduces artificial singularities at the poles that can
be ruled out either by switching to the chart that covers the sphere with the
equatorial half-circle removed or by taking smoothness into account – see
relation (2.7) below.

The $4$-dimensional tangent bundle $T\mathbb{S}^{2}$ of the $2$-sphere
$\mathbb{S}^{2}$ is not parallelizable as a consequence of the hairy ball
theorem. However, at every point $X$ of $\mathbb{S}^{2}\setminus\\{N,S\\}$,
having spherical coordinates
$(\varphi,\theta)\in[0,2\pi]\times(-\frac{\pi}{2},\frac{\pi}{2})$, the tangent
vectors

$\mathbf{e}_{\varphi}=\frac{1}{\cos\theta}\,\partial_{\varphi},\qquad\mathbf{e}_{\theta}=\partial_{\theta}\,,$

provide us with a basis of the tangent space $T_{X}\mathbb{S}^{2}$ at
$X\in\mathbb{S}^{2}$. In these coordinates, the Riemannian volume element is

${\rm d}\sigma=\cos\theta\,{\rm d}\varphi\,{\rm d}\theta\,,$

and the classical differential operators (gradient and Laplace-Beltrami for
scalar functions $\psi:\mathbb{S}^{2}\to{\mathbb{R}}$, divergence for vector
fields $F:\mathbb{S}^{2}\to T\mathbb{S}^{2}$) are given by

$\displaystyle\operatorname{grad}\psi=\partial_{\theta}\psi\,\mathbf{e}_{\theta}+\frac{\partial_{\varphi}\psi}{\cos\theta}\,\mathbf{e}_{\varphi}\,,$
$\displaystyle\operatorname{div}(F_{\varphi}\,\mathbf{e}_{\varphi}+F_{\theta}\,\mathbf{e}_{\theta})=\frac{1}{\cos\theta}\,[\partial_{\varphi}F_{\varphi}+\partial_{\theta}(\cos\theta\,F_{\theta})]\,,$
$\displaystyle\operatorname{\Delta}\psi=\operatorname{div}\operatorname{grad}\psi=\partial_{\theta}^{2}\psi-\tan\theta\,\partial_{\theta}\psi+\frac{1}{(\cos\theta)^{2}}\,\partial_{\varphi}^{2}\psi\,,$

with the formula for $\operatorname{grad}$ following from the definition
$\partial_{s}\psi(\gamma(s))=\operatorname{grad}\psi\cdot\gamma^{\prime}(s)$
for a path $\gamma$ on $\mathbb{S}^{2}$, while the formula for
$\operatorname{div}$ follows by duality (see [46]). The covariant derivatives
have the form

$\nabla_{\mathbf{e}_{\theta}}\mathbf{e}_{\theta}=\nabla_{\mathbf{e}_{\theta}}\mathbf{e}_{\varphi}=0\,,\qquad\nabla_{\mathbf{e}_{\varphi}}\mathbf{e}_{\theta}=-\tan\theta\,\mathbf{e}_{\varphi}\,,\qquad\nabla_{\mathbf{e}_{\varphi}}\mathbf{e}_{\varphi}=\tan\theta\,\mathbf{e}_{\theta}\,,$

and can be computed by projecting Euclidean derivatives on the tangent space
to the $2$-sphere: for instance
$\nabla_{\mathbf{e}_{\varphi}}\mathbf{e}_{\theta}=P\partial_{\varphi}\mathbf{e}_{\theta}$,
where $P$ is the projection operator. Note also that the $2$ sphere
$\mathbb{S}^{2}$ admits the complex structure $J$ (corresponding to a rotation
in the tangent space) defined by

$J\mathbf{e}_{\varphi}=\mathbf{e}_{\theta},\qquad
J\mathbf{e}_{\theta}=-\mathbf{e}_{\varphi}\,.$

The Laplace-Beltrami operator $\Delta$ on $\mathbb{S}^{2}$, operating in the
Hilbert space $L^{2}({\mathbb{S}}^{2})$ obtained as the completion of the
smooth functions $f:{\mathbb{S}}^{2}\to{\mathbb{C}}$ of zero mean (i.e., with
$\iint_{{\mathbb{S}}^{2}}f\,{\rm d\sigma}=0$) with respect to the inner
product

$\langle
f_{1},f_{2}\rangle=\iint_{{\mathbb{S}}^{2}}f_{1}\,\overline{f_{2}}\,{\rm
d\sigma}\,,$

(where the overbar denotes complex conjugation) is negative, self-adjoint and
its spectrum is the discrete set of eigenvalues $\bigcup_{j\geq
1}\\{-j(j+1)\\}$, the spherical harmonics $\\{Y_{j}^{m}\\}_{j\geq 1,\,|m|\leq
j}$ being an orthonormal basis of eigenfunctions in $L^{2}({\mathbb{S}}^{2})$,
with

$\Delta Y_{j}^{m}=-j(j+1)Y_{j}^{m}\,,\qquad j\geq 1\,,\ m\in\\{-j,\dots,j\\}$

(see the Appendix).

###  2.2. Euler equation

We start from the stream function $\psi$, from which the velocity field $U$ is
obtained as

$U=J\operatorname{grad}\psi\quad\text{or}\quad
U=u\mathbf{e}_{\varphi}+v\mathbf{e}_{\theta}\,,$

with the geostrophic relations

(2.1) $\begin{cases}\displaystyle u=-\partial_{\theta}\psi\,,\\\ \displaystyle
v=\tfrac{1}{\cos\theta}\,\partial_{\varphi}\psi\,.\end{cases}$

The vorticity is then given by

(2.2) $\Omega=\Delta\psi=-\operatorname{div}JU\,,$

while the material derivative, describing the transport by the velocity field
$U$, can be expressed in the form

$D_{t}=\partial_{t}+\nabla_{U}=\partial_{t}+u\nabla_{\mathbf{e}_{\varphi}}+v\nabla_{\mathbf{e}_{\theta}}\,.$

The material derivative can be applied to scalars or vectors; when applied to
scalars, it becomes

$D_{t}=\partial_{t}+\frac{1}{\cos\theta}\left[-\partial_{\theta}\psi\,\partial_{\varphi}+\partial_{\varphi}\psi\,\partial_{\theta}\right]\,.$

The Euler equation on a sphere rotating at speed $\omega$ about the polar axis
can then be written for either of the variables $\psi$, $U$, or $\Omega$. At
the level of the stream function, it is equation (${\mathcal{E}}_{\omega}$):

$\partial_{t}\,\Delta\psi+\frac{1}{\cos\theta}\,\left[-\partial_{\theta}\psi\partial_{\varphi}+\partial_{\varphi}\psi\partial_{\theta}\right]\left(\Delta\psi+2\omega\sin\theta\right)=0\,.$

To express the Euler equation in terms of the velocity field we have to add
the divergence-free condition to the evolution equation
(${\mathcal{E}}_{\omega}$), together with an auxiliary scalar pressure field
$p$ (that arises as a Lagrange multiplier for the divergence-free constraint)

(2.3)
$\left\\{\begin{array}[]{l}D_{t}U+2\omega\sin\theta\,JU=-\operatorname{grad}p\,,\\\
\operatorname{div}U=0\,.\end{array}\right.$

At the level of the vorticity, the evolution equation becomes111Sometimes
termed the two-dimensional baroclinic Ertel equation for the material
conservation of potential vorticity – see [31].

(2.4) $D_{t}(\Omega+2\omega\sin\theta)=0\,,$

and has to be complemented with the Biot-Savart law, which recovers at every
instant $t$ the stream function (and thus the velocity field) from the
vorticity:

$\psi(\xi_{0})=\iint_{\mathbb{S}^{2}}{\mathcal{G}}(\xi,\xi_{0})\Omega(\xi)\,{\rm
d}\sigma(\xi)\,,$

where (see [23, 42])

${\mathcal{G}}(\xi,\xi_{0})=\frac{1}{2\pi}\,\ln\Big{(}\frac{|\xi-\xi_{0}|}{2}\Big{)}\quad\text{with}\quad|\xi-\xi_{0}|\quad\text{the
distance in}\ {\mathbb{R}}^{3}\ \text{between}\ \xi\neq\xi_{0}\ \text{on}\
\mathbb{S}^{2}$

is the Green function, satisfying

(2.5) $\Delta\mathcal{G}(\xi,\xi_{0})=\delta(\xi-\xi_{0})-\frac{1}{4\pi}\,,$

with $\delta$ the Dirac delta distribution corresponding to a point vortex
located at $\xi_{0}\in\mathbb{S}^{2}$. Note that since the velocity field is
divergence-free, an immediate consequence of the divergence theorem is the
validity of the Gauss constraint

(2.6) $\iint_{\mathbb{S}^{2}}\Omega\,{\rm d}\sigma=0\,,$

so that the factor $-\frac{1}{4\pi}$ in (2.5) plays the role of a compensating
uniform vorticity distribution on $\mathbb{S}^{2}$ to guarantee the validity
of (2.6). Regarding the apparent singularity of the meridional velocity
component $v$ in (2.1) at the poles, let us point out that for any
$C^{1}$-function $\psi:{\mathbb{S}}^{2}\to{\mathbb{R}}$ the continuity of the
gradient with respect to the spherical coordinates $(\varphi,\theta)$ implies

(2.7)
$\lim_{\theta\to\pm\frac{\pi}{2}}\partial_{\varphi}\psi(\varphi,\theta)=0$

since on any genuine circle of latitude
$\theta\in\big{(}-\tfrac{\pi}{2},\tfrac{\pi}{2}\big{)}$ the periodicity of
$\psi$ in the longitudinal direction ensures the existence of a point where
$\partial_{\varphi}\psi$ vanishes.

Generally, insight in the flow dynamics is more readily available working with
the stream function $\psi$, rather than with the vorticity
$\Omega=\Delta\psi$. Equation (${\mathcal{E}}_{\omega}$) is the barotropic
vorticity equation, describing the motion of an inviscid, unforced,
incompressible, homogeneous fluid on a rotating sphere (see [25]). Note that
with respect to the symplectic structure on $\mathbb{S}^{2}$, whose Poisson
bracket is given in spherical coordinates by

(2.8)
$\\{f,h\\}=\frac{1}{\cos\theta}\,\big{(}\partial_{\theta}h\,\partial_{\varphi}f-\partial_{\theta}f\,\partial_{\varphi}h\big{)}\,,$

the vorticity equation (${\mathcal{E}}_{\omega}$) can be expressed as the
Hamiltonian flow

$\partial_{t}(\Delta\psi+2\omega\sin\theta)=\\{\Delta\psi+2\omega\sin\theta,\psi\\}\,.$

###  2.3. Symmetries

The Euler equations (${\mathcal{E}}_{\omega}$) with different rotation speeds
$\omega$ are related through the change-of-frame transformation

(2.9)
$\psi_{0}(\varphi,\theta,t)\longleftrightarrow\psi_{\omega}(\varphi,\theta,t)=\psi_{0}(\varphi+\omega
t,\theta,t)+\omega\sin\theta\,.$

More precisely, $\psi_{\omega}$ solves (${\mathcal{E}}_{\omega}$) if and only
if $\psi_{0}$ solves the Euler equation on a fixed sphere,
(${\mathcal{E}}_{0}$).

The classical scaling of the Euler equation in a fixed frame has to be
modified to take $\omega$ into account: namely, if $\psi(\varphi,\theta,t)$
solves (${\mathcal{E}}_{\omega}$), then $\lambda\psi(\varphi,\theta,\lambda
t)$ solves (${\mathcal{E}}_{\lambda\omega}$), for any $\lambda>0$.

Another invariance is related to the symmetries of the 2-sphere, given by the
orthogonal group $\mathbb{O}(3)$, a compact Lie group of dimension 3,
consisting of the isometries of ${\mathbb{R}}^{3}$ which fix the origin: one
can think of $\mathbb{O}(3)$ as the group of orthogonal real $3\times 3$
matrices or as a group of transformations of ${\mathbb{R}}^{3}$. The action of
${\mathbb{O}}(3)$ is defined by

$Gf(X)=f(GX)\,,\qquad X\in{\mathbb{S}}^{2}\,,\quad G\in{\mathbb{O}}(3)\,,$

for a scalar function $f:{\mathbb{S}}^{2}\to{\mathbb{R}}$ and the following
transformations leave the set of solutions of (${\mathcal{E}}_{\omega}$)
invariant:

$\begin{cases}&\psi(X,t)\mapsto\psi(GX,t)\,,\\\ &U(X,t)\mapsto GU(GX,t)\,,\\\
&\Omega(X,t)\mapsto\Omega(GX,t)\,,\end{cases}\qquad X\in{\mathbb{S}}^{2}\,.$

Note that the non-abelian subgroup of $\mathbb{O}(3)$ of all orthogonal
$3\times 3$ real matrices $R$ with $\text{det}(R)=1$, itself a compact Lie
group of dimension 3, is called the rotation group $\mathbb{SO}(3)$ since each
transformation $X\mapsto RX$ with $R\in\mathbb{SO}(3)$ can be obtained by
first choosing a fixed direction through the origin and subsequently rotating
the coordinate system through a suitable angle about this direction as an axis
(see [46] and the Appendix).

###  2.4. Conservation laws

The identification of integrals of motion provides insight into the flow
dynamics.

######  Proposition 1 (Integrals of motion for a fixed sphere).

The following quantities are conserved by smooth solutions of the vorticity
equation (${\mathcal{E}}_{0}$):

  * (i)

(kinetic energy)
$\displaystyle\frac{1}{2}\iint_{{\mathbb{S}}^{2}}|U|^{2}\,{\rm d}\sigma$,

  * (ii)

(Casimir invariants) $\displaystyle\iint_{{\mathbb{S}}^{2}}F(\Omega)\,{\rm
d}\sigma$ for any differentiable function $F$,

  * (iii)

(first eigenspace of the Laplace-Beltrami operator) the vorticity components
in the direction of each of the three spherical harmonics of degree $1$.

###### Proof.

(i) Taking the time derivative of the energy gives, with the help of (2.3),

$\frac{d}{dt}\Big{(}\frac{1}{2}\iint_{{\mathbb{S}}^{2}}|U|^{2}\,{\rm
d}\sigma\Big{)}=\iint_{{\mathbb{S}}^{2}}[-\nabla_{U}U-2\omega\sin\theta
JU-\operatorname{grad}p]\cdot U\,{\rm d}\sigma\,.$

It is immediate to see that the second term on the right side vanishes (due to
the antisymmetry of $J$), and the third term as well (since $U$ is divergence-
free). But the first term on the right side is also zero since

$\iint_{{\mathbb{S}}^{2}}\nabla_{U}U\cdot U\,{\rm
d}\sigma=\iint_{{\mathbb{S}}^{2}}U\cdot\operatorname{grad}\frac{|U|^{2}}{2}\,{\rm
d}\sigma=-\iint_{{\mathbb{S}}^{2}}\operatorname{div}U\,\frac{|U|^{2}}{2}\,{\rm
d}\sigma=0\,.$

(ii) This is an immediate consequence of the vorticity equation (2.4), due to
a simple change of variables in phase-space since the flow-map is area
preserving ($U$ being divergence-free). Note that these Casimir functionals
reflect the underlying noncanonical Hamiltonian structure induced by (2.8):
their Poisson bracket with any other functional vanishes.

(iii) Since the spherical harmonics $Y_{1}^{\pm 1}$ can be obtained from
$Y_{1}^{0}$ by a rotation, due to the invariance of the vorticity equation
under the action of $\mathbb{SO}(3)$, it suffices to show that
$\iint_{{\mathbb{S}}^{2}}\Omega\sin\theta\,{\rm d}\sigma$ is conserved. Since
$\operatorname{grad}\sin\theta=\cos\theta\,\mathbf{e}_{\varphi}$, an
integration by parts reveals that

$\iint_{{\mathbb{S}}^{2}}\Omega\sin\theta\,{\rm
d}\sigma=-\iint_{{\mathbb{S}}^{2}}\operatorname{div}JU\sin\theta\,{\rm
d}\sigma=\iint_{{\mathbb{S}}^{2}}U\cdot\mathbf{e}_{\varphi}\cos\theta\,{\rm
d}\sigma\,.$

Taking the time derivative of this quantity and using the velocity evolution
equation (2.3) gives

$\frac{d}{dt}\iint_{{\mathbb{S}}^{2}}U\cdot\mathbf{e}_{\varphi}\cos\theta\,{\rm
d}\sigma=\iint_{{\mathbb{S}}^{2}}\left[-\nabla_{U}U-2\omega\sin\theta
JU-\operatorname{grad}p\right]\cdot\mathbf{e}_{\varphi}\cos\theta\,{\rm
d}\sigma\,.$

The third term on the above right-hand side is zero, since
$\operatorname{div}(\cos\theta\,\mathbf{e}_{\varphi})=0$. The second term is
also zero, as can be seen by using the defition of $U$ in terms of $\psi$, and
the fact that
$\operatorname{div}(\sin\theta\cos\theta\,\mathbf{e}_{\varphi})=0$:

$\iint_{{\mathbb{S}}^{2}}\sin\theta
JU\cdot\mathbf{e}_{\varphi}\,\cos\theta\,{\rm
d}\sigma=\iint_{{\mathbb{S}}^{2}}\operatorname{grad}\psi\cdot(\sin\theta\cos\theta\,\mathbf{e}_{\varphi})\,{\rm
d}\sigma=0\,.$

Finally, using the identity
$\nabla_{U}(\cos\theta\,\mathbf{e}_{\varphi})=JU\sin\theta$, the fact that $U$
is divergence-free, and the antisymmetry of $J$, one obtains

$\iint_{{\mathbb{S}}^{2}}\nabla_{U}U\cdot\mathbf{e}_{\varphi}\cos\theta\,{\rm
d}\sigma=-\iint_{{\mathbb{S}}^{2}}U\cdot\nabla_{U}(\mathbf{e}_{\varphi}\cos\theta)\,{\rm
d}\sigma=-\iint_{{\mathbb{S}}^{2}}U\cdot JU\sin\theta\,{\rm d}\sigma=0\,,$

so that the third term vanishes as well.∎

Due to the transformation (2.9), Proposition 1 gives invariants for
(${\mathcal{E}}_{\omega}$) with $\omega\neq 0$.

######  Corollary 1 (Integrals of motion for a rotating sphere).

The following quantities are constant in time for smooth solutions of
(${\mathcal{E}}_{\omega}$) with $\omega\neq 0$:

  * (i)

(kinetic energy) $\displaystyle\int|U|^{2}\,dx$,

  * (ii)

(Casimir invariants) $\displaystyle\int F(\Omega+2\omega\sin\theta)\,dx$ for
any differentiable function $F$,

  * (iii)

(first eigenspace of the Laplace-Beltrami operator) ${\rm e}^{{\rm i}m\omega
t}c_{1}^{m}(t)$ for the coefficients

$c_{1}^{m}=\iint_{{\mathbb{S}}^{2}}\Omega\,\overline{Y_{1}^{j}}\,{\rm
d}\sigma\,,\qquad m\in\\{-1,0,1\\}\,,$

of the $L^{2}({\mathbb{S}}^{2})$-expansion of $\Omega$ in terms of the
spherical harmonics $\\{Y_{j}^{m}\\}_{j\geq 1,\,|m|\leq j}$. In particular,
the real number $c_{1}^{0}(t)$ and the absolute values of the complex numbers
$c_{1}^{\pm 1}(t)$ are flow-invariants.

###### Proof.

The proof of (i)-(ii) being rather straightforward, we only discuss that of
(iii). Since $\sin\theta=2\sqrt{\frac{\pi}{3}}\,Y_{1}^{0}(\theta)$, using
(2.9) and the explicit dependence of the spherical harmonics $Y_{1}^{m}$ on
the longitude angle $\varphi$, for a solution
$\psi_{\omega}(\varphi,\theta,t)$ of (${\mathcal{E}}_{\omega}$) we get

$\displaystyle{\rm e}^{{\rm i}m\omega t}c_{1}^{m}(t)$ $\displaystyle={\rm
e}^{{\rm i}m\omega
t}\iint_{{\mathbb{S}}^{2}}\psi_{\omega}(\cdot,\cdot,t)\,\overline{Y_{1}^{m}}\,{\rm
d}\sigma$ $\displaystyle={\rm e}^{{\rm i}m\omega
t}\iint_{{\mathbb{S}}^{2}}\psi_{0}(\cdot+\omega
t,\cdot,t)\,\overline{Y_{1}^{m}}\,{\rm
d}\sigma+2\omega\sqrt{\tfrac{\pi}{3}}\,{\rm e}^{{\rm i}m\omega
t}\iint_{{\mathbb{S}}^{2}}Y_{1}^{0}\,\overline{Y_{1}^{m}}\,{\rm d}\sigma$
$\displaystyle=\iint_{{\mathbb{S}}^{2}}\psi_{0}(\cdot,\cdot,t)\,\overline{Y_{1}^{m}}\,{\rm
d}\sigma+2\omega\sqrt{\tfrac{\pi}{3}}\,{\rm e}^{{\rm i}m\omega
t}\iint_{{\mathbb{S}}^{2}}Y_{1}^{0}\,\overline{Y_{1}^{m}}\,{\rm d}\sigma$

and we can conclude by Proposition 1, the spherical harmonics being
orthonormal.∎

###  2.5. Stationary solutions

Stationary solutions of (${\mathcal{E}}_{\omega}$) satisfy

(2.10)
$\left[-\partial_{\theta}\psi\,\partial_{\varphi}+\partial_{\varphi}\psi\,\partial_{\theta}\right]\left(\Delta\psi+2\omega\sin\theta\right)=0$

In view of the transformation (2.9), these solutions, which are stationary for
a fixed $\omega$, correspond to uniform rotation around the polar axis for
other values of $\omega$.

Geometrically, (2.10) means that the gradients of the stream function $\psi$
and of the potential vorticity $\Delta\psi+2\omega\sin\theta$ are parallel.
Since the gradient is orthogonal to the level set, in regions of
${\mathbb{S}}^{2}$ where $\text{grad}\,\psi\neq(0,0)$ the rank theorem ensures
that (2.10) is locally equivalent to the elliptic problem (1.1), namely

$\Delta\psi+2\omega\sin\theta=F(\psi)$

holds for some $C^{1}$-function $F$; see the discussion in [16]. It is easy to
check that any solution of the elliptic problem (1.1) on ${\mathbb{S}}^{2}$
will also solve (2.10), but the converse is not true in general. For example,
any zonal function $\psi$ solves (2.10) but does not have to be a solution of
(1.1), as shown by the case of constant functions.

Two classes of explicit solutions of (2.10) are known:

  * •

zonal solutions $\psi(\theta)$;

  * •

Rossby-Haurwitz waves of the form

(2.11)
$\psi(\varphi,\theta)=\tfrac{2\omega}{2-j(j+1)}\,\sin\theta+\beta\,Y(\varphi,\theta)\,,\qquad
j\geq 2\,,\quad\beta\in{\mathbb{R}}\,,\quad Y\in E_{j}\,,$

where $E_{j}$ is the $(2j+1)$-dimensional eigenspace of the Laplace-Beltrami
operator associated to the eigenvalue $-j(j+1)$. These are solutions of (2.10)
for $F(s)=-j(j+1)s$.

It is easy to verify that functions of this type solve (2.10). Using the
symmetry (2.9), one obtains explicit non-trivial travelling-wave solutions of
(${\mathcal{E}}_{\omega}$) of the form

(2.12) $\psi(\varphi-ct,\theta)=\alpha\,\sin\theta+\beta\,Y(\varphi-
ct,\theta)\quad\text{with}\quad\begin{cases}&Y\in E_{j}\ (j\geq
1)\,,\qquad\quad\alpha\in{\mathbb{R}}\,,\quad\beta\in{\mathbb{R}}\setminus\\{0\\}\,,\\\
&c=\frac{2\omega}{j(j+1)}+\alpha\,\frac{j(j+1)-2}{j(j+1)}\,,\end{cases}$

as one can easily check. Two particular cases are of great interest:

  * •

for $\alpha=\tfrac{2\omega}{2-j(j+1)}$ with $j\geq 2$ we obtain the stationary
waves (2.11) with wave speed $c=0$;

  * •

for $\alpha=\omega$ we obtain geostationary waves that propagate azimuthally
with wave speed $c=\omega$, which can be subsumed into 3D stratospheric flows
(see Section 6).

Some other explicit solutions are listed in the next section and an approach
providing us with further classes of (non-explicit) stationary solutions is
provided in Section 6.

##  3\. Rigidity results

Due to the important role played by (1.1) in the quest of stationary solutions
of the vorticity equation (${\mathcal{E}}_{\omega}$), and given the
accessibility of symmetries for ${\mathbb{S}}^{2}$, it is natural to expect
some concrete answers to the question about solutions capturing symmetries of
the spherical domain.

######  Theorem 4.

Consider classical solutions $\psi$ of (1.1), for a continuously
differentiable function $F:{\mathbb{R}}\to{\mathbb{R}}$.

  * (i)

If $\omega=0$ and $F^{\prime}>-2$, then $\psi$ is a constant. This statement
is optimal since the non-constant elements of the first eigenspace
$\mathbb{E}_{1}$ of the Laplace-Beltrami operator satisfy $\Delta\psi=-2\psi$.

  * (ii)

If $\omega=0$ and $F^{\prime}>-6$, then $\psi$ is zonal up to a transformation
in $\mathbb{SO}(3)$. This statement is optimal since the non-zonal elements of
the second eigenspace $\mathbb{E}_{2}$ of the Laplace-Beltrami operator
satisfy $\Delta\psi=-6\psi$.

  * (iii)

For $\omega\neq 0$ there exists $\alpha\in\mathbb{R}$ such that
$\mathbb{P}_{2}\psi=\alpha\sin\theta$.

  * (iv)

If $\omega\neq 0$ and $F^{\prime}>-6$, then $\psi$ is zonal up to a
transformation in $\mathbb{SO}(3)$. This statement is optimal since the
Rossby-Haurwitz wave (2.11) with $j=2$ and $Y$ a non-trivial linear
combination of all spherical harmonics $\\{Y_{2}^{m}\\}_{|m|\leq 2}$ satisfies
$\Delta\psi=-6\psi$.

###### Proof.

(i) Pairing (1.1) for $\omega=0$ with $\Delta\psi$ and integrating by parts
leads to the identity

$\iint_{{\mathbb{S}}^{2}}|\Delta\psi|^{2}\,{\rm
d}\sigma=-\iint_{{\mathbb{S}}^{2}}F^{\prime}(\psi)|\operatorname{grad}\psi|^{2}\,{\rm
d}\sigma\,.$

Assuming that $F^{\prime}>2$, the fact that the smallest nonzero eigenvalue of
$-\Delta$ is $2$ gives

$\iint_{{\mathbb{S}}^{2}}|\Delta\psi|^{2}\,{\rm
d}\sigma<2\iint_{{\mathbb{S}}^{2}}|\operatorname{grad}\psi|^{2}\,{\rm
d}\sigma\leq\iint_{{\mathbb{S}}^{2}}|\Delta\psi|^{2}\,{\rm d}\sigma\,$

which is a contradiction unless $\Delta\psi=0$, which forces
$|\operatorname{grad}\psi|=0$.

(ii) Up to the action of $\mathbb{O}(3)$, we can assume that
$\mathbb{P}_{2}\psi=\alpha\sin\theta$ for some $\alpha\in\mathbb{R}$.
Differentiating (1.1) for $\omega=0$ with respect to $\varphi$, and pairing
the result with $\partial_{\varphi}\psi$ leads to the identity

$\iint_{{\mathbb{S}}^{2}}|\nabla\partial_{\varphi}\psi|^{2}\,{\rm
d}\sigma=-\iint_{{\mathbb{S}}^{2}}F^{\prime}(\psi)|\partial_{\varphi}\psi|^{2}\,{\rm
d}\sigma\,.$

Under the assumption that $F^{\prime}>-6$, this implies that

$\iint_{{\mathbb{S}}^{2}}|\nabla\partial_{\varphi}\psi|^{2}\,{\rm
d}\sigma<6\iint_{{\mathbb{S}}^{2}}|\partial_{\varphi}\psi|^{2}\,{\rm
d}\sigma\,,$

unless $\partial_{\varphi}\psi\equiv 0$. Since
$\partial_{\varphi}\mathbb{P}_{2}\psi=0$, and $\partial_{\varphi}$ commutes
with the spectral projectors of the Laplacian, this gives
$\mathbb{P}_{2}\partial_{\varphi}\psi=0$. Therefore, we can use the fact that
the third eigenvalue of $-\Delta$ is $6$ to conclude that

$\iint_{{\mathbb{S}}^{2}}|\Delta\partial_{\varphi}\psi|^{2}\,{\rm
d}\sigma<6\iint_{{\mathbb{S}}^{2}}|\partial_{\varphi}\psi|^{2}\,dx\leq\int|\Delta\partial_{\varphi}\psi|^{2}\,{\rm
d}\sigma\,,$

which is a contradiction unless $\partial_{\varphi}\psi\equiv 0$.

(iii) If $\psi$ solves (1.1), then

$\left[-\partial_{\theta}\psi\partial_{\varphi}+\partial_{\varphi}\psi\partial_{\theta}\right]\left(\Delta\psi+2\omega\sin\theta\right)=0.$

Multiplying the above by by ${\rm e}^{{\rm i}\varphi}$ and integrating the
result on the sphere, we get

$0=\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\int_{0}^{2\pi}{\rm e}^{{\rm
i}\varphi}(-\partial_{\theta}\psi\,\partial_{\varphi}+\partial_{\varphi}\psi\,\partial_{\theta})\left[\partial_{\theta}^{2}\psi-\tan\theta\,\partial_{\theta}\psi+\frac{1}{(\cos\theta)^{2}}\,\partial_{\varphi}^{2}\psi+2\omega\sin\theta\right]\cos\theta\,{\rm
d\varphi}{\rm d}\theta=I+II+III+IV,$

where $I$, $II$, $III$ and $IV$ correspond to the four terms in the bracketed
expression. Integrating by parts repeatedly, and using the fact that all
boundary terms vanish since $\cos\theta=\partial_{\varphi}\psi=0$ at
$\theta=\pm\frac{\pi}{2}$, we obtain

$\displaystyle I=\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\int_{0}^{2\pi}{\rm
e}^{{\rm i}\varphi}\Big{(}{\rm
i}\sin\theta(\partial_{\theta}\psi)^{2}-\cos\theta\,\partial_{\varphi}\psi\,\partial_{\theta}\psi\Big{)}\,{\rm
d\varphi}{\rm d}\theta\,,$ $\displaystyle
II=\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\int_{0}^{2\pi}{\rm e}^{{\rm
i}\varphi}\left(-{\rm
i}\sin\theta(\partial_{\theta}\psi)^{2}-\frac{(\sin\theta)^{2}}{\cos\theta}\,\partial_{\varphi}\psi\,\partial_{\theta}\psi\right)\,{\rm
d\varphi}{\rm d}\theta$ $\displaystyle
III=\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\int_{0}^{2\pi}{\rm e}^{{\rm
i}\varphi}\frac{1}{\cos\theta}\,\partial_{\varphi}\psi\,\partial_{\theta}\psi\,{\rm
d\varphi}{\rm d}\theta\,,$

from which it follows that

$I+II+III=0\,.$

Turning to $IV$, repeated integration by parts show that

$IV=-2{\rm
i}\omega\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\int_{0}^{2\pi}\psi{\rm e}^{{\rm
i}\varphi}(\cos\theta)^{2}\,{\rm d\varphi}{\rm d}\theta\,.$

Furthermore the above implies that $IV=0$, or in other words

$\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\int_{0}^{2\pi}\psi(\cos\theta\sin\varphi)\cos\theta\,{\rm
d\varphi}{\rm
d}\theta=\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\int_{0}^{2\pi}\psi(\cos\theta\cos\varphi)\cos\theta\,{\rm
d\varphi}{\rm d}\theta\,,$

which is the desired result.

(iv) Differentiating (1.1) with respect to $\varphi$, and pairing the result
with $\partial_{\varphi}\psi$ leads to the identity

$\iint_{{\mathbb{S}}^{2}}|\nabla\partial_{\varphi}\psi|^{2}\,{\rm
d}\sigma=-\iint_{{\mathbb{S}}^{2}}F^{\prime}(\psi)|\partial_{\varphi}\psi|^{2}\,{\rm
d}\sigma\,.$

Assuming that $F^{\prime}>-6$, this implies that

$\iint_{{\mathbb{S}}^{2}}|\nabla\partial_{\varphi}\psi|^{2}\,{\rm
d}\sigma<6\iint_{{\mathbb{S}}^{2}}|\partial_{\varphi}\psi|^{2}\,{\rm
d}\sigma\,.$

However, we know from (iii) that $\partial_{\varphi}\mathbb{P}_{2}\psi=0$.
Since $\partial_{\varphi}$ commutes with the spectral projectors of the
Laplacian, this gives $\mathbb{P}_{2}\partial_{\varphi}\psi=0$. Therefore, the
third eigenvalue of $-\Delta$ being $6$, we conclude that

$\iint_{{\mathbb{S}}^{2}}|\Delta\partial_{\varphi}\psi|^{2}\,{\rm
d}\sigma<6\iint_{{\mathbb{S}}^{2}}|\partial_{\varphi}\psi|^{2}\,{\rm
d}\sigma\leq\iint_{{\mathbb{S}}^{2}}|\Delta\partial_{\varphi}\psi|^{2}\,{\rm
d}\sigma\,,$

which is a contradiction unless $\partial_{\varphi}\psi\equiv 0$. ∎

Finding explicit smooth solutions of equation (1.1) for nonlinear functions
$F$ is far from obvious. The case

$F(s)=a\,{\rm e}^{bs}+\frac{2}{b}\qquad\text{with}\quad
a\,,b\in{\mathbb{R}}\setminus\\{0\\}\,,$

is of great interest in statistical mechanics and Riemannian geometry (see
[12, 13, 33]), with the specific value of the additive constant (given the
exponential term) determined by the Gauss constraint. For $ab>0$ there are no
smooth solutions, while for $ab<0$ the general solution is available (see
[19]), and can be expressed in terms of a single analytic function by
performing the stereographic projection of the 2-sphere ${\mathbb{S}}^{2}$
which maps the North Pole to infinity and the South Pole to the origin of the
compactifed complex plane ${\mathbb{C}}\cup\\{\infty\\}$. Unfortunately,
singularities are ubiquituous222A notable exception being the solution
$\psi_{0}(\varphi,\theta)=\frac{2}{b}\,\ln(1+|\zeta|^{2})+\frac{1}{b}\,\ln\Big{(}\frac{2\alpha^{2}}{-ab(1+|\alpha\zeta+\beta|^{2})^{2}}\Big{)}$
with free parameters $\alpha\in{\mathbb{R}}\setminus\\{0\\}$ and
$\beta\in{\mathbb{C}}$, where
$\zeta=\tan\big{(}\frac{\theta}{2}+\frac{\pi}{4}\big{)}\,{\rm e}^{{\rm
i}\varphi}\in{\mathbb{C}}\cup\\{\infty\\}$. This solution is generated by the
polynomial $\zeta\mapsto\alpha\zeta+\beta$ and all streamlines are circles
(the streamline pattern being identical to that depicted in Fig. 3). A
suitable rotation in ${\mathbb{S}O}(3))$ transforms this solution into a zonal
one., so that this type of nonlinearity is not suitable for regular global
flows but rather for regular flows in regions bounded by a streamline.

Figure 3. Depiction of the streamlines for the solution (3.5).

Theorem 4 permits us to construct nonlinearities $F$ for which the equation
(1.1) admits classical non-zonal solutions. Denoting for $\alpha\in(0,1)$ and
$k\in{\mathbb{N}}\cup\\{0\\}$ by $C^{k,\alpha}({\mathbb{S}}^{2})$ the Banach
space of $k$-times Hölder continuously differentiable functions
$\psi:{\mathbb{S}}^{2}\to{\mathbb{R}}$ with exponent $\alpha$, we can write
(1.1) for $\omega=0$ and any continuously differentiable function
$F:{\mathbb{R}}\to{\mathbb{R}}$ in the form

(3.1) ${\mathcal{F}}(\psi)=0\quad\text{for the nonlinear
operator}\quad{\mathcal{F}}:C^{2,\alpha}({\mathbb{S}}^{2})\to
C^{0,\alpha}({\mathbb{S}}^{2})\,,\qquad{\mathcal{F}}(\psi)=\Delta\psi-F(\psi)\,.$

Note that (3.1) is ${\mathbb{O}}(3)$-equivariant:

${\mathcal{F}}(G\psi)=G{\mathcal{F}}(\psi)\quad\text{for all}\quad
G\in{\mathbb{O}}(3)\quad\text{and all}\quad\psi\in
C^{2,\alpha}({\mathbb{S}}^{2})\,,$

the natural action of ${\mathbb{O}}(3)$ on functions
$\psi:{\mathbb{S}}^{2}\to{\mathbb{R}}$ being given by

$(G\psi)(X):=\psi(G^{T}X)\,,\qquad G\in{\mathbb{O}}(3)\,,$

where $G^{T}$ is the transpose of the matrix $G$ in the matrix representation
of the orthogonal group ${\mathbb{O}}(3)$. Therefore, if $\psi\in
C^{2,\alpha}({\mathbb{S}}^{2})$ solves (3.1), then so does $G\psi$ for any
$G\in{\mathbb{O}}(3)$. In particular, since the zonal solutions of (3.1) are
those symmetric with respect to the polar axis, if $G\in{\mathbb{O}}(3)$
breaks this symmetry, then $G\psi$ is a non-zonal solution of (3.1) whenever
$\psi$ is a non-constant zonal solution. We may take $G$ to be the rotation
that transforms rotations about the polar axis into rotations about a fixed
horizontal axis. By computing the expression

$\frac{{\rm d}^{2}f}{{\rm d}\theta^{2}}-\tan(\theta)\,\frac{{\rm d}f}{{\rm
d}\theta}$

for a suitable function of the variable $\sin\theta$, we obtain specific zonal
solutions of (3.1). For example, one can check that for every $\varepsilon>0$
the function

$\psi_{0}(\theta)=\ln\Big{(}\frac{1+\varepsilon\sin\theta}{1-\varepsilon\sin\theta}\Big{)}\,,\qquad-\frac{\pi}{2}\leq\theta\leq\frac{\pi}{2}\,,$

is a zonal solution of (1.1) with $\omega=0$, for

(3.2) $F(\psi)=-\frac{1-\varepsilon^{2}}{2}\,[2\sinh(\psi)+\sinh(2\psi)]\,.$

Consequently,

(3.3)
$\psi(\varphi,\theta)=\ln[1+\varepsilon\cos^{2}(\theta)\sin^{2}(\varphi-\varphi_{0})]$

is a non-zonal solution of (3.1) for every fixed $\varphi_{0}\in[0,2\pi)$.
While the solution (3.3) is known (see [1]), the above symmetry approach
rather than the ad-hoc Ansatz made in [1] explains how it comes about and also
provides a procedure leading to the construction of other explicit non-zonal
solutions. Indeed, the zonal solution

$\psi_{0}(\theta)={\rm
e}^{\varepsilon\sin\theta}-1\,,\qquad-\frac{\pi}{2}\leq\theta\leq\frac{\pi}{2}\,,$

of (3.1) for

(3.4)
$F(\psi)=\varepsilon^{2}(1+\psi)-(1+\psi)\ln^{2}(1+\psi)-2(1+\psi)\ln(1+\psi)\,.$

leads to the non-zonal solution

(3.5) $\psi(\varphi,\theta)={\rm
e}^{\varepsilon\cos\theta\sin(\varphi-\varphi_{0})}-1$

for every fixed $\varphi_{0}\in[0,2\pi)$. Note that the streamlines of the
solutions (3.3)-(3.5) are circles with collinear centres along a segment lying
in the equatorial plane and passing through the centre of the sphere (obtained
by suitably rotating Figure 3), since by passing to spherical coordinates in
${\mathbb{R}}^{3}$ we see that the level sets
$[\cos\theta\sin(\varphi-\varphi_{0})=d]$ are precisely the points on the
sphere at distance $|d|$ from the line obtained rotating the $x$-axis by
$\varphi_{0}$-degrees in the equatorial plane.

##  4\. Stability of stationary solutions

To gain insight into the complicated dynamics of the atmosphere it is
advantageous to regard it as a background zonal flow presenting fluctuations,
some transient or at least unable to develop significantly while other are
enhanced in time (instabilities). One distinguishes between barotropic
instability, in which perturbations draw strength at the expense of the
kinetic energy of the background flow, and baroclinic instability, fed by the
release potential energy through the lifting (sinking) of warm, relatively
light (respectively of cold, relatively dense) fluid. Barotropic flows are
prevalent in the stratosphere of the outer planets, so that the investigation
of the stability of the explicitly known stationary solutions of
(${\mathcal{E}}_{\omega}$) is of great physical relevance.

###  4.1. Linear stability of zonal flows

The giant planet atmospheres feature alternating prograde (eastward) and
retrograde (westward) jets of different speeds and widths (see Fig. 2).
Moreover, observations of the zonal flow patterns of Jupiter and Saturn
indicate the development of eddies around the peaks of the westward jets. This
type of instability is reminiscent of Rayleigh stability criterion for shear
flows.

Consider a zonal flow with a stream function $\psi_{0}(\theta)$ of class
$C^{3}$, representing a stationary solution of the vorticity equation
(${\mathcal{E}}_{\omega}$). To study perturbations of this zonal flow, it is
convenient to use the variable $s=\sin\theta\in[-1,1]$ instead of the latitude
$\theta$. Note that the conjugate variable to the longitude $\varphi$ is not
the latitude $\theta$ but the axial component in Cartesian coordinates,
$s=\sin\theta$, for which Hamilton’s equations hold:

$\begin{cases}D_{t}\varphi=-\Psi_{s}\,,\\\
D_{t}s=\Psi_{\varphi}\,.\end{cases}$

Setting

(4.1) $\Psi(\varphi,s)=\psi(\varphi,\theta)\,,$

the expressions for the velocity and vorticity are

(4.2) $u=-\sqrt{1-s^{2}}\,\partial_{s}\Psi\,,\qquad
v=\frac{1}{\sqrt{1-s^{2}}}\,\partial_{\varphi}\Psi\,,\qquad\Omega=(1-s^{2})\,\partial_{s}^{2}\Psi-2s\,\partial_{s}\Psi+\frac{1}{1-s^{2}}\,\partial_{\varphi}^{2}\Psi:=\Upsilon(\varphi,s,t)\,,$

and (${\mathcal{E}}_{\omega}$) takes the form

(4.3)
$\partial_{t}\Upsilon+\partial_{\varphi}\Psi\,(\partial_{s}\Upsilon+2\omega)-\partial_{s}\Psi\,\partial_{\varphi}\Upsilon=0\,.$

The linearization of the vorticity equation (4.3) around the zonal flow with
stream function $\Psi_{0}(s)=\psi_{0}(\theta)$ and vorticity $\Upsilon_{0}$ is
the equation

(4.4)
$\partial_{t}\Upsilon-\Psi_{0}^{\prime}\,\partial_{\varphi}\Upsilon+(\Upsilon_{0}^{\prime}+2\omega)\,\partial_{\varphi}\Psi=0$

for the infinitesimal perturbation $\Psi(\varphi,s,t)$ with vorticity
$\Upsilon$, where $\Psi_{0}^{\prime}=\partial_{s}\Psi_{0}$ and
$\Upsilon_{0}^{\prime}=\partial_{s}\Upsilon_{0}$. In accordance to the
discussion in the preamble of relation (2.7), the boundary conditions
associated to (4.4) are

(4.5) $\partial_{\varphi}\Psi=0\quad\text{at}\quad s=\pm 1\,.$

We can write (4.4) in the form

(4.6) $\partial_{t}\Upsilon={\mathcal{L}}_{0}\Upsilon\,,$

with the linear operator

${\mathcal{L}}_{0}=\Psi_{0}^{\prime}\,\partial_{\varphi}-(\Upsilon_{0}^{\prime}+2\omega)\,\partial_{\varphi}\Delta^{-1}$

acting in the space

$L_{0}^{2}=\\{\Upsilon\in
L^{2}({\mathbb{S}}^{2}):\,\int_{0}^{2\pi}\int_{-1}^{1}\Upsilon(\varphi,s)\,{\rm
d}s{\rm d}\varphi=0\\}$

of vorticities subject to the Gauss constraint (2.6). Since the operator
${\mathcal{L}}_{0}$ commutes with translations in the azimuthal direction, by
means of the Fourier modes

$f(\varphi,s,t)=\sum_{k\in{\mathbb{Z}}\setminus\\{0\\}}f_{k}(s,t)\,{\rm
e}^{{\rm i}k\varphi}\,,$

for $f=\Delta^{-1}\Upsilon$, we can decompose (4.6) into

$\partial_{t}\Upsilon_{k}={\rm i}k\,{\mathcal{L}}_{0}^{k}\Upsilon_{k}\,,\qquad
k\in{\mathbb{Z}}\setminus\\{0\\}\,,$

where the operators

${\mathcal{L}}_{0}^{k}=\Psi_{0}^{\prime}-(\Upsilon_{0}^{\prime}+2\omega)\,\Delta_{k}^{-1}$

act in $L^{2}[-1,1]$ subject to the boundary conditions $\mu(s)=0$ at $s=\pm
1$ for $\mu=\Delta_{k}^{-1}\Upsilon_{k}$, which ensure that

$\Upsilon_{k}=\Delta_{k}\mu=(1-s^{2})\mu^{\prime\prime}-2s\mu^{\prime}-\tfrac{k^{2}}{1-s^{2}}\,\mu\,,\qquad
k\in{\mathbb{Z}}\setminus\\{0\\}\,,$

is regular at $s=\pm 1$. Since for every $k\in{\mathbb{Z}}\setminus\\{0\\}$
the operator ${\mathcal{L}}_{0}^{k}$ is a compact perturbation of the
multiplication operator $\Psi_{0}^{\prime}$ with the purely essential spectrum

$\Sigma_{0}=\Big{[}\inf_{s\in[-1,1]}\\{\Psi_{0}^{\prime}(s)\\}\,,\,\sup_{s\in[-1,1]}\\{\Psi_{0}^{\prime}(s)\\}\Big{]}\,,$

the essential spectrum of ${\mathcal{L}}_{0}^{k}$ coincides with the closed
real interval $\Sigma_{0}$, and the rest of the spectrum of
${\mathcal{L}}_{0}^{k}$ consists of at most countably many isolated
eigenvalues of finite multiplicity. The discrete spectrum of
${\mathcal{L}}_{0}^{k}$ is symmetric about the real axis, since the complex
conjugate of an eigenfunction for the eigenvalue
$\lambda\in{\mathbb{C}}\setminus{\mathbb{R}}$ is an eigenfunction for
$\overline{\lambda}$. The Fourier mode decomposition thus yields the spectrum
of the operator ${\mathcal{L}}_{0}$: the essential spectrum

$\\{\lambda\in{\mathbb{C}}:\ \lambda={\rm i}kr\quad\text{with}\quad
k\in{\mathbb{Z}}\setminus\\{0\\}\quad\text{and}\quad r\in\Sigma_{0}\\}$

is located on the imaginary axis, and the discrete spectrum is symmetric about
the imaginary axis and comprises at most countable many isolated eigenvalues
of finite multiplicity. Linear instability thus amounts to the existence of an
eigenvalue of ${\mathcal{L}}_{0}$ with non-zero real part. Indeed, due to the
symmetry of the discrete spectrum about the imaginary axis, this means the
existence of an eigenvalue $\xi={\rm i}k\lambda$ with strictly positive real
part, for some $k\in{\mathbb{Z}}\setminus\\{0\\}$ and some eigenvalue
$\lambda\not\in\Sigma_{0}$ of ${\mathcal{L}}_{0}^{k}$. If $\Upsilon_{k}(s)$ is
the corresponding eigenfunction of ${\mathcal{L}}_{0}^{k}$, then
$\Upsilon_{k}(s)\,{\rm e}^{{\rm i}k(\varphi+\lambda t)}$ solves (4.6) and its
amplitude grows indefinitely for $t\to\infty$.

The spectral analysis is very challenging, and only few thorough
investigations appear to be available:

  * •

For $\psi_{0}(\theta)=\alpha\sin\theta$ with
$\alpha\in{\mathbb{R}}\setminus\\{0\\}$, the essential spectrum of
${\mathcal{L}}_{0}^{k}=\alpha\partial_{\varphi}+2(\omega-\alpha)\partial_{\varphi}\Delta_{k}^{-1}$
is the single point $\alpha$ and by expanding in spherical harmonics we see
that the discrete spectrum consists of the real eigenvalues
$\Big{\\{}\Big{(}\alpha-\frac{2(\alpha-\omega)}{|k|(|k|+1)}\Big{)}\Big{\\}}$
of multplicity $2|k|$. This zonal flow is thus linearly stable.

  * •

The zonal flows

(4.7) $\psi_{0}(\theta)=\alpha Y_{1}^{0}(\theta)+\beta Y_{2}^{0}(\theta)\,,$

where $\alpha,\,\beta\in{\mathbb{R}}\setminus\\{0\\}$ and
$Y_{1}^{0}(\theta)=\sqrt{\tfrac{3}{4\pi}}\,\sin\theta$,
$Y_{2}^{0}(\theta)=\sqrt{\tfrac{5}{16\pi}}\,(3\sin^{2}\theta-1)$, are the
zonal spherical harmonics of degree $1$ and $2$ (see the Appendix), turn out
to be linearly stable (see [48, 52]).

  * •

The zonal flows

(4.8) $\psi_{0}(\theta)=\beta
Y_{3}^{0}(\theta)\,,\qquad\beta\in{\mathbb{R}}\setminus\\{0\\}\,,$

with three jets (defined as extrema of the azimuthal velocity,
$-\partial_{\theta}\psi_{0}$, and equal to the number of nodes of the zonal
spherical harmonics $Y_{3}^{0}$ of degree $3$) are known (see [48, 52]) to be
linearly unstable if their amplitude $|\beta|$ exceeds a critical value
$\beta^{\ast}(\omega)$.

  * •

Numerical simulations (see [8]) indicate that the zonal flows

(4.9) $\psi_{0}(\theta)=\beta
Y_{j}^{0}(\theta)\,,\qquad\beta\in{\mathbb{R}}\setminus\\{0\\}\,,$

with $j>3$ jets are linearly unstable if their amplitude $|\beta|$ is
sufficiently large.

Let us now briefly describe a classical approach – typically pursued within
the setting of flat-space geometry (see [42]) – that provides insight into the
linear stability without actually finding eigenvalues. Motivated by the
considerations made above, we consider normal mode perturbations of zonal
flows of the form

(4.10) $\Psi(\varphi,s,t)=\Phi(s)\,{\rm e}^{{\rm i}k(\varphi-\mu t)}\,,$

subject to the boundary condition (4.5), where
$k\in{\mathbb{Z}}\setminus\\{0\\}$ is the wave number and $\mu\in{\mathbb{C}}$
is the wave speed. The zonal flow with stream function $\Psi_{0}$ is linearly
unstable if there exists a nontrivial solution (4.10) to (4.4) and (4.5) with
$\mathfrak{Im}(k\mu)>0$, since in this case the amplitude of the perturbation
(4.10) grows indefinitely with time. For $\mathfrak{Im}(\mu)\neq 0$ the
equation (4.4) reduces to

(4.11)
$\Big{(}(1-s^{2})\Phi^{\prime}\Big{)}^{\prime}-\Big{\\{}\frac{k^{2}}{1-s^{2}}+\frac{\Upsilon_{0}^{\prime}+2\omega}{\Psi_{0}^{\prime}+\mu}\Big{\\}}\,\Phi=0\,,\qquad
s\in(-1,1)\,,$

for the amplitude $\Phi$. Multiplying (4.11) by the complex-conjugate
$\overline{\Phi}$ and integrating the outcome on $[-1,1]$, an integration by
parts yields

(4.12)
$\int_{-1}^{1}\Big{\\{}(1-s^{2})|\Phi^{\prime}(s)|^{2}+\frac{k^{2}}{1-s^{2}}\,|\Phi(s)|^{2}\Big{\\}}\,{\rm
d}\theta+\int_{-1}^{1}\frac{\Big{(}\Upsilon_{0}^{\prime}(s)+2\omega\Big{)}\Big{(}\Psi_{0}^{\prime}(s)+\overline{\mu}\Big{)}}{|\Psi_{0}^{\prime}(s)+\mu|^{2}}\,|\Phi(s)|^{2}\,{\rm
d}s=0\,.$

Since the imaginary part of (4.12) is

$-\mathfrak{Im}(\mu)\,\int_{-1}^{1}\frac{\Upsilon_{0}^{\prime}(s)+2\omega}{|\Psi_{0}^{\prime}(s)+\mu|^{2}}\,|\Phi(s)|^{2}\,{\rm
d}s\,,$

we obtain the Rayleigh necessary condition for linear instability

(4.13)
$\int_{-1}^{1}\frac{\Upsilon_{0}^{\prime}(s)+2\omega}{|\Psi_{0}^{\prime}(s)+\mu|^{2}}\,|\Phi(s)|^{2}\,{\rm
d}s=0\,,$

which requires $\Big{(}\Upsilon_{0}^{\prime}(s)+2\omega\Big{)}$ to change sign
on $(-1,1)$. A further necessary condition, due to Fjortoft, is obtained by
taking also into account the real part of (4.12): if (4.13) holds, then (4.12)
yields

(4.14)
$\int_{-1}^{1}\frac{\Big{(}\Upsilon_{0}^{\prime}(s)+2\omega\Big{)}\Psi_{0}^{\prime}(s)}{|\Psi_{0}^{\prime}(s)+\mu|^{2}}\,|\Phi(s)|^{2}\,{\rm
d}s=-\int_{-1}^{1}\Big{\\{}(1-s^{2})|\Phi^{\prime}(s)|^{2}+\frac{k^{2}}{1-s^{2}}\,|\Phi(s)|^{2}\Big{\\}}\,{\rm
d}\theta<0\,.$

For any $\gamma\in{\mathbb{R}}$, adding (4.13) multiplied by $\gamma$ to
(4.14), we get

$\int_{-1}^{1}\frac{\Big{(}\Upsilon_{0}^{\prime}(s)+2\omega\Big{)}\Big{(}\Psi_{0}^{\prime}(s)+\gamma\Big{)}}{|\Psi_{0}^{\prime}(s)+\mu|^{2}}\,|\Phi(s)|^{2}\,{\rm
d}s=-\int_{-1}^{1}\Big{\\{}(1-s^{2})|\Phi^{\prime}(s)|^{2}+\frac{k^{2}}{1-s^{2}}\,|\Phi(s)|^{2}\Big{\\}}\,{\rm
d}\theta<0\,.$

Consequently, we must have
$\Big{(}\Upsilon_{0}^{\prime}(s_{0})+2\omega\Big{)}\Big{(}\Psi_{0}^{\prime}(s_{0})+\gamma\Big{)}<0$
at some $s_{0}\in(-1,1)$. In terms of latitude, these necessary conditions for
the linear instability of a zonal flow with azimuthal velocity $U_{0}$ and
vorticity $\Omega_{0}$ read:

  * •

(Rayleigh’s criterion) the meridional gradient of the total vorticity of the
zonal flow,

$\Omega_{0}^{\prime}(\theta)+2\omega\cos\theta=\Big{(}\Upsilon_{0}^{\prime}(\sin\theta)+2\omega\Big{)}\cos\theta\,,$

changes sign on the interval $\big{(}-\tfrac{\pi}{2},\tfrac{\pi}{2}\big{)}$;

  * •

(Fjortoft’s criterion) for every $\gamma\in{\mathbb{R}}$ we must have

$\Big{(}\Omega_{0}^{\prime}(\theta)+2\omega\cos\theta\Big{)}\Big{(}\frac{U_{0}(\theta)}{\cos\theta}-\gamma\Big{)}>0$

at some $\theta_{0}\in\big{(}-\tfrac{\pi}{2},\tfrac{\pi}{2}\big{)}$.

Note that Fjortoft’s criterion implies that of Rayleigh since the existence of
real numbers $\gamma_{1}>\gamma_{2}$ with

$\frac{U_{0}(\theta)}{\cos\theta}-\gamma_{1}<0<\frac{U_{0}(\theta)}{\cos\theta}-\gamma_{2}\,,\qquad\theta\in\big{(}-\tfrac{\pi}{2},\tfrac{\pi}{2}\big{)}\,,$

combined with Fjortoft’s constraint for $\gamma=\gamma_{1}$ and
$\gamma=\gamma_{2}$ ensures that the expression
$\Omega_{0}^{\prime}(\theta)+2\omega\cos\theta$ has opposite signs at two
points in $\big{(}-\tfrac{\pi}{2},\tfrac{\pi}{2}\big{)}$. Both criteria fail
for the linearly stable zonal flow $\psi_{0}(\theta)=\alpha\theta$ with
$\alpha\in{\mathbb{R}}\setminus\\{0\\}$ since in this case

$\Omega_{0}^{\prime}(\theta)+2\omega\cos\theta=2(\omega-\alpha)\cos\theta\,,\qquad\frac{U_{0}(\theta)}{\cos\theta}-\gamma=-\alpha-\gamma\,.$

However, they are generally far from sufficient to ensure linear instability:
e.g., both hold for the linearly stable zonal flow
$\psi_{0}(\theta)=\frac{\omega}{3}\,\sin^{2}\theta$, with

$\Omega_{0}^{\prime}(\theta)+2\omega\cos\theta=2\omega\cos\theta(1-2\sin\theta)\,,\qquad\frac{U_{0}(\theta)}{\cos\theta}-\gamma=-\frac{2\omega}{3}\,\Big{(}\sin\theta+\frac{3\gamma}{2\omega}\Big{)}\,.$

Rayleigh’s criterion appears not to be sufficient for linear instability: it
holds for the latitudinal profile of the persistent zonal jets on Jupiter and
Saturn (see the data in [45]).

###  4.2. Nonlinear stability of stationary solutions: the Arnold approach

Due to the intricate nature of the investigation of linear stability and to
its rather inconclusive nature333The question whether linearization is
conclusive with respect to the stability of solutions to the nonlinear
vorticity equation is open., we pursue the issue of nonlinear stability
directly, relying on linear results (whenever available) to make an informed
guess.

We consider a solution $\psi_{0}$ of (2.10), which is thus a stationary
solution of (${\mathcal{E}}_{\omega}$). We implement Arnold’s method (see [3,
4]) by defining the functional

$\mathcal{E}(\psi)=\iint_{{\mathbb{S}}^{2}}\left[\frac{1}{2}|U_{0}+U|^{2}+K(\Omega_{0}+\Omega+2\omega\sin\theta)+A\sin\theta\Omega+B(\mathbb{P}_{1}\psi)^{2}\right]\,{\rm
d}\sigma\,,$

which is a sum of conserved quantities for the flow of
(${\mathcal{E}}_{\omega}$). Our aim is to choose the function $K$ and the
constants $A$ and $B$ so that $\psi=0$ is a critical point of $\mathcal{E}$,
and, furthermore, the second variation of $\mathcal{E}$ at $0$ is a definite
quadratic form.

With this in mind, we now compute the first variation of $\mathcal{E}$ at $0$:

$\displaystyle d\mathcal{E}_{0}(\delta{\psi})$
$\displaystyle=\iint_{{\mathbb{S}}^{2}}\left[U_{0}\cdot\delta
U+K^{\prime}(\Omega_{0}+2\omega\sin\theta)\delta\Omega+A\sin\theta\,\delta\Omega\right]\,{\rm
d}\sigma\,.$

Expressing $\delta\Omega$ in terms of $\delta U$ and integrating by parts
using that $\operatorname{grad}\sin\theta=\cos\theta\,\mathbf{e}_{\theta}$,
this becomes

$\displaystyle d\mathcal{E}_{{\psi}}(\delta{\psi})$
$\displaystyle=\iint_{{\mathbb{S}}^{2}}\left[U_{0}\cdot\delta
U-K^{\prime}(\Omega_{0}+2\omega\sin\theta)\operatorname{div}J\delta
U-A\sin\theta\operatorname{div}J\delta U\right]\,{\rm d}\sigma$
$\displaystyle=\iint_{{\mathbb{S}}^{2}}\left[U_{0}\cdot\delta
U-K^{\prime\prime}(\Omega_{0}+2\omega\sin\theta)J\operatorname{grad}(\Omega_{0}+2\omega\sin\theta)\cdot\delta
U+A\cos\theta\,\mathbf{e}_{\varphi}\cdot\delta U\right]\,{\rm d}\sigma\,.$

The stationary solution $\psi_{0}$ being such that
$\Omega_{0}+2\omega\sin\theta=F(\psi_{0})$ implies, after taking the gradient,
that

$\operatorname{grad}(\Omega_{0}+2\omega\sin\theta)=F^{\prime}(\psi_{0})\operatorname{grad}\psi_{0}=-F^{\prime}(\psi_{0})(JU_{0}+\omega\cos\theta\,\mathbf{e}_{\theta}).$

Therefore

$d\mathcal{E}_{\psi_{0}}(\delta{\psi})=\iint_{{\mathbb{S}}^{2}}\delta{U}\cdot
U_{0}\left[1-K^{\prime\prime}(F(\psi_{0}))F^{\prime}(\psi_{0})\right]\,dx+\int[A-\omega
K^{\prime\prime}(F(\psi_{0}))F^{\prime}(\psi_{0})]\cos\theta\,\mathbf{e}_{\varphi}\cdot\delta{U}\,{\rm
d}\sigma\,.$

It remains to choose $K^{\prime\prime}(F(x))F^{\prime}(x)=1$ and $A=\omega$ to
obtain that the first variation is zero.

Turning to the second variation,

$d^{2}\mathcal{E}_{{\psi_{0}}}(\delta{\psi})=\iint_{{\mathbb{S}}^{2}}\left[|\delta
U|^{2}+\frac{1}{F^{\prime}(\psi_{0})}(\delta\Omega)^{2}+2B(\mathbb{P}_{1}\delta\psi)^{2}\right]\,{\rm
d}\sigma\,,$

we expand $\psi$ in spherical harmonics: letting $\mathbb{P}_{k}$ denote the
projector on the $k$-th eigenspace of $\Delta$ associated to the eigenvalue
$-k(k+1)$, we write

$\delta\psi=\sum_{k\geq 1}\alpha_{k}\,\mathbb{P}_{k}\delta\psi\,,$

since we can subtract a constant from $\psi$ to ensure that
$\mathbb{P}_{0}\psi=0$. Then

$\displaystyle d^{2}\mathcal{E}_{0}(\delta\psi)$
$\displaystyle=\iint_{{\mathbb{S}}^{2}}\left[\sum_{k\geq
1}k(k+1)\,|\mathbb{P}_{k}\delta\psi|^{2}+\frac{1}{F^{\prime}(\psi_{0})}\sum_{k\geq
1}k^{2}(k+1)^{2}|\mathbb{P}_{k}\delta\psi|^{2}+2B\,|\mathbb{P}_{1}\delta\psi|^{2}\right]\,{\rm
d}\sigma$
$\displaystyle=\iint_{{\mathbb{S}}^{2}}\left[\left(2+\frac{4}{F^{\prime}(\psi_{0})}+2B\right)(\mathbb{P}_{1}\delta\psi)^{2}+\sum_{k\geq
2}\left(k(k+1)+\frac{k^{2}(k+1)^{2}}{F^{\prime}(\psi_{0})}\right)(\mathbb{P}_{k}\delta\psi)^{2}\right]\,{\rm
d}\sigma\,.$

The question of the coercivity of $d^{2}\mathcal{E}_{0}$ reduces to
determining eigenvalues of the Schrödinger operator
$-\frac{\Delta}{F^{\prime}(\psi)}+1$ on the orthogonal complement of
$\mathbb{E}_{1}$. We distinguish several cases, according to the the range of
$F^{\prime}(\psi_{0})$:

  * •

If $F^{\prime}>0$, the quadratic form is positive-definite if one chooses
$B=0$, for instance. This corresponds to Type II Arnold stable solutions (see
Subsection 1.3 for the definition). However, the condition $F^{\prime}>0$ is
only satisfied by constant solutions (see Theorem 4).

  * •

If $-6<F^{\prime}(\psi_{0})<0$, then the quadratic form is negative-definite:
indeed, $k(k+1)+\frac{k^{2}(k+1)^{2}}{F^{\prime}(\psi_{0})}>0$ for all $k\geq
2$, and the mode $k=1$ can be handled by choosing $B=-10$. This corresponds to
Type I Arnold stable solutions, with the difference that the first eigenvalue
of the Laplacian is replaced by the second eigenvalue – a modification related
to conservation laws which are specific to the sphere. Note that, due to
Theorem 4, the condition $-6<F^{\prime}(\psi_{0})<0$ forces solutions to be
zonal up to a rotation.

Combining the above considerations with standard arguments leads to the
following result.

######  Theorem 5.

For $0>F^{\prime}>-6$ any solution $\psi_{0}$ of (2.10) is stable in
$H^{2}({\mathbb{S}}^{2})$: if $\widehat{\psi}(t)$ is the solution of
(${\mathcal{E}}_{\omega}$) with initial data $\psi(0)=\tilde{\psi}_{0}$, then
$\|\widehat{\psi}(t)-\psi_{0}\|_{H^{2}}\lesssim\|\widehat{\psi}_{0}-\psi_{0}\|_{H^{2}}\,.$

The limiting case in the above theorem is given by $F^{\prime}=-6$, which
corresponds to Rossby-Haurwitz solutions in $\mathbb{E}_{1}+\mathbb{E}_{2}$.
They will be the focus of Section 5 below.

Theorem 5 applies to the explicit stationary solutions discussed in Section 3.
Indeed, since (3.2)-(3.3) yield

$\frac{1}{F^{\prime}(\psi)}=-\,\frac{1}{2(1-\varepsilon^{2})}\,\frac{[1-\varepsilon^{2}\cos^{2}(\theta)\sin^{2}(\varphi-\varphi_{0})]^{2}}{1+3\,\varepsilon^{2}\cos^{2}(\theta)\sin^{2}(\varphi-\varphi_{0})}\,,$

while (3.4)-(3.5) lead to

$\frac{1}{F^{\prime}(\psi_{0})}=-\,\frac{1}{2-\varepsilon^{2}+4\varepsilon\cos(\theta)\sin(\varphi-\varphi_{0})+\varepsilon^{2}\cos^{2}(\theta)\sin^{2}(\varphi-\varphi_{0})}\,,$

we see that for $\varepsilon>0$ small enough the solutions (3.3) and (3.5) are
stable. Regarding the physical relevance of these stable stationary solutions,
note that with the exception of the equatorial regions containing a broad
eastward zonal jet, vortices are generally found on Jupiter and Saturn at all
latitudes, preferentially in regions of westward zonal flow (see the data in
[57]).

###  4.3. Nonlinear stability of zonal flows

####  4.3.1. Arnold’s approach

We now investigate the nonlinear stability of smooth zonal flows
$\psi_{0}=f(\theta)$, with associated azimuthal velocity and vorticity given
by

$U_{0}=-f^{\prime}(\theta)\mathbf{e}_{\varphi}\,,\qquad\Omega_{0}=f^{\prime\prime}(\theta)-\tan\theta
f^{\prime}(\theta)=g(\theta)\,,$

respectively. The main result (Theorem 6) was proved in [52] in greater
generality (for rotating surfaces that are not necessarily spherical), but we
include it here for ease of reference, along with a more straightforward proof
in the case of the sphere.

######  Theorem 6.

If there exists $\epsilon,\,A\in{\mathbb{R}}$ such that

(4.15)
$\left|\frac{f^{\prime}(\theta)-A\cos\theta}{g^{\prime}(\theta)}\right|>\epsilon>0\qquad\mbox{on
$\quad\displaystyle\left(-\frac{\pi}{2},\frac{\pi}{2}\right)$}\,,$

stable in $H^{2}({\mathbb{S}}^{2})$: if $\widehat{\psi}(t)$ is the solution of
(${\mathcal{E}}_{\omega}$) with initial data $\psi(0)=\widehat{\psi}_{0}$,
then
$\|\widehat{\psi}(t)-\psi_{0}\|_{H^{2}}\lesssim\|\widehat{\psi}_{0}-\psi_{0}\|_{H^{2}}$.

Note that the condition (4.15) is satisfied if, for instance,

  * •

$|g^{\prime}(\theta)|>0$ on $(-\frac{\pi}{2},\frac{\pi}{2})$,

  * •

$|f^{\prime}(\theta)|\lesssim|\theta\pm\frac{\pi}{2}|$ for $\theta$ close to
$\mp\frac{\pi}{2}$,

  * •

$|g^{\prime}(\theta)|\gtrsim|\theta\pm\frac{\pi}{2}|$ for $\theta$ close to
$\mp\frac{\pi}{2}$,

hold simultaneously.

###### Proof.

We set

$\mathcal{E}(\psi)=\iint_{{\mathbb{S}}^{2}}\left[\frac{1}{2}|U|^{2}+K(\Omega)+A\sin\theta\Omega\right]\,{\rm
d}\sigma\,.$

Then, since
$\operatorname{grad}\Omega_{0}=g^{\prime}(\theta)\mathbf{e}_{\theta}$, we have

$\displaystyle d\mathcal{E}_{{\psi_{0}}}(\delta{\psi})$
$\displaystyle=\iint_{{\mathbb{S}}^{2}}\left[U_{0}\cdot\delta
U+K^{\prime}(\Omega_{0})\delta\Omega+A\sin\theta\delta\Omega\right]\,{\rm
d}\sigma$ $\displaystyle=\iint_{{\mathbb{S}}^{2}}\left[U_{0}\cdot\delta
U-K^{\prime\prime}(\Omega_{0})J\operatorname{grad}\Omega_{0}\cdot\delta
U+A\operatorname{grad}\sin\theta\cdot J\delta U\right]\,{\rm d}\sigma$
$\displaystyle=\iint_{{\mathbb{S}}^{2}}\left[-f^{\prime}(\theta)+K^{\prime\prime}(g(\theta))g^{\prime}(\theta)+A\cos\theta\right]\mathbf{e}_{\varphi}\cdot\delta
U_{0}\,{\rm d}\sigma\,.$

In order for this first variation to vanish, we choose $K$ such that
$-f^{\prime}(\theta)+K^{\prime\prime}(g(\theta))g^{\prime}(\theta)=-A\cos\theta$.
The second variation of $\mathcal{E}$ is then

$d^{2}\mathcal{E}_{{\psi_{0}}}(\delta{\psi})=\iint_{{\mathbb{S}}^{2}}\left[|\delta
U|^{2}+\frac{f^{\prime}(\theta)-A\cos\theta}{g^{\prime}(\theta)}(\delta\Omega)^{2}\right]\,{\rm
d}\sigma\,,$

which ensures stability. ∎

####  4.3.2. Application to the atmospheres of Uranus and Neptune

To illustrate the applicability of Theorem 6, let us now investigate the zonal
wind profiles of Uranus and Neptune, consisting of one broad retrograde
equatorial jet flanked by two prograde jets at higher latitudes (see Figure
2). The zonal flow is symmetric about the Equator for both planets, but there
are noticeable differences of the latitudinal flow profiles:

  * •

on Uranus the equatorial jet is located within the latitude band between 30∘N
and 30∘S, while on Neptune it extends over 50∘;

  * •

the prograde/retrograde (eastward/westward) zonal flows on Uranus, measured
relative to the planet’s rotation speed about its polar axis, peak at about
200 m/s, respectively at 80 m/s, the corresponding values for Neptune being
about 200 m/s and 400 m/s, respectively.

Recalling (2.9), if the latitudinal profile of the zonal flow with respect to
the rotation about the planet’s polar axis (with zonal velocity
$\theta\mapsto\omega\,\cos\theta$), depicted in Figure 2, is given by the
function

(4.16)
$U_{0}(\theta)=\alpha\,\cos^{5}\theta+\beta\,\cos^{3}\theta+\gamma\,\cos\theta\,,\qquad\theta\in\Big{(}-\frac{\pi}{2},\frac{\pi}{2}\Big{)}\,,$

for some real constants $\alpha\neq 0$, $\beta$ and $\gamma$, then

$U_{0}(\theta)=-f^{\prime}(\theta)-\omega\,\cos\theta\,,\qquad\theta\in\Big{[}-\frac{\pi}{2},\frac{\pi}{2}\Big{]}\,,$

with the notation of Theorem 6. We can now compute

(4.17)
$\frac{-f^{\prime}(\theta)+A\,\cos\theta}{g^{\prime}(\theta)}=\frac{\cos^{4}\theta+\frac{\beta}{\alpha}\,\cos^{2}\theta+\frac{\gamma-\omega+A}{\alpha}}{30\Big{(}\cos^{4}\theta+\frac{2(\beta-2\alpha)}{5\alpha}\,\cos^{2}\theta+\frac{\gamma-\omega-8\beta}{15\alpha}\Big{)}}\,,\qquad\theta\in\Big{(}-\frac{\pi}{2},\frac{\pi}{2}\Big{)}\,.$

so that the stability criterion provided in Theorem 6 applies if the quadratic
polynomials

(4.18)
$x^{2}+\frac{\beta}{\alpha}\,x+\frac{\gamma-\omega+A}{\alpha}\quad\text{and}\quad
x^{2}+\frac{2(\beta-2\alpha)}{5\alpha}\,x+\frac{\gamma-\omega-8\beta}{15\alpha}$

have the same roots in the interval $(0,1)$.

  * •

Since the unit of the non-dimensional zonal speed scaling for Uranus
corresponds to 150 m/s (see the first table in Section 7), the latitudinal
profile of the zonal flow on Uranus depicted in Figure 2 is well-approximated
by a function of the form (11) if we require $U_{0}$ to have the minimum
$U_{0}(0)=-\frac{8}{15}$ and the maximum $U_{0}(\frac{\pi}{3})=\frac{4}{3}$ on
$\big{[}0,\frac{\pi}{2}\big{]}$, these being the non-dimensional counterparts
of an eastward speed of 80 m/s and a westward speed of 200 m/s, respectively.
The condition that $\theta=\frac{\pi}{3}$ is a critical point and the above
two specific values of the function $U_{0}$ yield the linear system

$\begin{cases}&\alpha+\beta+\gamma=-\frac{8}{15}\,,\\\
&\tfrac{5}{16}\,\alpha+\tfrac{3}{4}\,\beta+\gamma=0\,,\\\
&\tfrac{1}{32}\,\alpha+\tfrac{1}{8}\,\beta+\tfrac{1}{2}\,\gamma=\tfrac{4}{3}\,,\end{cases}$

whose unique solution is

$\alpha=\frac{64}{45}\,,\quad\beta=-\frac{272}{45}\,,\quad\gamma=\tfrac{184}{45}\,.$

With $\omega=18$ the relevant value for Uranus (see the first table in Section
7), we obtain that the second quadratic polynomial from (4.18), expressed in
terms of $x=\cos^{2}\theta$, is

$x\mapsto x^{2}-\tfrac{5}{2}\,x+\tfrac{155}{96}\,,$

with no real roots. Choosing $A\in{\mathbb{R}}$ so that the first quadratic
polynomial in (4.18) has no roots, we can apply Theorem 6 and we conclude that
the zonal flow pattern of Uranus is stable.

  * •

For Neptune the non-dimensional unit for the zonal speed corresponds to 200
m/s (see the first table in Section 7), the latitudinal profile of the zonal
flow depicted in Figure 2 is well-approximated by a function of the form (11)
if we require $U_{0}$ to have the minimum $U_{0}(0)=-2$ and the maximum
$U_{0}(\frac{5\pi}{12})=1$ on $\big{[}0,\frac{\pi}{2}\big{]}$. The condition
that $\theta=\frac{5\pi}{12}$ with $\cos\theta\approx\tfrac{1}{4}$ is a
critical point and the above two specific values of the function $U_{0}$ yield
the linear system

$\begin{cases}&\alpha+\beta+\gamma=-2\,,\\\
&\tfrac{5}{4^{4}}\,\alpha+\tfrac{3}{4^{2}}\,\beta+\gamma=0\,,\\\
&\tfrac{1}{4^{5}}\,\alpha+\tfrac{1}{4^{3}}\,\beta+\tfrac{1}{4}\,\gamma=1\,,\end{cases}$

whose unique solution is $\alpha=\frac{2048}{75}$, $\beta=-\frac{2656}{75}$,
$\gamma=\tfrac{458}{75}$. With $\omega=13$ the relevant value for Neptune (see
the first table in Section 7), we obtain that the second quadratic polynomial
in $x=\cos^{2}\theta$ from (4.18) is

$x\mapsto x^{2}-\tfrac{211}{160}\,x+\tfrac{20731}{30720}\,,$

with no real roots as $\tfrac{211}{160}\approx 1.31$ and
$\tfrac{20731}{30720}\approx 0.674$. Consequently, choosing $A\in{\mathbb{R}}$
so that the first quadratic polynomial in (4.18) has no roots, Theorem 6
yields that the zonal flow pattern of Neptune is also stable.

Unfortunately this approach is not applicable to the likely stability of the
zonal jets on Jupiter and Saturn to breaking up into meanders and vortex-like
eddies, but both zonal jet patterns (see [45] for their detailed profiles) are
not far off from entering the framework of Theorem 6. In contrast to this, the
profiles of terrestrial stratospheric jets derived from observational data
(see [41]) are well beyond condition (4.15), as is to be expected since the
Earth’s polar jet stream is known to be unstable.

##  5\. Stability results for degree 2 Rossby-Haurwitz waves

Due to the considerable physical relevance of the largest-scale Rossby-
Haurwitz waves (2.12) (that is, those with low degree), their stability
properties are of great interest.

######  Theorem 7.

(i) (Stability of zonal Rossby-Haurwitz flows of degree $n\leq 2$) Zonal
solutions $\psi_{0}$ to (${\mathcal{E}}_{\omega}$) of the form

(5.1) $\psi_{0}(\theta)=\alpha\,\sin\theta+\beta
Y_{2}^{0}(\theta)\,,\qquad\alpha\in{\mathbb{R}}\,,\quad\beta\in{\mathbb{R}}\setminus\\{0\\}\,,$

are stable in $H^{2}({\mathbb{S}}^{2})$ under perturbations with bounded
vorticity.

(ii) (Instability of non-zonal Rossby-Haurwitz waves) Non-zonal Rossby-
Haurwitz waves $\psi_{0}$ of the form (2.12) are unstable. To be more
specific, there exists $\epsilon>0$ and a sequence
$\widehat{\psi}^{n}_{0}\to\psi_{0}$ in $H^{2}({\mathbb{S}}^{2})$, so that for
the solutions $\widehat{\psi}^{n}(t)$ of (${\mathcal{E}}_{\omega}$) with
initial data $\widehat{\psi}^{n}(0)=\widehat{\psi}^{n}_{0}$ we have

$\sup_{t>0}\|\psi^{n}(t)-\psi_{0}(t)\|_{L^{2}({\mathbb{S}}^{2})}>\epsilon>0.$

(iii) (Stability of $\mathbb{E}_{1}+\mathbb{E}_{2}$) The instability of a
Rossby-Haurwitz wave of degree 2 can only occur by energy transfers between
spherical harmonics of degree $2$. More precisely, if we consider a
perturbation $\widehat{\psi}_{0}\in H^{2}({\mathbb{S}}^{2})$ of

(5.2)
$\psi_{0}(\varphi,\theta)=-\frac{\omega}{2}\,\sin\theta+\beta_{0}Y(\varphi,\theta)\,,$

where $Y\in\mathbb{E}_{2}$ with $\|Y\|_{L^{2}({\mathbb{S}}^{2})}=1$ and
$\beta_{0}\in{\mathbb{R}}\setminus\\{0\\}$, the solution $\widehat{\psi}(t)$
of (${\mathcal{E}}_{\omega}$) with initial data
$\widehat{\psi}(0)=\widehat{\psi}_{0}$ can be written as

$\widehat{\psi}(t)=-\alpha\,\sin\theta+c_{1}^{-1}\,{\rm e}^{{\rm i}\omega
t}\,Y_{1}^{-1}+c_{1}^{1}\,{\rm e}^{-{\rm i}\omega
t}\,Y_{1}^{1}+\beta(t)\,Y(t)+\tilde{\psi}(t)\,,\qquad t\geq 0\,,$

where

$\left\\{\begin{array}[]{l}c_{1}^{\pm
1}\in{\mathbb{C}}\quad\text{and}\quad\alpha,\,\beta(t)\in\mathbb{R}\,,\\\
Y(t)\in\mathbb{E}_{2}\quad\text{with}\quad\|Y(t)\|_{L^{2}({\mathbb{S}}^{2})}=1\,,\\\
\mathbb{P}_{1}\tilde{\psi}(t)=\mathbb{P}_{2}\tilde{\psi}(t)=0\quad\text{for
all}\quad t\geq 0\,,\\\
|\alpha-\tfrac{\omega}{2}|+|c_{1}^{-1}|+|c_{1}^{1}|+|\beta(t)-\beta_{0}|+\|\tilde{\psi}(t)\|_{H^{2}({\mathbb{S}}^{2})}\lesssim\|\psi_{0}-\widehat{\psi}_{0}\|_{H^{2}({\mathbb{S}}^{2})}\,.\end{array}\right.$

######  Remark 1.

It seems natural to conjecture that all Rossby-Haurwitz waves of degree 2 are
orbitally stable. It should be possible to adapt the proof of assertion $(i)$
to prove this conjecture, though this would entail considerable technical
complications.

###### Proof.

$(i)$ We begin with a reduction: by using the symmetries of the problem (see
Subsection 2.3), more specifically the scaling and change-of-frame symmetries,
it suffices to prove $(i)$ in the case $\omega=0$, $\beta=1$.

Consider a smooth perturbation

(5.3)
$\psi(\varphi,\theta,t)=\sum_{l=1}^{\infty}\Big{\\{}\sum_{m=-l}^{l}c_{l}^{m}(t)\,Y_{l}^{m}(\varphi,\theta)\Big{\\}}$

of the zonal flow (5.1), expressed in terms of the spherical harmonics
$Y_{l}^{m}$ by means of the time-dependent coefficients
$c_{l}^{m}(t)\in{\mathbb{C}}$. Since $\psi$ is real-valued, (A.1) yields

(5.4)
$c_{l}^{-m}(t)=\int_{{\mathbb{S}}^{2}}\psi\,\overline{Y_{l}^{-m}}\,d\sigma=(-1)^{m}\int_{{\mathbb{S}}^{2}}\psi\,Y_{l}^{m}\,d\sigma=(-1)^{m}\overline{\int_{{\mathbb{S}}^{2}}\psi\,\overline{Y_{l}^{m}}\,d\sigma}=(-1)^{m}\,\overline{c_{l}^{m}(t)}\,,\quad|m|\leq
l\,.$

Furthermore, we know that

(5.5) $c_{1}^{0}(t)=c_{1}^{0}(0)\,,\quad|c_{1}^{\pm 1}(t)=|c_{1}^{\pm
1}(0)|\,,\qquad t\geq 0\,.$

Assuming that initially (at time $t=0$) the solution (5.3) of
(${\mathcal{E}}_{\omega}$) is $\varepsilon$-close to the Rossby flow (5.1) in
the $H^{2}$-norm, with $\varepsilon>0$ small, we have

(5.6) $\displaystyle\|\Delta[\psi_{0}-\psi(\cdot,0)]\|_{L^{2}}^{2}=$
$\displaystyle 4\,|\alpha-
c_{1}^{0}(0)|^{2}+4\,|c_{1}^{-1}(0)|^{2}+4\,|c_{1}^{1}(0)|^{2}+36\,|\beta-
c_{2}^{0}(0)|^{2}$ $\displaystyle\quad+36\sum_{0<|m|\leq
2}^{l}|c_{2}^{m}(0)|^{2}+\sum_{l=3}^{\infty}l^{2}(l+1)^{2}\Big{\\{}\sum_{m=-l}^{l}|c_{l}^{m}(0)|^{2}\Big{\\}}<\varepsilon^{2}\,.$

On the other hand, the conservation of energy for the solution (5.3) to
(${\mathcal{E}}_{\omega}$) reads

$\sum_{l=1}^{\infty}l(l+1)\Big{\\{}\sum_{m=-l}^{l}|c_{l}^{m}(t)|^{2}\Big{\\}}=\sum_{l=1}^{\infty}l(l+1)\Big{\\{}\sum_{m=-l}^{l}|c_{l}^{m}(0)|^{2}\Big{\\}}\,,\qquad
t\geq 0\,,$

which, using (5.5), can be re-written in the form

(5.7)
$\sum_{l=2}^{\infty}l(l+1)\Big{\\{}\sum_{m=-l}^{l}|c_{l}^{m}(t)|^{2}\Big{\\}}=\sum_{l=2}^{\infty}l(l+1)\Big{\\{}\sum_{m=-l}^{l}|c_{l}^{m}(0)|^{2}\Big{\\}}\,,\qquad
t\geq 0\,.$

The time-invariance of the integral
$\int_{{\mathbb{S}}^{2}}|\Delta\psi|^{2}\,{\rm d}\sigma$ gives the equality

$\displaystyle
4\sum_{m=-1}^{1}|c_{1}^{m}(t)|^{2}+\sum_{l=2}^{\infty}l^{2}(l+1)^{2}\Big{\\{}\sum_{m=-l}^{l}|c_{l}^{m}(t)|^{2}\Big{\\}}=4\sum_{m=-1}^{1}|c_{1}^{m}(0)|^{2}+\sum_{l=2}^{\infty}l^{2}(l+1)^{2}\Big{\\{}\sum_{m=-l}^{l}|c_{l}^{m}(0)|^{2}\Big{\\}}\,,\qquad
t\geq 0\,.$

Using (5.5), we infer that

(5.8)
$\sum_{l=2}^{\infty}l^{2}(l+1)^{2}\Big{\\{}\sum_{m=-l}^{l}|c_{l}^{m}(t)|^{2}\Big{\\}}=\sum_{l=2}^{\infty}l^{2}(l+1)^{2}\Big{\\{}\sum_{m=-l}^{l}|c_{l}^{m}(0)|^{2}\Big{\\}}\,,\qquad
t\geq 0\,.$

The identities (5.7) and (5.8) ensure

$\sum_{l=2}^{\infty}[l^{2}(l+1)^{2}-6l(l+1)]\Big{\\{}\sum_{m=-l}^{l}|c_{l}^{m}(t)|^{2}\Big{\\}}=\sum_{l=2}^{\infty}[l^{2}(l+1)^{2}-6l(l+1)]\Big{\\{}\sum_{m=-l}^{l}|c_{l}^{m}(0)|^{2}\Big{\\}}\,,\qquad
t\geq 0\,,$

which, since the $l=2$ coefficient vanishes, can be written

(5.9)
$\sum_{l=3}^{\infty}[l^{2}(l+1)^{2}-6l(l+1)]\Big{\\{}\sum_{m=-l}^{l}|c_{l}^{m}(t)|^{2}\Big{\\}}<\varepsilon^{2}\,,\qquad
t\geq 0\,.$

Recalling (5.5), we conclude that the instability of the zonal Rossby flow
(5.1) can only be caused by a substantial energy transfer between the
spherical harmonic components of mode $l=2$. To rule this out, we rely on the
time-invariance of the integrals

(5.10)
$I_{k}(\psi(\cdot,t))=\int_{{\mathbb{S}}^{2}}\Big{(}\Delta\psi(\cdot,t)\Big{)}^{k}\,{\rm
d}\sigma\,,\qquad k\in\\{2,3,5\\}\,,$

which ensures, using the Cauchy-Schwarz inequality and the boundedness of the
vorticity,

(5.11)
$|I_{k}(\psi(\cdot,t))-I_{k}(\psi_{0})|=|I_{k}(\psi(\cdot,0))-I_{k}(\psi_{0})|\\\
\lesssim\|\Delta\psi(\cdot,0)-\Delta\psi_{0}\|_{L^{2}}<\varepsilon\,.$

We now take advantage of (5.9) and of the specific structure of the spherical
harmonics to elucidate the leading order of the integrals in (5.10) as
$\varepsilon\to 0$. For this, note that integration by parts yields the
recurrence formula

$\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos^{2k+1}\theta\,{\rm
d}\theta=\frac{2k}{2k+1}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos^{2k-1}\theta\,{\rm
d}\theta\,,\qquad k\geq 1\,,$

which, since $\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos\theta\,{\rm
d}\theta=2$, yields the value of the Wallis integrals

$\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos^{3}\theta\,{\rm
d}\theta=\tfrac{4}{3}\,,\quad\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos^{5}\theta\,{\rm
d}\theta=\tfrac{16}{15}\,,\quad\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos^{7}\theta\,{\rm
d}\theta=\tfrac{32}{35}\,,\quad\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos^{9}\theta\,{\rm
d}\theta=\tfrac{256}{315}\,,\quad\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos^{11}\theta\,{\rm
d}\theta=\tfrac{512}{693}\,.$

Taking into account the explicit formulas for the spherical harmonics (see the
Appendix), we can now compute

$\displaystyle\int_{{\mathbb{S}}^{2}}(Y_{2}^{0})^{2}\,{\rm
d}\sigma=1\,,\quad\int_{{\mathbb{S}}^{2}}Y_{2}^{-1}Y_{2}^{1}\,{\rm
d}\sigma=-1\,,\quad\int_{{\mathbb{S}}^{2}}Y_{2}^{-2}Y_{2}^{2}\,{\rm
d}\sigma=1\,,\quad\int_{{\mathbb{S}}^{2}}(Y_{2}^{0})^{3}\,{\rm
d}\sigma=\tfrac{\sqrt{5}}{7\sqrt{\pi}}\,,$
$\displaystyle\int_{{\mathbb{S}}^{2}}Y_{2}^{0}Y_{2}^{-1}Y_{2}^{1}\,{\rm
d}\sigma=-\tfrac{\sqrt{5}}{14\sqrt{\pi}}\,,\quad\int_{{\mathbb{S}}^{2}}Y_{2}^{0}Y_{2}^{-2}Y_{2}^{2}\,{\rm
d}\sigma=-\tfrac{\sqrt{5}}{7\sqrt{\pi}}\,,\quad\int_{{\mathbb{S}}^{2}}Y_{2}^{2}(Y_{2}^{-1})^{2}\,{\rm
d}\sigma=\tfrac{\sqrt{15}}{7\sqrt{2\pi}}\,,\quad\int_{{\mathbb{S}}^{2}}Y_{2}^{-2}(Y_{2}^{1})^{2}\,{\rm
d}\sigma=\tfrac{\sqrt{15}}{7\sqrt{2\pi}}\,,$
$\displaystyle\int_{{\mathbb{S}}^{2}}(Y_{2}^{0})^{4}\,{\rm
d}\sigma=\tfrac{15}{28\pi}\,,\quad\int_{{\mathbb{S}}^{2}}(Y_{2}^{0})^{2}Y_{2}^{-2}Y_{2}^{2}\,{\rm
d}\sigma=\tfrac{5}{28\pi}\,,\quad\int_{{\mathbb{S}}^{2}}(Y_{2}^{0})^{2}Y_{2}^{-1}Y_{2}^{1}\,{\rm
d}\sigma=-\tfrac{5}{28\pi}\,,\quad\int_{{\mathbb{S}}^{2}}Y_{2}^{0}Y_{2}^{2}(Y_{2}^{-1})^{2}\,{\rm
d}\sigma=0\,,$
$\displaystyle\int_{{\mathbb{S}}^{2}}Y_{2}^{0}Y_{2}^{-2}(Y_{2}^{1})^{2}\,{\rm
d}\sigma=0\,,\quad\int_{{\mathbb{S}}^{2}}(Y_{2}^{2})^{2}(Y_{2}^{-2})^{2}\,{\rm
d}\sigma=\tfrac{5}{14\pi}\,,\quad\int_{{\mathbb{S}}^{2}}(Y_{2}^{1})^{2}(Y_{2}^{-1})^{2}\,{\rm
d}\sigma=\tfrac{5}{14\pi}\,,\quad\int_{{\mathbb{S}}^{2}}Y_{2}^{-2}Y_{2}^{2}Y_{2}^{-1}Y_{2}^{1}\,{\rm
d}\sigma=-\tfrac{5}{28\pi}\,,$
$\displaystyle\int_{{\mathbb{S}}^{2}}(Y_{2}^{0})^{5}\,{\rm
d}\sigma=\tfrac{5^{2}\cdot
199\sqrt{5}}{154\pi\sqrt{\pi}}\,,\quad\int_{{\mathbb{S}}^{2}}(Y_{2}^{0})^{3}Y_{2}^{-2}Y_{2}^{2}\,{\rm
d}\sigma=-\tfrac{5\sqrt{5}}{154\pi\sqrt{\pi}}\,,\quad\int_{{\mathbb{S}}^{2}}(Y_{2}^{0})^{3}Y_{2}^{-1}Y_{2}^{1}\,{\rm
d}\sigma=-\tfrac{5^{2}\sqrt{5}}{4\cdot 154\pi\sqrt{\pi}}\,.$

We now use (5.10) and the multinomial formula

$(x_{1}+\dots+x_{n})^{k}=\sum_{k_{1}+\dots+k_{n}=k}\frac{k!}{k_{1}!\cdots
k_{n}!}\prod_{j=1}^{n}x_{j}^{k_{j}}\,,$

under the assumption (5.6), which ensures (5.9). Noticing that the
$\varphi$-dependence of the integrand shows that the integral
$\int_{{\mathbb{S}}^{2}}\prod_{j=1}^{n}(Y_{2}^{m_{j}})^{k_{j}}\,{\rm d}\sigma$
vanishes if $\sum_{j=1}^{n}m_{j}k_{j}\neq 0$, from (5.4) we obtain

(5.12) $\displaystyle I_{2}((\psi(\cdot,t))$
$\displaystyle=4|c_{1}^{0}(t)|^{2}+8\,|c_{1}^{1}(t)|^{2}+\sum_{l=2}^{\infty}l^{2}(l+1)^{2}\Big{\\{}|c_{l}^{0}(t)|^{2}+2\sum_{m=1}^{l}|c_{l}^{m}(t)|^{2}\Big{\\}}\,,$
(5.13) $\displaystyle I_{2}(\psi_{0})$ $\displaystyle=4\alpha^{2}+36\,.$

Since (5.5) and (5.6) ensure

(5.14)
$\Big{|}|c_{1}^{0}(t)|^{2}-\alpha^{2}\Big{|}=|c_{1}^{0}(t)-\alpha|\cdot\Big{|}c_{1}^{0}(t)+\alpha\Big{|}=|c_{1}^{0}(0)-\alpha|\cdot\Big{|}c_{1}^{0}(0)+\alpha\Big{|}\lesssim\varepsilon\,,$

we see that (5.9) and (5.12)-(5.13) yield

$I_{2}((\psi(\cdot,t))-I_{2}(\psi_{0})=36\,\Big{\\{}[c_{2}^{0}(t)]^{2}+2|c_{2}^{2}(t)|^{2}+2|c_{2}^{1}(t)|^{2}-1\Big{\\}}+{\rm
O}(\varepsilon)\,.$

From (5.11) we therefore get

(5.15) $[c_{2}^{0}(t)]^{2}+2\,|c_{2}^{2}(t)|^{2}+2\,|c_{2}^{1}(t)|^{2}=1+{\rm
O}(\varepsilon)\,,$

We now show that

(5.16)
$\displaystyle\tfrac{7}{15}\alpha^{2}c_{2}^{0}(t)+[c_{2}^{0}(t)]^{3}+6\,c_{2}^{0}(t)|c_{2}^{2}(t)|^{2}-3\,c_{2}^{0}(t)|c_{2}^{1}(t)|^{2}+3\sqrt{6}\,{\mathfrak{R}e}\Big{[}c_{2}^{2}(t)\,\Big{(}c_{2}^{1}(t)\Big{)}^{2}\Big{]}=\tfrac{7}{15}\alpha^{2}+1+{\rm
O}(\varepsilon)\,,$

where we denote by ${\mathfrak{R}e}(z)$ the real part of the complex number
$z$. For this, note that $c_{1}^{\pm 1}(t)$ and $c_{l}^{m}(t)$ for $l\geq 3$
are ${\rm O}(\varepsilon)$, in view of (5.5), (5.6) and (5.9). Since
$\int_{{\mathbb{S}}^{2}}\prod_{j=1}^{3}(Y_{2}^{m_{j}})^{k_{j}}\,{\rm d}\sigma$
vanishes if $\sum_{j=1}^{3}m_{j}k_{j}\neq 0$, and

$\int_{{\mathbb{S}}^{2}}(Y_{1}^{0})^{3}\,{\rm
d}\sigma=\int_{{\mathbb{S}}^{2}}Y_{1}^{0}(Y_{2}^{0})^{2}\,{\rm
d}\sigma=\int_{{\mathbb{S}}^{2}}Y_{1}^{0}Y_{2}^{-2}Y_{2}^{2}\,{\rm
d}\sigma=\int_{{\mathbb{S}}^{2}}Y_{1}^{0}Y_{2}^{-1}Y_{2}^{1}\,{\rm d}\sigma=0$

because in each case we integrate an odd function of $\theta$ over
$\big{[}-\frac{\pi}{2},\,\frac{\pi}{2}\big{]}$, while

$\int_{{\mathbb{S}}^{2}}(Y_{1}^{0})^{2}Y_{2}^{0}\,{\rm
d}\sigma=\tfrac{3\sqrt{5}}{8\sqrt{\pi}}\,\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}(3\sin^{2}\theta-1)\sin^{2}\theta\cos\theta\,{\rm
d}\theta=\tfrac{1}{\sqrt{5\pi}}\,,$

we obtain

$\displaystyle I_{3}(\psi(\cdot,t))$
$\displaystyle=-\tfrac{72}{\sqrt{5\pi}}|c_{1}^{0}(t)|^{2}c_{2}^{0}(t)$
$\displaystyle\qquad-\tfrac{216\sqrt{5}}{7\sqrt{\pi}}\,\Big{\\{}[c_{2}^{0}(t)]^{3}+6\,c_{2}^{0}(t)|c_{2}^{2}(t)|^{2}-3\,c_{2}^{0}(t)|c_{2}^{1}(t)|^{2}+3\sqrt{6}\,{\mathfrak{R}}e\Big{[}c_{2}^{2}(t)\,\Big{(}c_{2}^{1}(t)\Big{)}^{2}\Big{]}\Big{\\}}+{\rm
O}(\varepsilon)\,.$

Due to (5.5) and (5.11) for $k=3$, the relation (5.16) now emerges by
subtracting from the above

$I_{3}(\psi_{0})=-\tfrac{72}{\sqrt{5\pi}}\alpha^{2}-\tfrac{216\sqrt{5}}{7\sqrt{\pi}}$

and taking (5.14) into account.

We now investigate the leading order of $I_{4}((\psi(\cdot,t))$. Again, since
$c_{1}^{\pm 1}(t)$ and $c_{l}^{m}(t)$ for $l\geq 3$ are ${\rm
O}(\varepsilon)$, for this we need only to keep track of the integrals
$\int_{{\mathbb{S}}^{2}}\prod_{j=1}^{3}(Y_{1}^{0})^{k}(Y_{2}^{m_{j}})^{k_{j}}\,{\rm
d}\sigma$ with $\sum_{j=1}^{3}m_{j}k_{j}=0$ and $k+\sum_{j=1}^{3}k_{j}=4$. We
compute

$\displaystyle\int_{{\mathbb{S}}^{2}}(Y_{1}^{0})^{4}\,{\rm
d}\sigma=\tfrac{9}{20\pi}\,,\quad\int_{{\mathbb{S}}^{2}}Y_{1}^{0}(Y_{2}^{0})^{3}\,{\rm
d}\sigma=0\,,\quad\int_{{\mathbb{S}}^{2}}(Y_{1}^{0})^{2}(Y_{2}^{0})^{2}\,{\rm
d}\sigma=\tfrac{11}{28\pi}\,,$
$\displaystyle\int_{{\mathbb{S}}^{2}}(Y_{1}^{0})^{2}Y_{2}^{2}Y_{2}^{-2}\,{\rm
d}\sigma=\tfrac{3}{14\pi}\,,\qquad\int_{{\mathbb{S}}^{2}}(Y_{1}^{0})^{2}Y_{2}^{1}Y_{2}^{-1}\,{\rm
d}\sigma=\tfrac{9}{28\pi}\,,$

and infer that

$\displaystyle I_{4}((\psi(\cdot,t))$
$\displaystyle=\tfrac{36}{5\pi}\,|c_{1}^{0}(t)|^{4}+\tfrac{6^{3}\cdot
11}{7\pi}|c_{1}^{0}(t)|^{2}[c_{2}^{0}(t)]^{2}+\tfrac{6^{4}\cdot
4}{7\pi}\,|c_{1}^{0}(t)|^{2}|c_{2}^{2}(t)|^{2}-\tfrac{6^{4}\cdot
3}{7\pi}|c_{1}^{0}(t)|^{2}|c_{2}^{1}(t)|^{2}\,,$
$\displaystyle\qquad+\tfrac{6^{4}\cdot
15}{28\pi}\Big{\\{}[c_{2}^{0}(t)]^{2}+2\Big{(}|c_{2}^{2}(t)|^{2}+|c_{2}^{1}(t)|^{2}\Big{)}\Big{\\}}^{2}+{\rm
O}(\varepsilon)\,.$

Subtracting from this the relation

$I_{4}(\psi_{0})=\tfrac{36}{5\pi}\alpha^{4}+\tfrac{6^{3}\cdot
11}{7\pi}\alpha^{2}+\tfrac{6^{4}\cdot 15}{28\pi}\,,$

invoking (5.11) with $k=4$ and (5.5), we get

(5.17)
$22\alpha^{2}[c_{2}^{0}(t)]^{2}+48\alpha^{2}|c_{2}^{2}(t)|^{2}-36\alpha^{2}|c_{2}^{1}(t)|^{2}+45\Big{\\{}[c_{2}^{0}(t)]^{2}+2\Big{(}|c_{2}^{2}(t)|^{2}+|c_{2}^{1}(t)|^{2}\Big{)}\Big{\\}}^{2}=22\alpha^{2}+45+{\rm
O}(\varepsilon)\,,$

using (5.14) and the fact that (5.6) ensures

$\displaystyle\Big{|}|c_{1}^{0}(t)|^{4}-\alpha^{4}\Big{|}\lesssim\varepsilon\,.$

We can now prove the stability of the zonal solution (5.1) for $\alpha\neq 0$.
In this case, from (5.15) and (5.17) we get

(5.18)
$[c_{2}^{0}(t)]^{2}+\tfrac{24}{11}\,|c_{2}^{2}(t)|^{2}-\tfrac{18}{11}\,|c_{2}^{1}(t)|^{2}=1+{\rm
O}(\varepsilon)\,,$

which, subtracted from (5.15), yields

(5.19) $|c_{2}^{2}(t)|^{2}-15\,|c_{2}^{1}(t)|^{2}={\rm O}(\varepsilon)\,.$

In view of the continuous dependence on data guaranteed by the well-posedness
of (${\mathcal{E}}_{\omega}$), invoking (5.5) and the fact that (5.9) is
guaranteed by (5.6), we see that the zonal Rossby flow $\psi_{0}$ given by
(5.1) is unstable only if there exists some $\delta_{0}>0$ and a sequence of
initial data $\\{\psi_{N}(\cdot,0)\\}_{N\geq 1}$ converging to $\Psi$ in
$H^{2}$ and such that for every $N_{0}\geq 1$ large enough there exists a time
$t_{N_{0}}>0$ with

(5.20)
$|c_{N_{0},2}^{2}(t_{N_{0}})|^{2}+|c_{N_{0},2}^{1}(t_{N_{0}})|^{2}=2\delta_{0}^{2}\,,$

where $\\{c_{N_{0},l}^{m}(t)\\}_{l\geq 1,\,|m|\leq l}$ are the coefficients of
the expansion of the solution $\psi_{N_{0}}(\cdot,t)$ to
(${\mathcal{E}}_{\omega}$) with initial data $\psi_{N_{0}}(\cdot,0)$. The
validity of (5.20) for some $\delta_{0}>0$ ensures a similar relation for any
$\delta\in(0,\delta_{0})$, at some time $T>0$. We can thus analyse the limit
$\delta\to 0$: from (5.15), (5.19) and (5.20) we get

(5.21) $\begin{cases}&[c_{N,2}^{0}(t_{N})]^{2}=1-4\delta^{2}+{\rm
O}(\varepsilon)\,,\\\ &|c_{N,2}^{2}(t_{N})|^{2}=\tfrac{15}{8}\,\delta^{2}+{\rm
O}(\varepsilon)\,,\\\ &|c_{N,2}^{1}(t_{N})|^{2}=\tfrac{1}{8}\,\delta^{2}+{\rm
O}(\varepsilon)\,.\end{cases}$

Thus

$c_{N,2}^{0}(t_{N})=(1-2\delta^{2}-\tfrac{1}{2}\,\delta^{4})+{\rm
O}(\delta^{6})+{\rm O}(\varepsilon)\,,$

which, together with the last two relations in (5.21) makes the matching at
${\rm O}(\delta^{3})$ in (5.16) impossible. This contradiction proves the
stability of (5.1) if $\alpha\neq 0$.

It remains to deal with the case $\alpha=0$, a setting in which (5.16)
simplifies to

(5.22)
$[c_{2}^{0}(t)]^{3}+6\,c_{2}^{0}(t)|c_{2}^{2}(t)|^{2}-3\,c_{2}^{0}(t)|c_{2}^{1}(t)|^{2}+3\sqrt{6}\,{\mathfrak{R}}e\Big{[}c_{2}^{2}(t)\,\Big{(}c_{2}^{1}(t)\Big{)}^{2}\Big{]}=1+{\rm
O}(\varepsilon)\,.$

but (5.17) does not provide additional information with respect to (5.15), so
that we can not rely on (5.19). To compensate for the ineffectiveness of
(5.17) we have to take advantage of (5.11) with $k=5$. Assuming instability,
we find for every $\delta>0$ small enough an initial data $\varepsilon$-close
to (5.1) such that for the corresponding solution $\psi_{N}(\cdot,t)$ at some
time $T_{N}>0$ we have

(5.23) $|c_{N,2}^{2}(t_{N})|^{2}+|c_{N,2}^{1}(t_{N})|^{2}=2\delta^{2}\,,$

From (5.15) we then get

$[c_{N,2}^{0}(t_{N})]^{2}=(1-4\delta^{2})+{\rm O}(\varepsilon)\,,$

so that

(5.24) $c_{2}^{0}(t_{N})=(1-2\delta^{2})+{\rm O}(\delta^{4})+{\rm
O}(\sqrt{\varepsilon})\,,$

and by writing

$6\,c_{2}^{0}(t)|c_{2}^{2}(t)|^{2}-3\,c_{2}^{0}(t)|c_{2}^{1}(t)|^{2}=3\,c_{2}^{0}(t)\,\Big{\\{}|c_{2}^{2}(t)|^{2}+|c_{2}^{1}(t)|^{2}\Big{\\}}+3\,c_{2}^{0}(t)\,\Big{\\{}|c_{2}^{2}(t)|^{2}-2|c_{2}^{1}(t)|^{2}\Big{\\}}$

the relations (5.22)-(5.23) yield

(5.25) $|c_{N,2}^{2}(t_{N})|^{2}-2\,|c_{N,2}^{1}(t_{N})|^{2}={\rm
O}(\delta^{3})+{\rm O}(\sqrt{\varepsilon})\,.$

We now notice that $I_{5}(\psi_{N}(\cdot,t_{N})$ is asymptotically an additive
${\rm O}(\varepsilon)$-correction of linear combinations of integrals of the
type $\int_{{\mathbb{S}}^{2}}\prod_{j=1}^{5}(Y_{2}^{m_{j}})^{k_{j}}\,{\rm
d}\sigma$ with $\sum_{j=1}^{5}m_{j}k_{j}=0$, so that by computing the eleven
integrals

$\displaystyle\int_{{\mathbb{S}}^{2}}(Y_{2}^{0})^{5}\,{\rm
d}\sigma\,,\quad\int_{{\mathbb{S}}^{2}}(Y_{2}^{0})^{3}Y_{2}^{-2}Y_{2}^{2}\,{\rm
d}\sigma\,,\quad\int_{{\mathbb{S}}^{2}}(Y_{2}^{0})^{3}Y_{2}^{-1}Y_{2}^{1}\,{\rm
d}\sigma\,,\quad\int_{{\mathbb{S}}^{2}}(Y_{2}^{0})^{2}Y_{2}^{2}(Y_{2}^{-1})^{2}\,{\rm
d}\sigma\,,\quad\int_{{\mathbb{S}}^{2}}(Y_{2}^{0})^{2}Y_{2}^{-2}(Y_{2}^{1})^{2}\,{\rm
d}\sigma\,,$
$\displaystyle\int_{{\mathbb{S}}^{2}}Y_{2}^{0}(Y_{2}^{-2})^{2}(Y_{2}^{2})^{2}\,{\rm
d}\sigma\,,\quad\int_{{\mathbb{S}}^{2}}Y_{2}^{0}(Y_{2}^{-1})^{2}\,{\rm
d}\sigma\,,\quad\int_{{\mathbb{S}}^{2}}Y_{2}^{0}(Y_{2}^{-1})^{2}(Y_{2}^{1})^{2}\,{\rm
d}\sigma\,,\quad\int_{{\mathbb{S}}^{2}}Y_{2}^{0}Y_{2}^{-2}Y_{2}^{2}Y_{2}^{-1}Y_{2}^{1}\,{\rm
d}\sigma\,,$
$\displaystyle\int_{{\mathbb{S}}^{2}}(Y_{2}^{2})^{2}Y_{2}^{-2}(Y_{2}^{-1})^{2}\,{\rm
d}\sigma\,,\quad\int_{{\mathbb{S}}^{2}}(Y_{2}^{-2})^{3}Y_{2}^{2}(Y_{2}^{1})^{2}\,{\rm
d}\sigma\,,$

we can determine the leading order with respect to $\varepsilon$. If we also
bring $\delta\gg\varepsilon$ into play, then (5.23) and (5.24) ensure that
$I_{5}(\psi_{N}(\cdot,t_{N})$ consists of a linear combination of the first
three integrals (already computed) and an additive ${\rm
O}(\delta^{3},\varepsilon)$-correction:

(5.26) $I_{5}(\psi_{N}(\cdot,t_{N})=-\tfrac{6^{5}\cdot 5^{2}\cdot
199\sqrt{5}}{154\pi\sqrt{\pi}}\,\Big{\\{}[c_{2}^{0}(t_{N})]^{5}+\frac{1}{199}\,[c_{2}^{0}(t_{N})]^{3}[c_{N,2}^{1}(t_{N})]^{2}-\frac{4}{199}\,[c_{2}^{0}(t_{N})]^{3}|c_{N,2}^{2}(t_{N})|^{2}\Big{\\}}+{\rm
O}(\delta^{3})+{\rm O}(\varepsilon)\,,$

while

(5.27) $I_{5}(\Psi)=-\tfrac{6^{5}\cdot 5^{2}\cdot
199\sqrt{5}}{154\pi\sqrt{\pi}}\,.$

Due to (5.11) with $k=5$, and taking into account (5.24) and (5.23), from
(5.26)-(5.27) we get

(5.28) $[c_{N,2}^{1}(t_{N})]^{2}-4\,|c_{N,2}^{2}(t_{N})|^{2}=80\cdot
199\delta^{2}+{\rm O}(\delta^{3})+{\rm O}(\sqrt{\varepsilon})\,.$

But (5.23), (5.25) and (5.28) clearly cannot hold simultaneously. The obtained
contradiction proves that the zonal flow (5.1) is stable also for $\alpha=0$.

$(ii)$ Let

$\widehat{\psi}^{n}_{0}=\psi_{0}+\tfrac{1}{n}\,\sin\theta=\big{(}\alpha+\tfrac{1}{n}\big{)}\,\sin\theta+\beta\,Y(\varphi,\theta)$

with $Y\in{\mathbb{E}}_{j}$. According to (2.12), the solution
$\widehat{\psi}^{n}(t)$ of (${\mathcal{E}}_{\omega}$) with initial data
$\widehat{\psi}^{n}(0)=\widehat{\psi}^{n}_{0}$ is the travelling wave

$\widehat{\psi}^{n}(t)=\big{(}\alpha+\tfrac{1}{n}\big{)}\,\sin\theta+\beta\,Y(\varphi-\widehat{c}\,t,\theta)$

with propagation speed
$\widehat{c}=\frac{2\omega}{j(j+1)}+\left(\alpha+\frac{1}{n}\right)\,\frac{j(j+1)-2}{j(j+1)}$.
On the other hand,

$\psi_{0}(t)=\alpha\,\sin\theta+\beta\,Y(\varphi-ct,\theta)$

with $c=\frac{2\omega}{j(j+1)}+\alpha\,\frac{j(j+1)-2}{j(j+1)}$. If (see the
Appendix)

$Y(\varphi,\theta)=\sum_{|m|\leq j}c_{j}\,P_{j}^{m}(\sin\theta)\,{\rm e}^{{\rm
i}m\varphi}\,,$

we get

$\displaystyle\sup_{t>0}\|\widehat{\psi}^{n}(t)-\psi(t)\|_{L^{2}({\mathbb{S}}^{2},{\rm
d}\sigma)}^{2}$
$\displaystyle=\sup_{t>0}\left\|\tfrac{1}{n}\,\sin\theta+\beta\,Y(\varphi-\widehat{c}\,t,\theta)-\beta\,Y(\varphi-
ct,\theta)\right\|_{L^{2}({\mathbb{S}}^{2},{\rm d}\sigma)}^{2}$
$\displaystyle=\sup_{t>0}\left\|\tfrac{1}{n}\,\sin\theta+\beta\,\sum_{|m|\leq
j}c_{j}\,P_{j}^{m}(\sin\theta)\,\Big{(}{\rm e}^{{\rm
i}m(\varphi-\widehat{c}\,t)}-{\rm e}^{{\rm
i}m(\varphi-c\,t)}\Big{)}\right\|_{L^{2}({\mathbb{S}}^{2},{\rm d}\sigma)}^{2}$
$\displaystyle=\sup_{t>0}\Big{\\{}\frac{4\pi}{3n}+\beta^{2}\sum_{0<|m|\leq
j}|c_{j}|^{2}\Big{(}\int_{-\tfrac{\pi}{2}}^{\tfrac{\pi}{2}}\Big{|}P_{j}^{m}(\sin\theta)\Big{|}^{2}\cos\theta\,{\rm
d}\theta\Big{)}\Big{(}\int_{0}^{2\pi}\Big{|}{\rm e}^{{\rm
i}m(\varphi-\widehat{c}\,t)}-{\rm e}^{{\rm i}m(\varphi-c\,t)}\Big{|}^{2}{\rm
d}\varphi\Big{)}\Big{\\}}$
$\displaystyle=\sup_{t>0}\Big{\\{}\frac{4\pi}{3n}+2\pi\beta^{2}\sum_{0<|m|\leq
j}|c_{j}|^{2}\Big{(}\int_{-\tfrac{\pi}{2}}^{\tfrac{\pi}{2}}\Big{|}P_{j}^{m}(\sin\theta)\Big{|}^{2}\cos\theta\,{\rm
d}\theta\Big{)}\Big{|}1-{\rm e}^{{\rm
i}m(\widehat{c}-c)\,t}\Big{|}^{2}\Big{\\}}$
$\displaystyle=\sup_{t>0}\Big{\\{}\frac{4\pi}{3n}+4\pi\beta^{2}\sum_{0<|m|\leq
j}|c_{j}|^{2}\Big{(}\int_{-\tfrac{\pi}{2}}^{\tfrac{\pi}{2}}\Big{|}P_{j}^{m}(\sin\theta)\Big{|}^{2}\cos\theta\,{\rm
d}\theta\Big{)}\Big{(}1-\cos\big{(}m\,\tfrac{j(j+1)-2}{nj(j+1)}\,t\big{)}\Big{)}\Big{\\}}$
$\displaystyle=\frac{4\pi}{3n}+4\pi\beta^{2}\sum_{0<|m|\leq
j}|c_{j}|^{2}\Big{(}\int_{-\tfrac{\pi}{2}}^{\tfrac{\pi}{2}}\Big{|}P_{j}^{m}(\sin\theta)\Big{|}^{2}\cos\theta\,{\rm
d}\theta\Big{)}$ $\displaystyle>4\pi\beta^{2}\sum_{0<|m|\leq
j}|c_{j}|^{2}\Big{(}\int_{-\tfrac{\pi}{2}}^{\tfrac{\pi}{2}}\Big{|}P_{j}^{m}(\sin\theta)\Big{|}^{2}\cos\theta\,{\rm
d}\theta\Big{)}$

for all $n\geq 1$.

$(iii)$ Inspired by the approach used in the proof of Theorem 6, we define

$\mathcal{E}(\psi)=\iint_{{\mathbb{S}}^{2}}\Big{[}\tfrac{1}{2}\,|U|^{2}-\tfrac{1}{12}\,|\Omega+2\omega\sin\theta|^{2}+\omega\sin\theta\,\Omega-|\mathbb{P}_{1}\psi|^{2}\Big{]}\,{\rm
d}\sigma\,.$

This functional is constant along solutions $\psi$ of
(${\mathcal{E}}_{\omega}$). Expanding a solution $\psi=\psi_{0}+\delta\psi$ of
(${\mathcal{E}}_{\omega}$) in terms of spherical harmonics,

$\begin{cases}\psi_{0}(\varphi,\theta)=-\omega\sqrt{\tfrac{\pi}{3}}\,Y_{1}^{0}(\theta)+\beta_{0}\displaystyle\sum_{|m|\leq
2}a_{m}\,Y_{2}^{m}(\varphi,\theta)\quad\text{with}\quad\sum_{|m|\leq
2}|a_{m}|^{2}=1\,,\\\ \delta\psi(\varphi,\theta,t)=\displaystyle\sum_{j\geq
1}\Big{\\{}\sum_{|m|\leq
j}c_{j}^{m}(t)\,Y_{j}^{m}(\varphi,\theta)\Big{\\}}\quad\text{with}\quad
c_{1}^{0}(t)=c_{1}^{0}\quad\text{and}\quad c_{1}^{\pm 1}(t)=c_{1}^{\pm
1}(0)\,{\rm e}^{\mp{\rm i}\omega t}\quad\text{for}\quad t\geq 0\,,\end{cases}$

we can write (see the Appendix)

$\displaystyle\mathcal{E}(\psi)$
$\displaystyle=-\tfrac{1}{3}\,\Big{\\{}\Big{(}c_{1}^{0}(0)-\omega\sqrt{\tfrac{\pi}{3}}\Big{)}^{2}+|c_{1}^{-1}(0)|^{2}+|c_{1}^{1}(0)|^{2}\Big{\\}}+\tfrac{1}{12}\sum_{j\geq
3}\Big{\\{}\sum_{|m|\leq j}|c_{j}^{m}(t)|^{2}[6-j(j+1)]\,j(j+1)\Big{\\}}$
$\displaystyle\leq-\tfrac{1}{3}\,\Big{\\{}\Big{(}c_{1}^{0}(0)-\omega\sqrt{\tfrac{\pi}{3}}\Big{)}^{2}+|c_{1}^{-1}(0)|+|c_{1}^{1}(0)|\Big{\\}}-\tfrac{1}{2}\sum_{j\geq
3}j(j+1)\Big{\\{}\sum_{|m|\leq j}|c_{j}^{m}(t)|^{2}\Big{\\}}\,.$

Exploiting the conservation of the kinetic energy, which ensures that the
expression

$\Big{(}c_{1}^{0}(0)-\omega\sqrt{\tfrac{\pi}{3}}\Big{)}^{2}+|c_{1}^{-1}(0)|^{2}+|c_{1}^{1}(0)|^{2}+3\sum_{|m|\leq
2}|c_{2}^{m}(t)+\beta_{0}a_{m}|^{2}+\sum_{j\geq
3}\tfrac{j(j+1)}{2}\Big{\\{}\sum_{|m|\leq j}|c_{j}^{m}(t)|^{2}\Big{\\}}$

is time-independent, setting

$\alpha=\tfrac{1}{2}\,c_{1}^{0}(0)\,\sqrt{\tfrac{\pi}{3}}\,,\quad\beta(t)=\sqrt{\sum_{|m|\leq
2}|c_{2}^{m}(t)+\beta_{0}a_{m}|^{2}}\,,\quad Y(t)=\sum_{|m|\leq
2}\tfrac{c_{2}^{m}(t)+\beta_{0}a_{m}}{\beta(t)}\,Y_{2}^{m}\,,\quad\tilde{\psi}(t)=\sum_{j\geq
3}\Big{\\{}\sum_{|m|\leq j}c_{j}^{m}(t)\,Y_{j}^{m}\Big{\\}}\,,$

proves the claim.∎

##  6\. Bifurcation from Rossby-Haurwitz waves

This section is dedicated to constructing stationary and travelling-wave
solutions of (${\mathcal{E}}_{\omega}$) which are different from the explicit
solutions of (1.1) studied above. For this, we seek non-zonal solutions to a
suitably modified form of equation (1.1) which bifurcate from Rossby-Haurwitz
waves.

###  6.1. The case $\omega=0$.

To implement a bifurcation approach, it is convenient to introduce the
parameter $\lambda\in{\mathbb{R}}$, seeking solutions $\psi\in
C^{2,\alpha}(\mathbb{S}^{2})$ of the nonlinear elliptic equation

(6.1) $-\Delta\psi+F(\lambda,\psi)=0$

for functions $F\in C^{2}({\mathbb{R}}^{2},{\mathbb{R}})$. Note that a
solution $\psi\in C^{2,\alpha}(\mathbb{S}^{2})$ to (6.1) satisfies the Gauss
constraint

(6.2) $\iint_{{\mathcal{S}}^{2}}F(\lambda,\psi)\,{\rm d}\sigma=0\,.$

On the other hand, any solution $\psi\in C^{2,\alpha}(\mathbb{S}^{2})$ of

(6.3)
$-\Delta\psi+F(\lambda,\psi)=\frac{1}{4\pi}\iint_{{\mathcal{S}}^{2}}F(\lambda,\psi)\,{\rm
d}\sigma$

provides us with a stationary solution of (${\mathcal{E}}_{\omega}$). Note
that the zero function $\psi\equiv 0$ solves (6.3), and the vorticity
$\Delta\psi$ of a solution $\psi\in C^{2,\alpha}(\mathbb{S}^{2})$ to (6.3)
satisfies the Gauss constraint (2.6). For this reason, rather than solving
(6.1) with the constraint (6.2), we will seek solutions to (6.3).

Since the Laplace-Beltrami operator is bijective from
$C^{2,\alpha}(\mathbb{S}^{2})$ onto

$C^{0,\alpha}_{0}(\mathbb{S}^{2})=\Big{\\{}f\in C^{0,\alpha}(\mathbb{S}^{2}):\
\iint_{{\mathcal{S}}^{2}}f\,{\rm d}\sigma=0\Big{\\}}\,,$

with a compact inverse that we denote by $\Delta^{-1}$, we can recast equation
(6.3) in the form

(6.4)
${\mathbb{F}}(\lambda,f)=0\quad\text{with}\quad{\mathbb{F}}:{\mathbb{R}}\times
C^{0,\alpha}_{0}(\mathbb{S}^{2})\to
C^{0,\alpha}_{0}(\mathbb{S}^{2})\,,\quad{\mathbb{F}}(\lambda,f)=f-\Delta^{-1}\Big{\\{}F(\lambda,f)-\frac{1}{4\pi}\iint_{{\mathcal{S}}^{2}}F(\lambda,f)\,{\rm
d}\sigma\Big{\\}}\,.$

Invariant spherical harmonic basis functions, tabulated by their respective
subgroups of ${\mathbb{O}}(3)$, are listed in [26]. For a finite subgroup
${\mathbb{G}}$ of ${\mathbb{O}}(3)$ with the property that the subspace of
${\mathbb{G}}$-invariant spherical harmonics of degree $n\geq 1$ is one-
dimensional, the fact that (6.3) is equivariant with respect to the natural
action of the orthogonal group enables us to consider this problem restricted
to ${\mathbb{G}}$-equivariant functions. This way, we can take advantage of
the symmetries to analyse the formation of regular flow patterns using the
Rabinowitz global bifurcation approach (see [35, 44]).

######  Lemma 1 (The Rabinowitz global bifurcation theorem).

Let ${\mathbb{Y}}$ be a real Banach space and let ${\mathbb{F}}\in
C^{2}({\mathbb{R}}\times{\mathbb{Y}},{\mathbb{Y}})$ be such that $f\mapsto
f-{\mathbb{F}}(\lambda,f)$ is compact operator from
${\mathbb{R}}\times{\mathbb{Y}}$ to ${\mathbb{Y}}$ and:

(i) ${\mathbb{F}}(\lambda,0)=0$ for all
$(\lambda,0)\in{\mathbb{R}}\times{\mathbb{Y}}$;

(ii)
$\partial_{f}{\mathbb{F}}(\lambda^{\ast},0)\in{\mathcal{L}}({\mathbb{Y}},{\mathbb{Y}})$
is a linear Fredholm operator of index zero and one-dimensional kernel
${\mathcal{N}}(\partial_{f}{\mathbb{F}}(\lambda^{\ast},0))$ generated by some
$f^{\ast}\in{\mathbb{Y}}\setminus\\{0\\}$;

(iii) the transversality condition holds, in the sense that
$[\partial_{\lambda,f}^{2}\,{\mathbb{F}}(\lambda^{\ast},0)]\,(1,f^{\ast})$
does not belong to the range of the operator
$\partial_{f}{\mathbb{F}}(\lambda^{\ast},0))$, where
$\partial^{2}_{\lambda,f}{\mathbb{F}}(\lambda^{\ast},0)=\partial_{\lambda}[\partial_{f}{\mathbb{F}}(\lambda,0)]\big{|}_{\lambda=\lambda^{\ast}}\in{\mathcal{L}}({\mathbb{R}},{\mathcal{L}}({\mathbb{Y}},{\mathbb{Y}}))={\mathcal{L}}({\mathbb{R}}\times{\mathbb{Y}},{\mathbb{Y}})$.

Then there exist $\varepsilon>0$, an open set
${\mathcal{O}}\subset{\mathbb{R}}\times{\mathbb{Y}}$ with
$(\lambda^{\ast},0)\in{\mathcal{O}}$ and a branch of solutions

$\\{(\lambda,f)=(\lambda(s),\,s\,\chi(s)):\
s\in{\mathbb{R}},\,|s|<\varepsilon\\}\subset{\mathbb{R}}\times{\mathbb{Y}}$

of ${\mathbb{F}}(\lambda,f)=0$ with $\lambda(0)=\lambda^{\ast}$,
$\chi(0)=f^{\ast}$, and such that $s\mapsto\lambda(s)\in{\mathbb{R}}$,
$s\mapsto s\chi(s)\in{\mathbb{Y}}$ are continuously differentiable on
$(-\varepsilon,\varepsilon)$, with

$\\{(\lambda,f)\in{\mathcal{O}}:\ {\mathbb{F}}(\lambda,f)=0,\ f\neq
0\\}=\\{(\lambda(s),\,s\,\chi(s)):\ 0<|s|<\varepsilon\\}.$

Furthermore, if ${\mathcal{S}}$ is the closure of the set of nontrivial
solutions of ${\mathbb{F}}(\lambda,f)=0$ in ${\mathbb{R}}\times{\mathbb{Y}}$,
then the connected component ${\mathcal{S}}^{\ast}$ of ${\mathcal{S}}$ to
which $(\lambda^{\ast},0)$ belongs has at least one of the following
properties:

(I) ${\mathcal{S}}^{\ast}$ is unbounded in ${\mathbb{R}}\times{\mathbb{Y}}$;

(II) there exists some $\lambda\neq\lambda^{\ast}$ such that
$(\lambda,0)\in{\mathcal{S}}^{\ast}$.

We now prove the following existence result.

######  Theorem 8.

Let $G$ be a finite subgroup ${\mathbb{G}}$ of ${\mathbb{O}}(3)$ with the
property that the subspace of ${\mathbb{G}}$-invariant spherical harmonics of
some specific degree $\ell\geq 3$ is one-dimensional and non-zonal. If
$F:{\mathbb{R}}^{2}\to{\mathbb{R}}$ is twice continuously differentiable and
such that

$F_{f}(\lambda^{*},0)=\ell(\ell+1),\qquad F_{f\lambda}(\lambda^{*},0)\neq 0,$

then there exists a maximal connected component of nontrivial solutions of
(6.3) such that all solutions $(\lambda,\psi)\in{\mathcal{S}}^{\ast}$ close
enough to $(\lambda^{\ast},0)$ are non-zonal.

Figure 4. Visualization of the non-axial symmetries of specific spherical
harmonics, using different colours to keep track of the geometry of different
level sets. Note that the subspace of spherical harmonics of degree 3 that are
invariant under the finite tetrahedral subgroup (the symmetry group of the
methane molecule) is generated by $Y_{3}^{-2}$, while the subspace of
spherical harmonics of degree 5 that are invariant under the finite subgroup
$D_{4}^{d}$ (the symmetry group of the bicapped square antiprism, describing
the shape of the octasulfur molecule) is generated by $Y_{5}^{-2}$.

###### Proof.

Consider the map ${\mathbb{F}}:{\mathbb{R}}\times{\mathbb{Y}}\to{\mathbb{Y}}$
defined as in (6.4) by

(6.5)
${\mathbb{F}}(\lambda,f)=f-\Delta^{-1}\Big{\\{}F(\lambda,f)-\frac{1}{4\pi}\iint_{{\mathcal{S}}^{2}}F(\lambda,f)\,{\rm
d}\sigma\Big{\\}}\,$

where the Banach space

${\mathbb{Y}}=\\{f\in C^{0,\alpha}_{0}({\mathbb{S}}^{2}):\ Gf=f\quad\text{for
all}\quad G\in{\mathbb{G}}\\}$

captures the symmetries associated with the group ${\mathbb{G}}$. The
linearization of the operator ${\mathbb{F}}$ about the trivial solution $f=0$
of ${\mathbb{F}}(\lambda,f)=0$ is

(6.6)
$\partial_{f}{\mathbb{F}}(\lambda,0)[\psi]=\psi-\Delta^{-1}\Big{\\{}F_{f}(\lambda,0)\psi-\frac{1}{4\pi}\iint_{{\mathcal{S}}^{2}}F_{f}(\lambda,0)\psi\,{\rm
d}\sigma\Big{\\}}=\psi-
F_{f}(\lambda,0)\,\Delta^{-1}\psi\,,\qquad\psi\in\mathbb{Y}\,.$

For $F_{f}(\lambda,0)=-\ell(\ell+1)$ this operator acting on
$C^{0,\alpha}_{0}({\mathbb{S}}^{2})$ has a nontrivial kernel given by the
$(2\ell+1)$-dimensional space of spherical harmonics of order $\ell$.
Consequently, if the finite subgroup ${\mathbb{G}}$ of ${\mathbb{O}}(3)$ has
the property that the subspace of ${\mathbb{G}}$-invariant spherical harmonics
of degree $\ell$ is one-dimensional, then for $\lambda$ such that
$F_{f}(\lambda,0)=-\ell(\ell+1)$, it follows that the kernel
${\mathcal{N}}(\partial_{f}{\mathbb{F}}(\lambda,0))$ of the operator
${\mathbb{F}}$ defined in (6.5) is one-dimensional, being generated by some
$f^{*}\in\mathbb{E}_{\ell}$. Note that for a simple eigenvalue, the condition
that $\partial_{f}{\mathbb{F}}(\lambda^{\ast},0)$ is a Fredholm operator of
index zero means that the range of this operator is closed and has a one-
dimensional complement. If $\lambda^{\ast}$ is such that
$F_{f}(\lambda^{\ast},0)=-\ell(\ell+1)$, then, due to elliptic regularity and
to the self-adjointness of the Laplace-Beltrami operator in
$L^{2}(\mathbb{S}^{2})$, we know that $\psi\in C^{0,\alpha}(\mathbb{S}^{2})$
belongs to the range of $\partial_{f}{\mathbb{F}}(\lambda^{\ast},0)$ if and
only if it belongs to the orthogonal complement in $L^{2}(\mathbb{S}^{2})$ of
the spherical harmonics of degree $\ell$. From this it follows at once that
$\partial_{f}{\mathbb{F}}(\lambda^{\ast},0)$ acting on ${\mathbb{Y}}$ has a
closed range with a one-dimensional complement in ${\mathbb{Y}}$. Furthermore,
since

$[\partial^{2}_{\lambda,f}{\mathbb{F}}(\lambda^{\ast},0)]\,(1,f^{\ast})=-\,F_{f\lambda}(\lambda^{\ast},0)\,\Delta^{-1}f^{\ast}\,,$

we see that the transversality condition is equivalent to
$F_{f\lambda}(\lambda^{\ast},0)\neq 0$. All the hypotheses in Lemma 1 hold,
and we deduce the existence of a curve $\\{(\lambda(s),\,s\,\chi(s)):\
s\in{\mathbb{R}},\,|s|<\varepsilon\\}$ of nontrivial solutions that bifurcates
at $(\lambda^{\ast},0)$ from the curve $\\{(\lambda,0):\
\lambda\in{\mathbb{R}}\\}$ of trivial solutions. Since the tangent vector of
this nontrivial solution curve at the bifurcation point is given by
$(\lambda^{\prime}(0),f^{\ast})$ and $f^{\ast}$ is non-zonal, near the
bifurcation point all nontrivial solutions are non-zonal. ∎

The next result describe settings in which it is possible to reveal structural
properties of the continuum of solutions found in Theorem 8.

######  Theorem 9.

If $P:{\mathbb{R}}\to{\mathbb{R}}$ is twice continuously differentiable and

  * •

$\lim_{\lambda\to\infty}P(\lambda)=\infty$ and
$\lim_{\lambda\to-\infty}P(\lambda)=-\infty$,

  * •

there exists $\lambda^{\ast}\in{\mathbb{R}}$ with
$P^{\prime}(\lambda^{\ast})=-l(l+1)$ and $P^{\prime\prime}(\lambda^{\ast})\neq
0$,

  * •

there exists $a>0$ with $P^{\prime}(\lambda)>0$ for $|\lambda|>a$,

then the maximal connected component built up in Theorem 8 for
$F(\lambda,f)=P(\lambda+f)-P(\lambda)$ has a closure ${\mathcal{S}}^{\ast}$
which is bounded in ${\mathbb{R}}\times{\mathbb{Y}}$, contains
$(\lambda^{\ast},0)$ as well as some $(\lambda,0)$ with
$\lambda\neq\lambda^{\ast}$.

###### Proof.

We have to prove that all nontrivial solutions $(\lambda,\psi)$ of
${\mathbb{F}}(\lambda,\psi)=0$ with ${\mathbb{F}}$ defined by (6.5) are
bounded a priori in ${\mathbb{R}}\times{\mathbb{Y}}$, or, equivalently, that
all nontrivial solutions of (6.3) are bounded a priori in ${\mathbb{R}}\times
C^{2,\alpha}_{0}({\mathbb{S}}^{2})$. Note that for $\psi\in
C^{2,\alpha}({\mathbb{S}}^{2})$, $\psi\not\equiv 0$, we have

$\psi_{m}=\inf_{\mathbb{S}^{2}}\\{\psi\\}<0<\sup_{\mathbb{S}^{2}}\\{\psi\\}=\Psi_{M}$

since $\iint_{\mathbb{S}^{2}}\psi\,{\rm d}\sigma=0$. For a function $P$ with
the properties specified in the statement one can easily see that there exist
$a_{+}>0$ and $a_{-}<0$ such that

  * •

$P(a_{+})>0>P(a_{-})$ and $P^{\prime}(\lambda)>0$ for
$\lambda\in(-\infty,a_{-})\cup(a_{+},\infty)$,

  * •

$\sup_{\lambda\in[a_{-},a_{+}]}\\{|P(\lambda)|\\}\leq
A=\max\\{-P(a_{-}),\,P(a_{+})\\}$.

If $(\lambda,\psi)\in{\mathbb{R}}\times{\mathbb{Y}}$ with $\psi\not\equiv 0$
lies on the continuum of solutions, the weak maximum principle applied to
(6.3) yields

$F(\lambda,\psi_{M})=P(\lambda+\psi_{M})-P(\lambda)\leq\frac{1}{4\pi}\iint_{{\mathcal{S}}^{2}}F(\lambda,\psi)\,{\rm
d}\sigma\leq P(\lambda+\psi_{m})-P(\lambda)=F(\lambda,\psi_{m})\,.$

Thus $P(\lambda+\psi_{m})\geq P(\lambda+\psi_{M})$ and since
$\lambda+\psi_{m}<\lambda<\lambda+\psi_{M}$ we deduce that all three points
belong to the interval $[a_{-},a_{+}]$. Consequently we have the following
$L^{\infty}(\mathbb{S}^{2})$ a priori bounds for any nontrivial solution
$\psi\in C^{2,\alpha}_{0}({\mathbb{S}}^{2})$ of equation (6.3):

(6.7) $|\lambda|+|\psi|\leq 2(a_{+}-a_{-})\quad\text{and}\quad|\Delta\psi|\leq
2A\,.$

We now invoke for $k>2$ the $L^{k}$-estimates for elliptic equations on smooth
compact manifolds without boundary (see [6]): there are constants $c_{0}>0$
and $c_{2}>0$ such that

$\|\psi\|_{W^{2,k}(\mathbb{S}^{2})}\leq
c_{2}\|\Delta\psi\|_{L^{k}(\mathbb{S}^{2})}+c_{0}\|\psi\|_{L^{1}(\mathbb{S}^{2})}$

for every $\psi\in W^{2,k}(\mathbb{S}^{2})$. In conjunction with (6.7) and
with the Sobolev embedding $W^{2,k}(\mathbb{S}^{2})\subset
C^{1}(\mathbb{S}^{2})$, we obtain a priori bounds in $C^{1}(\mathbb{S}^{2})$
for the nontrivial solutions $\psi\in C^{2,\alpha}_{0}({\mathbb{S}}^{2})$ of
(6.3). But then (6.3) yields by differentiation a priori bounds for
$\|\Delta\psi\|_{C^{0,\alpha}(\mathbb{S}^{2})}$. We now conclude the
boundedness of the nontrivial solutions of (6.3) in ${\mathbb{R}}\times
C^{2,\alpha}_{0}({\mathbb{S}}^{2})$ from the Schauder estimates (see [6]):
there are constants $C_{0}>0$ and $C_{2}>0$ such that

$\|\psi\|_{C^{2,\alpha}(\mathbb{S}^{2})}\leq
C_{2}\|\Delta\psi\|_{C^{0,\alpha}(\mathbb{S}^{2})}+C_{0}\|\psi\|_{C(\mathbb{S}^{2})}$

for every $\psi\in C^{2,\alpha}(\mathbb{S}^{2})$.∎∎

######  Remark 2.

Simple examples of functions satisfying the hypotheses of Theorem 8 are
provided by the polynomials

$P(\lambda)=\mu_{1}\lambda^{3}-[\mu+l(l+1)]\lambda$

with parameters $\mu>0$ and $\mu_{1}>0$.∎

###  6.2. The case $\omega>0$

We now establish the existence of a global continuum of solutions to (1.1) for
$\omega>0$ and suitable twice continuously differentiable functions
$F:{\mathbb{R}}\to{\mathbb{R}}$. To take advantage of the fact that the
spherical harmonics provide a representation of the orthogonal group
${\mathbb{O}}(3)$, we write (1.1) in the form

(6.8) $\Delta\psi(\xi)=F(\psi(\xi))-2\omega z\,,\qquad\xi\in\mathbb{S}^{2}\,,$

where $z$ is the distance of $\xi\in\mathbb{S}^{2}$ to the equatorial plane.
An adequate class of nonlinear functions $F$ is obtained by modifying outside
a neighbourhood of zero the linear functions $\psi\mapsto-l(l+1)\psi$ with
$l\in{\mathbb{N}}$, aiming at replacing the solution set

$\\{\omega
z+\sum_{m=-l}^{l}c_{m}Y_{l}^{m}(\varphi,\theta):\,c_{m}\in{\mathbb{R}}\
\text{for}\ -l\leq m\leq l\\}$

which features functions with gradients of all possible sizes by a continuum
of nontrivial solutions with an a priori bound on the corresponding velocity
fields and vorticities.

######  Theorem 10.

Let $G$ be a finite subgroup ${\mathbb{G}}$ of ${\mathbb{O}}(3)$ with the
property that the subspace of ${\mathbb{G}}$-invariant spherical harmonics of
some specific degree $l\geq 3$ is one-dimensional and non-zonal, being
generated by some $f^{\ast}\neq 0$. Given $\beta>0$, if
$P:{\mathbb{R}}\to{\mathbb{R}}$ is twice continuously differentiable and

  * •

there exists $\mu>\frac{2\nu}{l(l+1)}$ with
$P(\lambda)=-\frac{2\nu}{\mu}\,\lambda$ for $|\lambda|\leq 2\mu$, where
$\nu=\frac{\beta l(l+1)}{l(l+1)-2}>\beta$,

  * •

$\lim_{\lambda\to-\infty}P(\lambda)=-\infty$,
$\lim_{\lambda\to\infty}P(\lambda)=\infty$ and there exists $a>2\mu$ with
$P(\lambda)>0$ for $|\lambda|>a$,

then there exists a maximal connected component ${\mathcal{S}}_{\beta}$ of
nontrivial solutions of the vorticity equation

(6.9)
$\Delta\psi=P((1+\lambda^{2})\psi)-2\Big{(}\nu-\frac{\mu}{1+\lambda^{2}}\Big{)}\,z-\frac{1}{4\pi}\,\iint_{{\mathbb{S}}^{2}}P((1+\lambda^{2})\psi)\,{\rm
d}\sigma\,,$

such that the corresponding velocity fields are uniformly bounded a priori.
Moreover, the continuum ${\mathcal{S}}_{\beta}$ comprises non-zonal Rossby
solutions of the form

$\psi=c^{\ast}f^{\ast}-\frac{\mu}{1+\lambda^{2}}\,z\quad\textit{with}\quad
c^{\ast}\neq 0$

close to the zonal solution $\psi_{0}=-\frac{2\nu}{l(l+1)}\,z$ of the linear
equation

$\Delta\psi-l(l+1)\psi+2\beta z=0\,.$

###### Proof.

Let us first note that if $f$ solves

(6.10) $\Delta f-P\Big{(}(1+\lambda^{2})f-\mu z\Big{)}+2\nu
z+\frac{1}{4\pi}\,\iint_{{\mathbb{S}}^{2}}P\Big{(}(1+\lambda^{2})f-\mu
z\Big{)}\,{\rm d}\sigma=0\,,$

then

(6.11) $\psi=f-\frac{\mu}{1+\lambda^{2}}\,z$

solves (6.9). We will therefore develop a global bifurcation approach to
establish the existence of nontrivial solutions of (6.10). Using the Banach
space ${\mathbb{Y}}$ introduced in the proof of Theorem 8, we transform (6.10)
to ${\mathbb{F}}(\lambda,f)=0$, where

${\mathbb{F}}:{\mathbb{R}}\times{\mathbb{Y}}\to{\mathbb{Y}}\,,\qquad{\mathbb{F}}(\lambda,f)=f-\Delta^{-1}\Big{\\{}P\Big{(}(1+\lambda^{2})f-\mu
z\Big{)}-2\nu z\Big{\\}}\,.$

Note that ${\mathbb{F}}(\lambda,0)=0$. Since $|z|\leq 1$ we have

$\partial_{f}{\mathbb{F}}(\lambda,0)[f_{0}]=f_{0}-(1+\lambda^{2})\Delta^{-1}\Big{\\{}P^{\prime}(-\mu
z)f_{0}-\frac{1}{4\pi}\,\iint_{{\mathbb{S}}^{2}}P^{\prime}(-\mu z)f_{0}\,{\rm
d}\sigma\Big{\\}}=f_{0}+\frac{2\nu(1+\lambda^{2})}{\mu}\,\Delta^{-1}f_{0}\,,$

so that a necessary condition for the kernel of
$\partial_{f}{\mathbb{F}}(\lambda,0)$ to comprises more than $\\{0\\}$ is

$\frac{2\nu(1+\lambda^{2})}{\mu}=n(n+1)$

for some $n\in{\mathbb{N}}$. The case $n=l$ corresponds to

(6.12) $\pm\lambda^{\ast}=\sqrt{\frac{\mu l(l+1)}{2\nu}-1}$

and in this case the kernel is one-dimensional, being generated by $f^{\ast}$.
Since

$[\partial^{2}_{\lambda
f}{\mathbb{F}}(\lambda^{\ast}),0](1,f^{\ast})=\frac{4\lambda^{\ast}\nu}{\mu}\Delta^{-1}f^{\ast}\quad\text{with}\quad\frac{4\lambda^{\ast}\nu}{\mu}\neq
0\,,$

a reasoning analogous to that in the proof of Theorem 8 ensures the existence
of a global continuum of nontrivial solutions that bifurcate at
$(\lambda^{\ast},0)$ from the curve of trivial solutions. Performing the
transformation (6.11) we obtain a corresponding continuum of non-trivial
solutions to (6.9). Close to the bifurcation point, as long as
$|(1+\lambda^{2})f-\mu z|\leq 2\mu$ throughout $\mathbb{S}^{2}$, the equation
(6.10) takes the form

(6.13) $\Delta f+\frac{2\nu}{\mu}\,(1+\lambda^{2})f=0\,.$

Equation (6.13) has a nontrivial solution only if
$\frac{2\nu}{\mu}\,(1+\lambda^{2})=n(n+1)$ for some $n\in{\mathbb{N}}$. But
for $\lambda=\lambda^{\ast}$, equation (6.13) takes the form $\Delta
f+l(l+1)f=0$, due to (6.12). Consequently, for as long as
$|(1+\lambda^{2})f-\mu z|\leq 2\mu$ throughout $\mathbb{S}^{2}$, we must have
$\lambda=\lambda^{\ast}$ for the corresponding non-trivial solutions in the
continuum and the non-trivial function $\psi$ is a spherical harmonic of
degree $l$ and thus a multiple of $f^{\ast}$, as claimed. Moreover, taking
into account the two alternatives for the continuum of non-trivial bifurcating
solutions, we infer the existence of non-trivial solutions in this continuum
satisfying $\|(1+\lambda^{2})f-\mu z\|_{L^{\infty}(\mathbb{S}^{2})}>2\mu$,
solutions for which the nonlinear adjustment of $P$ comes into play.

It remains to prove that the gradients of solutions to (6.9) within the above
continuum are uniformly bounded a priori. As in the proof of Theorem 8, this
follows at once from elliptic $L^{k}$ estimates with $k>2$ in conjunction with
the Sobolev embedding $W^{2,k}(\mathbb{S}^{2})\subset C^{1}(\mathbb{S}^{2})$,
once we establish an $L^{\infty}(\mathbb{S}^{2})$ a priori bound for the non-
trivial solutions $\psi\in C^{2,\alpha}_{0}({\mathbb{S}}^{2})$ of (6.9). For
this, choose $b>a$ such that $P(-b)<0<P(b)$ and

$\min\\{P(b),\,|P(-b)|\\}\geq\max_{|\lambda|\leq a}\\{|P(\lambda)|\\}\,.$

Let us now note that if $(\lambda,\psi)\in{\mathbb{R}}\times
C^{2,\alpha}_{0}({\mathbb{S}}^{2})$ with $\psi\not\equiv 0$ solves (6.9), then

$\psi_{m}=\inf_{\mathbb{S}^{2}}\\{\psi\\}<0<\sup_{\mathbb{S}^{2}}\\{\psi\\}=\Psi_{M}$

since $\iint_{\mathbb{S}^{2}}\psi\,{\rm d}\sigma=0$. From the weak maximum
principle in conjunction with (6.10) we obtain

$P((1+\lambda^{2})\psi_{M}-\mu z)-2\nu
z\leq\frac{1}{4\pi}\,\iint_{{\mathbb{S}}^{2}}P\Big{(}(1+\lambda^{2})f-\mu
z\Big{)}\,{\rm d}\sigma\leq P((1+\lambda^{2})\psi_{m}-\mu z)-2\nu z\,,$

so that

$P((1+\lambda^{2})\psi_{M}-\mu z)\leq P((1+\lambda^{2})\psi_{m}-\mu z)\,.$

Since $\psi_{m}<\psi_{M}$, we must have

$-b\leq(1+\lambda^{2})\psi_{m}-\mu z\leq(1+\lambda^{2})\psi_{M}-\mu z\leq
b\,.$

But then $|z|\leq 1$ yields

$-(b+\mu)\leq(1+\lambda^{2})\psi_{m}\leq(1+\lambda^{2})\psi_{M}\leq b+\mu\,,$

and therefore $\|\psi\|_{L^{\infty}(\mathbb{S}^{2})}\leq\mu+b$. This completes
the proof.∎∎

##  7\. Relevance for stratospheric flows

In this section we show that some of the inviscid flows on a rotating sphere
studied hitherto are building blocks for the leading-order dynamics of 3D
stratospheric flows.

In the stratosphere the atmospheric flow is practically inviscid (see [11]),
being thus governed by the components of the Euler equation (see [25])

(7.1a)