A sufficient condition for inviscid shear instability: hurdle theorem and its application to alternating jets

Abstract

We propose a simple method to identify unstable parameter regions in general inviscid unidirectional shear flow stability problems. The theory is applicable to a wide range of basic flows, including those that are non-monotonic. We illustrate the method using a model of Jupiter's alternating jet streams based on the quasi-geostrophic equation. The main result is that the flow is unstable if there is an interval in the flow domain for which the reciprocal Rossby Mach number (a quantity defined in terms of the zonal flow and potential vorticity distribution), surpasses a certain threshold or ‘hurdle’. The hurdle height approaches unity when we can take the hurdle width to greatly exceed the atmosphere's intrinsic deformation length, as holds on gas giants. In this case, the Kelvin–Arnol’d sufficient condition of stability accurately detects instability. These results improve the theoretical framework for explaining the stable maintenance of Jupiter and Saturn's jets over decadal time scales.

1. Introduction

The stability of inviscid, parallel shear flows has applications in geophysical fluid dynamics, astrophysics, plasma physics and engineering thermofluid sciences. For the purely hydrodynamic cases, sufficient conditions for stability can be traced back to the Rayleigh (1880) inflection-point theorem, the refinement of which by Fjørtoft (1950) was later revealed to be one of the two conditions derived using Arnol'd's method (Arnol'd 1966). The existence of two distinct sufficient conditions was anticipated by Kelvin on energetic grounds (Thomson 1880), the same year as Rayleigh's inflection point theorem. Today they are called the Kelvin–Arnol'd first and second shear-stability theorems (KA-I and KA-II). KA-I corresponds to the Rayleigh–Fjørtoft condition, whereas KA-II is relevant in planetary physics.

A breach of sufficient conditions for stability does not necessarily imply instability. However, in applied fields, these conditions are sometimes treated as if they were sharp stability criteria, meaning that they accurately detect the stability boundary, contributing to potential confusion. The reason for this may be that the known necessary conditions for stability, i.e. conditions guaranteeing the existence of instability, require rather complex assumptions about the basic flow, $U(y)$. Tollmien (1935) was the first to find a class of basic flows where the KA-I condition becomes a sufficient and necessary condition for stability. Many assumptions are required of $U$ in his class, e.g. $U$ is symmetric and disappears on the wall. Howard (1964) also showed that KA-I becomes a sharp stability condition when the flow is considered in an unbounded domain and the basic flow is assumed to be in the class $\mathcal {H}$ he defined. Nevertheless the above results are only special cases and in general KA-I is not sharp. In a bounded domain, it is in fact possible to construct counterexamples that demonstrate that KA-I is not a necessary condition of stability (see Tollmien 1935; Drazin & Howard 1966; Drazin & Reid 1981).

The more general the basic flow considered, the more complex the argument required to derive necessary conditions for stability. Using the Nyquist method, Rosenbluth & Simon (1964) derived a sharp stability condition using an integral quantity that can be calculated from the basic flow. The assumption imposed by them on $U(y)$ is that it is monotonic and has only one inflection point in the domain. The latter assumption was removed by Balmforth & Morrison (1999), who also used the Nyquist method. Hirota, Morrison & Hattori (2014) used the variational method under almost the same assumptions and demonstrated that the definiteness of the quadratic form introduced by Barston (1991) is a necessary and sufficient condition for stability. In order to derive those conditions, it is critical to assume the monotonicity of $U(y)$. The special nature of this type of basic flow is that all neutral modes must possess only one critical layer, where the phase speed of the mode matches $U(y)$, and, moreover, the position of this layer is fixed at one of the inflection points.

Demonstrating the existence of a neutral mode first and then perturbing the wavenumber to find unstable modes has been a common approach since Tollmien (1935) and Howard (1964). The assumptions employed by Howard (1964) were aimed at enabling the use of Sturm–Liouville theory for neutral modes (see also Morse & Feshbach 1953), allowing the assertion of the existence of neutral modes for non-monotonic $U$. For the bounded flows, Tung (1981) investigated how much the assumptions about $U$ can be relaxed to demonstrate the presence of unstable modes around neutral modes. Although his analysis had some flaws, later, Lin (2003) independently completed a mathematically rigorous theory.

One notable example where the non-monotonicity of the basic velocity field becomes crucial is the alternating jet streams of Jupiter and Saturn. Their evolution is governed by the behaviour of Rossby waves, which are analogous to drift waves in plasma physics and arise when there is a gradient in the potential vorticity (PV), a conserved fluid-dynamical quantity formed by the combination of conservation of mass, momentum and thermal energy (Vallis 2017). Following standard practice, this article works within the quasi-geostrophic framework, which admits Rossby waves but filters out sound waves and inertia-gravity waves. The non-rotating, non-stratified KA-I result by Rayleigh and Fjørtoft was followed by extensions to rotating, non-stratified flow (Kuo 1949) and to rotating, stratified flow (Charney & Stern 1962). The latter result, widely known as the Charney–Stern stability criterion, established that a sufficient condition for stability is the absence of a point where the PV gradient changes sign, which we refer to as PV extremum. Note that in geophysics, a PV extremum is sometimes referred to as a critical latitude because the observed phase speed of the mode often matches $U(y)$ there. In this paper, however, a clear distinction is made between the two.

In terms of the basic state streamfunction, $\varPsi$, defined by $U=-{\rm d} \varPsi /{{\rm d}y}$, and the associated PV, $Q$, the KA-I and KA-II conditions can be expressed as $-{\rm d} Q/{\rm d}\varPsi \leq 0$ and $0\leq -{\rm d} Q/{\rm d}\varPsi \leq l^{-2}$, respectively, when perturbations have the largest length scale $l$ (McIntyre & Shepherd 1987). Non-dimensionalising those conditions for the planetary atmosphere problem uncovers the role played by the Rossby Mach number defined by

(1.1)\begin{equation} M=\frac{\kappa_0^2(U-\alpha)}{Q'}, \end{equation}

where $\kappa _0$ is a constant determined by the length scale of the eddies and $\alpha$ is a suitable Galilean shift. Dowling (2014, 2020) introduced $M$ based on physical intuition applied to the net propagation of the fastest Rossby waves relative to the flow, which are the longest waves. In brief, KA-I and KA-II can be interpreted as establishing that unstable flows must become ‘subsonic’ (i.e. $0< M<1$) somewhere, but with respect to Rossby waves rather than sound waves. Getting KA-I and KA-II to concatenate together via $M^{-1}\leq 1$ is a compact way to look at the stability condition. The above simple stability condition can only be applied for basic flows belonging to a certain class where the sign of $M^{-1}(y)$ does not change for all $y$. This class, in fact, closely resembles that defined by Howard (1964) and Lin (2003), and we base our analysis on this ground.

Dowling (1993) discovered that Jupiter's atmosphere has $U$ and $Q'$ profiles strongly correlated. This analysis, based on Voyager observations of the cloud-top vorticity field (Dowling & Ingersoll 1989), was expanded for Jupiter and Saturn by Read et al. (2006), Read et al. (2009b) and Read, Dowling & Schubert (2009a), as illustrated in figure 1. The correlation implies that the $M^{-1}(y)$ profile is almost constant. If the profile is neutrally stable and the KA-II condition is sharp, this constant must be unity. However, as already noted, the stability condition based on the reciprocal Rossby Mach number is only a sufficient condition for stability, and the significance of it regarding necessary conditions for stability remains not well understood. Moreover, the observed $M^{-1}(y)$ profiles are not precisely constant, as can be seen from the fact that the correlation in figure 1 is not perfect. The neutrally stable hypothesis for Jupiter's jets has been studied in the KA-II context for more than 30 years, as summarised in the most recent review article by Read (2024), but a conclusive resolution has yet to be reached.

Figure 1. Example Jupiter profiles with respect to latitude of zonal wind, $U$ (solid line), and PV gradient, $Q' = {\rm d} Q/{{\rm d}y} = Q_y$ (dashed line), from the southern hemisphere, at 11 hPa pressure. A strong correlation between $U$ and $Q_y$ is evident. Here, ${{\rm d}y} = R \,{\rm d}\phi$, where $y$ is latitude in metres, $\phi$ is latitude in radians (planetographic) and $R(\phi )$ is the local radius of curvature (see (5) in Dowling & Ingersoll 1989). These are a sample of the profiles produced by Read et al. (2006) using the Cassini Composite Infrared Spectrometer (CIRS). Similar correlations to this figure are found at other pressure levels.

Particularly intriguing from a planetary atmospheric physics perspective is thus when KA-II gives sharp stability boundaries. There are interesting numerical results in this regard: Stamp & Dowling (1993) set up a sinusoidal model basic flow such that $M$ becomes a constant, and numerically found that the instability disappears at the KA-II stability boundary, $M=1$. Motivated by these empirical results, our mathematical goals are twofold: (i) to derive simple conditions that guarantee the presence of instability, and (ii) to determine under what conditions KA-II achieves sharp stability boundaries. We also attempt to emphasise the simplicity of the theoretical results, because the pursuit of sharpness of the conditions often makes them difficult to use. For example, from a practical standpoint, the condition proposed by Balmforth & Morrison (1999) is not easily applicable due to the requirement of solving a Fredholm integral equation. In addition, to demonstrate the existence of instability using the quadratic form in Hirota et al. (2014) or the Rayleigh quotient in Lin (2003), it is necessary to find a convenient test function, but no useful recipe is known for non-monotonic basic flows in a finite domain.

The rest of the article is organised as follows. In § 2, the quasi-geostrophic equation is linearised around the basic flow to yield an eigenvalue problem. Our sufficient condition for instability, in the form of a hurdle of the $M^{-1}(y)$ profile, is stated in § 3, together with the KA stability theorem. Basic states are classified in the same section to clarify when the stability conditions can be used. One of the classes is identified as having Jupiter-style PV extrema; it yields an almost sharp stability condition that aligns with observations of the jets on Jupiter and Saturn. In § 4, the stability criteria are illustrated with model basic flow profiles. Section 5 studies the implications of the theoretical results obtained in the previous two sections to planetary physics problems. Section 6 contains mathematical proofs of the article's theorems. Our strategy is to extend the Rayleigh quotient method, developed by Howard (1964), Tung (1981) and Lin (2003), to the quasi-geostrophic system and then utilise it to check the parameter dependence of eigenvalues. Section 7 concludes with a summary and discussion of the hurdle theorem concept.

2. Formulation of the problem

For readers not familiar with geophysics, we introduce some basic terminology before presenting our model. Our starting point is quasi-geostrophic conservation of PV on a beta-plane (Vallis 2017), ${{\rm D} Q}/{{\rm D} t}=Q_t-\varPsi _y Q_x+\varPsi _x Q_y=0$, where $\varPsi$ is the streamfunction for the predominantly horizontal flow and the quasi-geostrophic PV is written as $Q=f_0+\beta y+\varPsi _{xx}+\varPsi _{yy}+(\,{f_0^2}/{\rho }) (({\rho }/{N^2})\varPsi _z)_z$. Here, $t$ is time and $x$, $y$ and $z$ are spatial coordinates in the zonal (longitudinal), meridional (latitudinal) and vertical directions, respectively; when these appear as a subscript, the meaning is partial differentiation. The zonal and meridional wind components can be found as $U = -\varPsi _y$ and $V = \varPsi _x$, respectively. In the quasi-geostrophic framework, the static density and squared buoyancy frequency, $\rho$ and $N^2$, are given functions of $z$.

The quasi-geostrophic equation can be derived by first applying the shallow layer approximation to the Navier–Stokes equations and then taking the limit of strong stratification and rapid rotation. Conversely, the KA-I stability condition can be extended from the quasi-geostrophic framework to the primitive shallow-water framework (Ripa 1983). In the primitive context, note that admission of buoyancy (gravity) waves in addition to Rossby waves complicates the shear stability question; however, the quasi-geostrophic results continue to be useful (e.g. Stamp & Dowling 1993, figures 3 and 5).

In the above formulation, the terms $f_0+\beta y$ are the first two terms of the Taylor series of the planetary vorticity, $f$, which is also called the Coriolis parameter. The terms $\varPsi _{xx}+\varPsi _{yy}$ are the relative vorticity. The last term involving two vertical derivatives is the stretching vorticity. If separation of variables can be used in $z$, the stretching vorticity may be simplified via the vertical eigenvalue problem

(2.1)\begin{equation} \frac{f_0^2}{\rho}\left(\frac{\rho}{N^2}\varPhi_z\right)_z={-}L_d^{{-}2}\varPhi, \end{equation}

with suitable boundary conditions, in which case there is a different Rossby deformation length, $L_d$, for each vertical mode. Usually it is the smallest $L_d^{-2}>0$ (the first baroclinic mode) that is of interest. The corresponding first Rossby deformation length (hereafter simply denoted as $L_d$) is the length scale that separates potential-energy dominated structures, i.e. large-scale pressure highs and lows maintained by the Coriolis effect, from kinetic-energy dominated structures, i.e. small-scale pressure anomalies that get flattened by gravity. The physical role played by $L_d$ is analogous to the Larmor radius in plasma flows (Hasegawa 1985).

2.1. The $1\frac {3}{4}$ model and the eigenvalue problem

Our formulation is based on multi-layer quasi-geostrophic systems (see § 5.3.2 of Vallis 2017). It is common practice in geophysical fluid dynamics to simplify the problem by studying a multi-layer shallow-water model in which constant-density layers are arranged in a (hydro)statically stable manner, with low density overlying high density. A two-layer model of this type, with a fully dynamic weather layer that overlies a layer containing a deep jet profile, $U_{deep}$, was first applied to Jupiter by Ingersoll & Cuong (1981). In this case, the weather layer PV is written as

(2.2)\begin{equation} Q=f_0+\beta y+\varPsi_{xx}+\varPsi_{yy}-L_d^{{-}2}(\varPsi-\varPsi_{deep}), \end{equation}

where $L_d$ is the first Rossby deformation length and $\varPsi _{deep}$ is the deep-layer streamfunction, also known as the dynamic topography.

The ratio of the depth of the deep layer to the weather layer is taken to be very large, so that $\varPsi _{deep}$ may be treated as steady (Majda & Wang 2005). We further assume for simplicity that the deep-layer circulation does not vary with longitude, $x$. Such variations can affect the phase-locking of long Rossby waves and thereby can play an indirect role, however the focus here is on unidirectional shear instability, which is relevant given the predominantly zonal nature of gas giant circulations. Under those assumptions, $\varPsi _{deep}$ and $U_{deep} = -(\varPsi _{deep})_y$ are functions of $y$ only. The weather layer corresponds to the first-baroclinic-mode structure of the atmosphere of a gas giant, whereas the barotropic structure of the gas-giant interior is modelled by $\varPsi _{deep}$. This two-layer configuration has traditionally been called the ‘$1\frac {1}{2}$ layer’ model, but was recently rebranded as the ‘$1\frac {3}{4}$ layer’ model when $U_{deep}$ has an alternating jet profile rather than being still ($U_{deep} = 0$), to emphasise the effect on the weather-layer dynamics that the corresponding undulations of $\varPsi _{deep}$ have via stretching vorticity (Flierl, Morrison & VilasurSwaminathan 2019).

In summary, the $1\frac {3}{4}$ layer model is written as conservation of quasi-geostrophic PV ${{\rm D} Q}/{{\rm D} t}=0$ for (2.2) which now yields one nonlinear equation in one unknown, $\varPsi (x,y,t)$. As noted previously, gas giants are observed to be predominantly zonal, hence for the linear stability analysis the basic state is assumed to be only a function of latitude, $y$. The $x$ dimension on a gas giant is periodic, hence the $x$ and $t$ dependencies are represented with a Fourier component by replacing $\varPsi (x,y,t)$ in (2.2) with $\varPsi (y)+\delta \psi (y) \exp [{\rm i}k(x-ct)]$, where $k$ and $c$ are the zonal wavenumber and phase speed, respectively, and $\delta$ sets the scale of the amplitude. The $y$ dimension is assumed to span a channel of width $L$ centred on the origin, $y \in (-L/2,L/2)$. When considering Jupiter's atmosphere, we assume that the channel is contained inside the northern or southern extratropical domain (i.e. poleward of the equatorial jet), where quasi-geostrophic theory holds. The meridional boundary conditions on the perturbations are Dirichlet, $\psi (-L/2)= 0$ and $\psi (L/2)=0$, unless otherwise stated.

The problem is linearised by restricting to $|\delta | \ll 1$ and retaining only $O(\delta )$ terms. These standard assumptions yield a linear ordinary differential equation that governs the meridional structure of small-amplitude perturbations,

(2.3)\begin{equation} \psi''-(k^2+L_d^{{-}2})\psi +\frac{Q'}{U-c}\psi=0,\quad y \in \varOmega, \end{equation}

where $\varOmega =(-L/2,L/2)$ and a prime denotes ordinary differentiation with respect to $y$. Three length scales exist for this problem: the channel width, $L$, the deformation length, $L_d$, and the scale at which the basic flow varies, $L_U$. In what follows, dimensional expressions are retained when addressing physical phenomena, but in the mathematical case studies expressions are implicitly non-dimensionalised using the length scale $L_U$ (i.e. $L$ and $L_d$ are $L/L_U$ and $L_d/L_U$ in the dimensional form, respectively).

In the $1\frac {3}{4}$ layer model, $U(y)=-\varPsi '(y)$ and $Q'(y)$ are linked by

(2.4)\begin{equation} Q'(y)=\beta -U''(y)+L_d^{{-}2} (U(y)-U_{deep}(y)). \end{equation}

Note that in this model a Galilean transformation in $x$ is applied for both layers so that $U$, $U_{deep}$ are replaced by $U-\alpha$, $U_{deep}-\alpha$, respectively, implying that $Q'$ is unchanged (note that this is not the case in the $1\frac {1}{2}$ layer model, where $U_{deep}$ is related to bottom topography and, hence, not altered in the Galilean shift). The link (2.4) implies that two of the three basic flow profiles, $U(y)=-\varPsi '(y)$, $U_{deep}(y)$ and $Q'(y)$ may be specified independently, after which the third is fixed. This is important in planetary science, as observing $U$ and $Q'$ determines $U_{deep}$, providing new insights into the nature of deep jets (Dowling 1995b; Read et al. 2006, 2009b,a). Furthermore, as we shall see shortly, if there is a neutrality hypothesis that can potentially be used to constrain the reciprocal Rossby Mach number, to be defined in (3.3), then $U$ and $Q'$ can be related. This is useful as $U$ is easily observable, whereas precise measurements of $Q'$ require accurate temperature retrievals and, thus, pose a relatively greater challenge. Physically, the gas giant zonal winds are observed to be quite stable; hence, in planetary physics, determining neutral Rossby Mach number profiles has been a focus.

Given a wavenumber, $k\geq 0$, and deformation radius, $L_d> 0$, (2.3) and the boundary conditions constitute an eigenvalue problem for the complex phase speed, $c=c_r+ic_i$. For fixed $L_d$, if there is a $k$ that yields $c_i> 0$, the flow is unstable. It is easy to confirm that the complex conjugate of the eigenvalue is also an eigenvalue; however, when $|c_i|$ is small, only the mode with $c_i>0$ provides a good approximation to the viscous problem.

In the wider context, it is important to note that our mathematical analysis holds for the quasi-geostrophic equation linearised around a broad range of $U(y)$ and $Q'(y)$ profiles. In the case of rotation with a non-zero planetary vorticity gradient, $\beta \neq 0$, the necessary and sufficient conditions for stability have also been discussed in the context of the analogous problem in plasma physics (e.g. Numata, Ball & Dewar 2007; Zhu, Zhou & Dodin 2018). In the case of no rotation and no stratification, $\beta = 0$ and $L_d \rightarrow \infty$, the PV gradient (2.4) simplifies to $Q' = -U''$ and (2.3) reduces to the original Rayleigh equation, in which case the PV extrema are inflection points of $U(y)$.

3. The stability conditions

The stability conditions for the eigenvalue problem can be succinctly expressed using the reciprocal of the Rossby Mach number, as we shall see shortly. In the subsequent two subsections, we define two classes of basic flows: class (i) and class (ii). This classification is important as the strength of the stability conditions varies between them. It is convenient to define the set of the PV extrema $Y_Q=\{y \in \varOmega \,|\, Q'(y)=0\}$ and the set of the zonal wind zeros in a suitable Galilean frame $Y_{U,\alpha }=\{y \in \varOmega \,|\, U(y)=\alpha \}$; the number of elements of the latter set depends on $\alpha$. Hereafter, we assume that $Q'(y)$, $U(y)$ are $C^1$. Therefore, unless $(U-c)$ vanishes somewhere, the eigenfunctions belong to the function space

(3.1)\begin{equation} C_0^2=\{\,f\in C^2(\bar{\varOmega})\,|\, f({-}L/2)=f(L/2)=0\},\end{equation}

where $\bar {\varOmega }=[-L/2,L/2]$. For neutral modes, the critical latitudes are expressed as $Y_{U,c}$.

3.1. The reciprocal Rossby Mach number

The reciprocal Rossby Mach number, $M^{-1}(y)$, emerges as being central to this work and there are various ways to motivate it. For example, in § 1 we briefly commented on the physical interpretation by Dowling (2014, 2020). The following motivation is based on mathematical consideration applied to the second-order ordinary differential equation (2.3). If (2.3) possesses a non-trivial solution that meets the boundary conditions, then in order for the solution to exhibit oscillatory behaviour at some point, the factor $Q'/(U-c)$, needs to be sufficiently larger than the smallest eigenvalue, $\kappa _0^2$, of the eigenvalue problem formulated by the remaining terms with $k=0$:

(3.2)\begin{equation} \varphi''-L_d^{{-}2}\varphi={-}\kappa^2 \varphi,\quad y \in \varOmega ,\end{equation}

with $\varphi (-L/2)=\varphi (L/2)=0$. This leads to the introduction of the reciprocal Rossby Mach number,

(3.3)\begin{equation} M_{\alpha}^{{-}1}(y)=\frac{1}{\kappa_0^2}\frac{Q'}{U-\alpha},\end{equation}

where the subscript alpha explicitly shows the dependence on a reference-frame shift of the zonal wind, $\alpha \in \mathbb {R}$. At this stage, $\alpha$ is arbitrary. However, as we will explain later, in the Jupiter-style basic flows, there is only one optimal choice of $\alpha$, allowing us to remove the subscript. Note from (2.4) that $Q'$ is invariant with respect to $\alpha$. Both $\kappa _0^2$ and its associated eigenfunction, $\varphi _0$, can be computed explicitly:

(3.4a,b)\begin{equation} \kappa_0^2=\frac{{\rm \pi}^2}{L^2}+L_d^{{-}2},\quad \varphi_0=\cos\left (\frac{\rm \pi}{L}y \right).\end{equation}

3.2. A sufficient condition for stability: class (i) basic flows

The extant sufficient conditions for stability can be expressed in terms of $M_{\alpha }^{-1}$ as follows.

Theorem 3.1 Suppose there exists $\alpha$ such that $M_{\alpha }^{-1}\leq 0$ for all $y\in \varOmega$, or that $0\leq M_{\alpha }^{-1} \leq 1$ for all $y\in \varOmega$. Then there is no unstable mode for any $k$.

Here, the former and latter conditions correspond to KA-I and KA-II, respectively. Mathematically, the stability conditions can be used in the equality cases, although they are often omitted in the physics literature (Dowling 2014). Our system is classified as a non-canonical Hamiltonian system, and the theorem can be shown by Arnol'd's method; if a Hamiltonian $H$ and a Poisson bracket $\{,\}$ satisfy $\delta H:=\{F,H\}=0$ for all functionals $F$ at a steady solution, then the solution is stable if $\delta ^2 H:=\{F,\{F,H\}\}$ is strictly positive or negative definite for all $F$. The method has a wide range of applications and can also derive nonlinear stability with respect to finite amplitude disturbances. However, if we restrict ourselves to the linear eigenvalue problem, Theorem 3.1 can be proved by a more straightforward approach, as shown in Appendix A.

If $Q'(y)$ does not change sign over $\varOmega$, i.e. there are no PV extrema, we can always choose a large enough or small enough $\alpha$ to satisfy KA-I, which recovers the Charney–Stern condition (KA-I is also known as a sufficient condition for the absence of over-reflections of Rossby waves; see Lindzen & Tung 1978). In other words, the interesting case occurs when there is at least one PV extremum. We assume the following properties for both class (i) and class (ii) basic flows:

1. (a) the PV gradient $Q'(y)$ changes sign somewhere in $\varOmega$;

2. (b) the zeros of $Q'(y)$, the PV extrema, are isolated;

3. (c) $Q''(y_l)\neq 0$ for all $y_l \in Y_Q=\{y \in \varOmega \,|\, Q'(y)=0\}$.

We say a basic flow belongs to class (i) if there exists an $\alpha \in \mathcal {R}$ such that the following additional condition is satisfied:

1. (d) $M_{\alpha }^{-1}(y)$ is continuous and one-signed in $\varOmega$.

Here $\mathcal {R}$ is the range of the zonal flow

(3.5)\begin{equation} \mathcal{R}=\left(\min_{y\in \bar{\varOmega}}U,\max_{y\in \bar{\varOmega}} U\right ). \end{equation}

In the definition, ‘a function is one-signed’ means that it is non-negative or non-positive for all $y \in \varOmega$. Of course, from Theorem 3.1, the non-positive cases are stable.

Class (i) occurs when the PV extrema and the zeros of $U-\alpha$ coincide, i.e. $Y_Q=Y_{U,\alpha }$. For example, this is the case when $Q'$ in figure 2(a) is paired with $U$ in figure 2(b). On the other hand, the basic flow is not class (i) when $U$ is replaced by the profile shown in figures 2(c) or 2(d).

Figure 2. Classification of basic flow profiles. For all these examples, consider the PV gradient $Q'(y)$ profile shown in (a). The two PV extrema are indicated by the horizontal red dashed lines. The $U(y)$ profile shown in (b) is class (i). The $U(y)$ profile shown in (c) is class (ii) but not class (i). The $U(y)$ profile shown in (d) is not class (ii).

The importance of class (i) can be seen from the fact that it is the case that Theorem 3.1 may be able to show stability, when there is a PV extremum. The KA-I and KA-II stability conditions can be combined; the flow is stable if

(3.6)\begin{equation} \max_{y\in \bar{\varOmega}} M^{{-}1}\leq 1. \end{equation}

Here the reciprocal Rossby Mach number is simply written as $M^{-1}$, as the choice of $\alpha$ is trivial (see Theorem 6.1 in § 6.1). Note if $M_{\alpha }^{-1}(y)$ changes sign (i.e. the basic flow is not class (i)), mathematically the condition (3.6) cannot guarantee stability.

Class (i) is the geophysically interesting case to which the ‘Jupiter-style’ shear flow belongs, assuming a suitable $U_{\rm {deep}}$ (Dowling 2020; Afanasyev & Dowling 2022). As mentioned in § 1, observations support that $M^{-1}$ is around unity in Jupiter's and Saturn's atmosphere (Dowling 1993; Read et al. 2006, 2009b,a). This configuration can be inferred from the most accurate available observational data, as depicted in figure 3. The link between class (i) and Jupiter's atmosphere is further discussed in § 5.

Figure 3. Comparison of Jupiter profiles of $U$ (solid) and $\kappa _0^{-2}Q_y$ (dashed), using the same data as in figure 1. Since $L/L_d \gg 1$ on Jupiter, it follows that $\kappa _0^2 \approx L_d^{-2}$, such that the similarity between the two profiles implies $M^{-1} \approx 1$ with $\alpha \approx 0$. This, and further evidence discussed in the text, support the identification of class (i) as the Jupiter-style class. The deformation length is set here to $L_d = c_{gw}/|\,f_0| = 1750\,{\rm km}$, using the gravity wave speed, $c_{gw} = 454\,{\rm m}\,{\rm s}^{-1}$, inferred from the Comet Shoemaker-Levy 9 impact analysis (Hammel et al. 1995) and the Coriolis parameter, $f_0$, evaluated at the G-fragment impact site, latitude $\phi _0 = -47.5^\circ$ (planetographic); the figure is centred on $\phi _0$.

3.3. A sufficient condition for instability: class (ii) basic flows

Basic flows that are not class (i) may also be physically important. For example, $M_{\alpha }^{-1}$ changes sign on a seasonal basis in Earth's atmosphere (Du, Dowling & Bradley 2015). Our main result, Theorem 3.2, which we call the hurdle theorem, can show the existence of instability for a wider range of basic flows than class (i).

We say a basic flow belongs to class (ii) if there exists an $\alpha \in \mathcal {R}$ satisfying the following conditions, called class (ii) conditions:

1. (a) $M_{\alpha }^{-1}(y)$ is continuous in $\varOmega$;

2. (b) $M_{\alpha }^{-1}(y_j)$ is one-signed for all $y_j \in Y_{U,\alpha }$.

Note that if a basic flow is class (i), then it is class (ii). However, for class (ii) basic flows, in general there can be more than one $\alpha$ that make $M_{\alpha }^{-1}$ continuous. In figure 2(a,c), for $\alpha =0$ the lower PV extremum cancels the singularity and $M_{\alpha }^{-1}$ changes sign at the upper PV extremum. Alternatively, an appropriately negative $\alpha$ may be chosen such that $U-\alpha$ vanishes at the upper PV extremum.

It is useful to define the sets $\{U(y)\,|\,y\in \varOmega _+\}$ and $\{U(y)\,|\,y\in \varOmega _-\}$, decomposing $\varOmega$ into the three parts, $\varOmega _+=\{y\in \varOmega \,|\, Q'(y)>0 \}$, $\varOmega _-=\{y\in \varOmega \,|\, Q'(y)<0 \}$ and $Y_Q$. Then if the set

(3.7)\begin{equation} \mathcal{D}=\mathcal{R}\backslash (\{U(y)\,|\,y\in \varOmega_+\}\cup \{U(y)\,|\,y\in \varOmega_-\}) \end{equation}

is non-empty, we may be able to select $\alpha$ from this set. Furthermore, in the vicinity of the chosen point, there should be no points belonging to

(3.8)\begin{equation} \{U(y)\,|\,y\in \varOmega_+\}\cap \{U(y)\,|\,y\in \varOmega_-\}\end{equation}

to satisfy the second class (ii) condition.

To visually check whether the basic flow is class (ii), first find $\varOmega _+$ and $\varOmega _-$ using $Q'(y)$ and then illustrate $\{U(y)\,|\,y\in \varOmega _+\}$ and $\{U(y)\,|\,y\in \varOmega _-\}$. For example, with $U(y)$ shown in figure 4(b,c), the number of elements in $\mathcal {D}$ is one and two, respectively. The condition that the point at which $\{U(y)\,|\,y\in \varOmega _+\}$ and $\{U(y)\,|\,y\in \varOmega _-\}$ are separated in the range of $U(y)$ is not particularly restrictive, such that a wide array of basic flows belongs to class (ii). In particular, if $U(y)$ is monotonic and there is at least one PV extremum, $\mathcal {D}$ should be non-empty. On the other hand, for the basic flow $U(y)$ shown in figure 2(d), $\mathcal {D}$ is empty, and in fact the basic flow cannot be class (ii).

Figure 4. The PV gradient profile $Q'(y)$ defines the sets $\varOmega _+=\{y\in \varOmega \,|\, Q'(y)>0\}$, $\varOmega _-=\{y\in \varOmega \,|\, Q'(y)<0 \}$. For the zonal wind profiles $U(y)$ shown in (b,c), the set $\mathcal {D}=\mathcal {R}\backslash (\{U(y)\,|\,y\in \varOmega _+\}\cup \{U(y)\,|\,y\in \varOmega _-\})$, has two and one points, respectively (indicated by orange bullets and lines). Furthermore, in the vicinity of these points, there are no points belonging to the set (3.8). Therefore, both cases are class (ii).

Our main result is the following ‘hurdle theorem’, which provides a sufficient condition for instability (and the contrapositive necessary condition for stability). This theorem is proved in § 6.

Theorem 3.2 Suppose the basic flow is class (ii). Fix an $\alpha$ so that the class (ii) conditions are satisfied. The flow is unstable if there is an interval $[y_1,y_2] \subset \varOmega$ over which

(3.9)\begin{equation} h\equiv \frac{\dfrac{{\rm \pi}^2}{(y_2-y_1)^2}+L_d^{{-}2}}{\dfrac{{\rm \pi}^2}{L^2}+L_d^{{-}2}}< M_{\alpha}^{{-}1}\end{equation}

is satisfied.

In short, if the reciprocal Rossby Mach number profile surpasses the hurdle $h$, then the flow is unstable, as illustrated in figure 5. As can be seen from the proof in § 6, this result is actually independent of the conditions applied at the boundaries (though it depends on $L$).

Figure 5. Example of an unstable $M^{-1}_{\alpha }(y)$ profile. Observe that the blue curve overcomes the hurdle $h$.

Note that the height of the hurdle, $h$, is always higher than unity, and depends on the width of the hurdle. However, if $L_d^{-2} \gg {\rm \pi}^2/(y_2-y_1)^2$, the hurdle height is almost 1. Thus, in this limit the condition for the existence of instability asymptotes to

(3.10)\begin{equation} \max_{y\in \bar{\varOmega}} M_{\alpha}^{{-}1}>1.\end{equation}

This is the contrapositive of (3.6), implying that the KA-II stability condition is almost sharp if the basic flow is class (i) and $L_d$ is small (recall that for class (ii) the stability condition (3.6) is not valid).

4. Analysis of model profiles

In this section, we stress-test the hurdle theorem (Theorem 3.2) with model basic flow profiles. Stamp & Dowling (1993) used the assumption $Q'=aL_d^{-2}U$, such that $M^{-1}_0$ becomes a constant $a$, together with a sinusoidal zonal wind profile $U =\cos (2{\rm \pi} y/L)$, and applied periodicity in both meridional and zonal directions. The assertion of Theorem 3.1 remains the same with Neumann or periodic boundary conditions for (3.2), in which case (3.4a,b) is replaced by

(4.1a,b)\begin{equation} \kappa_0^2=L_d^{{-}2},\quad \varphi_0=1. \end{equation}

Their basic flow is clearly class (i), and the numerical results are consistent with this theoretical result. Moreover, their numerical computations suggest that whenever $M_0^{-1}\!>\!1$, unstable modes exist for some $k$, implying that KA-II may give a sharp stability boundary. Theorem 3.2 validates this numerical result, because when $M_0^{-1}$ is a constant, we are able to employ the full-width hurdle with unit height, $h=1$.

In general, the reciprocal Rossby Mach number will not be a constant but rather a profile that may be ‘supersonic’ in some regions, $M^{-1}_\alpha (y) < 1$, and ‘subsonic’ in others, $M^{-1}_\alpha (y) > 1$. To explore the ramifications of this generality to shear instability, we next develop a model with a two-parameter, variable $M^{-1}_\alpha (y)$ profile. We use this model to numerically explore first class (i) profiles in § 4.1, and then class (ii) profiles that are not class (i) in § 4.2. Section 5 delves into the connection between the model and the stability of Jupiter's alternating jet streams. As commented at the end of § 2, our problem covers the classical Rayleigh equation problem. In Appendix C we analyse one more model flow to clarify where our analysis stands in the long history of Rayleigh-equation research.

4.1. Class (i)

Consider the quasi-geostrophic problem (2.3). In setting up a model basic flow, we assume that the PV gradient profile has the form

(4.2)\begin{equation} Q'(y)=\kappa^2_0\{a+(b-a){\rm sech}^2 (y)\}U(y), \end{equation}

where $a,b> 0$ are two free parameters. The advantage of using this designed basic flow is that, regardless of $U(y)$, the reciprocal Rossby Mach number with the choice $\alpha =0$ becomes

(4.3)\begin{equation} M_{0}^{{-}1}(y) =a+(b-a){\rm sech}^2 (y). \end{equation}

The case $a=b$ corresponds to the PV profile considered in Stamp & Dowling (1993). When $b > a$, the $M_0^{-1}$ profile has a hump of height $b$ centred at $y=0$, and likewise a dip when $b < a$.

If $U(y)$ has a zero (i.e. $0 \in \mathcal {R}$), the existence of a PV extremum is guaranteed. Then, since both $a$ and $b$ are positive, the function (4.3) is one-signed for all $y$, implying that the model basic flow belongs to class (i), i.e. ‘Jupiter-style’. The two different zonal wind profiles

(4.4)$$\begin{gather} U(y)=\sin(2{\rm \pi} y), \end{gather}$$

(4.5)$$\begin{gather}U(y)=\tanh(y), \end{gather}$$

are considered. Here, lengths are implicitly scaled by $L_U$. The profile (4.4) may be regarded as a model for Jupiter's alternating jets (with $L_U$ being the jet peak-to-peak length scale), whereas (4.5) is a simpler case where there is only one PV extremum.

For class (i) the phase speed of the neutral modes is uniquely determined (§ 6, Theorem 6.1), and for our basic flows it is $0$. Thus, setting $c=0$ in (2.3), the equation satisfied by the neutral solutions can be found as

(4.6)\begin{equation} \psi''-(k^2+L_d^{{-}2}-\kappa_{0}^2M^{{-}1})\psi = 0, \end{equation}

simplifying the notation as $M^{-1}=M_0^{-1}$. This equation and the Dirichlet boundary condition constitute an eigenvalue problem with $\lambda =-k^2$ as the eigenvalue. Considering that the model flows are meant to emulate the atmosphere of Jupiter, contemplating the case of large $L$ is natural.

When $L$ is large, it is possible to study (4.6) analytically. To see the behaviour of a typical neutral mode, let us choose the parameter $L_d^{-2}(b-a) = 2$ that makes the analysis simple. First, we note that $\psi =\text {sech}y$ approximately satisfies (4.6) when $k=\sqrt {1+L_d^{-2}(a-1)}$, because $\kappa _0^2\approx L_d^{-2}$. This mode decays exponentially for large $|y|$ as shown by the red curve in figure 6(a). Such localised modes do not interact with the boundary, and they are robust to changes in $L$. However, there are also oscillatory modes that do not show evanescent behaviour near the boundaries; see the green curves in figure 6(a). Those curves can also be found theoretically. Let $\omega =\sqrt {L_d^{-2}(a-1)-k^2}$ be real. Then $\psi =\omega \cos (\omega y)-(\tanh y)\sin (\omega y)$ and $\psi =\omega \sin (\omega y)+(\tanh y)\cos (\omega y)$ satisfy (4.6). If $L$ is large, the former and latter solutions approximately satisfy the boundary conditions when $\tan (\omega L/2)=\omega$ and $\tan (\omega L/2)=-\omega ^{-1}$, respectively. There are many $\omega$ that meet these requirements, and whenever $k=\sqrt {L_d^{-2}(a-1)-\omega ^2}$ is real, (4.6) has a neutral mode that oscillates near the boundaries. The number of allowed values of $k$ increases with increasing $L$, as shown in figure 6(b), and in the limit of $L\rightarrow \infty$ the eigenvalues of the oscillatory modes form a continuous spectrum occupying $k \in (0,L_d^{-1}\sqrt {a-1})$.

Figure 6. (a) The eigenmodes of (4.6) for $L=16$, $L_d=1$, $a=1.5$ and $b=3.5$. There are four neutral modes for this parameter set. The solid curves are the theoretical results assuming large $L$ ($k\approx 1.225, 0.6707, 0.5494, 0.2476$). The symbols are the numerical results ($k\approx 1.269, 0.7178, 0.6016, 0.3564$). The red curve and the filled circles is the localised mode. The green curves and the open circles are the oscillatory modes. (b) The neutral $k$ for $L_d=1$, $a=1.5$ and $b=3.5$. The solid curves are the theoretical results, and the larger $L$, the better the agreement with the numerical results, shown by the symbols. In the limit $L\rightarrow \infty$, the oscillatory modes (green curves) form a continuous spectrum in the interval $k\in (0,\sqrt {0.5})$. There is only one localised mode (red curve) at $k=\sqrt {1.5}$.

For $L_d^{-2}(b-a) \neq 2$, the theoretical analysis is a bit more complicated (Appendix B). However, the final results are neat and clean as summarised in figure 7(a), which shows the parameter plane with $b-1$ and $a-1$ as abscissa and ordinate. Clearly the third quadrant must be stable as labelled because of the KA stability condition (3.6). Moreover, if $L$ is large, the first and second quadrants are unstable; the easiest way to see this is to note that the right-hand side of (6.18) is almost $a$, because the integral is largely unaffected by what happens in the hump or dip. The analysis just below (B2) also shows that even if $a$ is slightly above 1, many oscillatory modes will appear when $L$ is large. The stability of the flow is in fact non-trivial only in the fourth quadrant of figure 7(a), but somewhat surprisingly, the stability in this region can be found analytically (see (B4)). In summary, the theoretical neutral curve for $L\gg 1$ can be described as

(4.7)\begin{equation} a =\begin{cases} 1, & \text{if}\ b \in (0,1) ,\\ 1-L_d^{{-}2}(b-1)^2, & \text{if}\ b \geq 1. \end{cases} \end{equation}

This neutral curve delimits the stable and unstable parameter regions, as shown in figure 7(a). In the fourth quadrant, only localised modes appear, and this is the reason why the neutral curve is so simple. Applying the shooting or Chebyshev collocation method to (4.6), we can calculate the neutral curve for finite $L$. The computational results indeed approach the theoretical result (4.7) as $L$ is increased, as shown in figure 7(b).

Figure 7. (a) The theoretical neutral curve (4.7) for $L_d=1$. (b) Comparison of the theoretical neutral curve (green) and the numerical results (black) for $L_d=1$. The three numerical neutral curves are computed using $L=10,12,16$. (c) Comparison of the numerical neutral curve (black) and the unstable parameter region by Theorem 3.2 (red shaded). Here $L=16$, $L_d=1$. The $b$–$a$ plane is used to make it easier to see the relationship with figure 10. (d) The theoretical curve (4.7) for $L_d^{-1}=1,2,4,8$. The blue shaded region is stable from Theorem 3.1.

One of the novel features of Theorem 3.2 is that, without resorting to eigenvalue computations, applying a simple hurdle yields a useful estimate of the behaviour of neutral curves. In figure 7(c), the unstable region is depicted, which can be identified by setting the various hurdles for $M^{-1}(y)$. The numerically calculated neutral curve should be sandwiched between the KA stable region and the unstable region identified by the hurdle theory. The latter region has a piecewise smooth boundary because the behaviour of hurdles in each quadrant is different. The entire first quadrant (with respect to the pivot point (1,1)) is unstable, as can be found by considering the full width hurdle of height unity. In the second quadrant the hurdle giving the best results sits between one of the boundaries and the dip. The third quadrant is KA stable. In the fourth quadrant, when $b$ is large enough the hump will be able to hurdle over. For this to be seen, the value of $b$ must be at least larger than 6 for $L_d=1$, meaning that the hurdle estimation does not give sharp results. However, as noted just below (3.10), the smaller $L_d$ becomes, the more accurate is the hurdle at detecting instability. In line with this expectation, the theoretical neutral curve approaches the KA stability boundary with decreasing $L_d$ (figure 7d).

The neutral modes do not depend on $U(y)$, but the unstable modes stemming from them do. Figure 8 shows the eigenvalue $c$ of (2.3) calculated for the same parameters as in figure 6(a). For the oscillatory zonal flow (4.4), all instabilities persist down to $k=0$ (figure 8a). However, for the monotonic zonal flow (4.5), the first and second neutral modes, as well as the third and fourth neutral modes, are connected as seen in figure 8(b), resulting in a qualitatively completely different diagram. For both cases, the growth rate $kc_i$ is well below the analytically derived upper bound shown by the dotted curve (see the supplementary material where Pedlosky's bound is tightened using the approach by Deguchi (2021)). The streamfunction at the maximum growth rate is indicated by arrows. The fastest growing mode inherits the properties of the localised neutral mode and forms a strong vortex in the centre of the region. The second and subsequent modes spread over the whole domain, like the oscillatory neutral modes.

Figure 8. The black crosses are numerical solutions of the eigenvalue problem (2.3) with the PV gradient profile (4.2) and the parameters $L=16$, $L_d=1$, $a=1.5$ and $b=3.5$: (a) $U=\sin (2{\rm \pi} y)$; (b) $U=\tanh (y)$. The magenta squares are the neutral modes computed in figure 6(a). The dotted curves indicate the analytic upper bound of $c_i$ shown in the supplementary material available at https://doi.org/10.1017/jfm.2024.728. The panels on the right depict contour plots of the streamfunction $\psi (y)\,{\rm e}^{{\rm i}kx}+\psi ^*(y)\,{\rm e}^{-{\rm i}kx}$ at the wavenumber that maximises the growth rate $kc_i$. The sign of the contours is distinguished by red and blue.

Of particular interest from a planetary physics perspective is the fourth quadrant of figure 7(a). Figure 9 shows the growth rates computed for the two different parameter sets, $(L,L_d,a,b)=(16,1,0.5,2)$ and $(L,L_d,a,b)=(16, 1/8, 0.5, 1.2)$. The latter setting is somewhat close to the situation found in Jupiter's atmosphere, as we shall see in § 5. In the chosen base flow $U=\sin (2{\rm \pi} y)$, there exist many critical latitudes, and in § 6 it is shown that they must coincides with the PV extrema for class (i). For both eigenfunctions shown in figure 9, one can observe that disturbances are localised near the centre of the domain, i.e. where the critical latitudes are subsonic ($M^{-1}>1$). Since $a=0.5$, other critical latitudes are supersonic ($M^{-1}<1$). The occurrence of vortices at subsonic latitudes is not a phenomenon specific to the current model. Mathematically, this expectation arises from the fact that, according to Sturm's oscillation theorem, neutral modes exhibit oscillatory behaviour only when $M^{-1}$ well exceeds unity. When $L_d$ is small, as soon as the hump of $M^{-1}$ exceeds 1, an unstable mode occurs, as noted at the end of § 3.3. This is why the eigenfunction for the $L_d=1/8$ case shown in figure 9 has smaller vortices.

Figure 9. Stability analysis in the fourth quadrant of figure 7(a), with $U=\sin (2{\rm \pi} y)$. Red plus symbols: $L=16$, $L_d=1$, $a=0.5$ and $b=2$. Green cross symbols: $L=16$, $L_d=1/8$, $a=0.5$ and $b=1.2$. In both cases, the most unstable eigenfunctions are indicated by arrows. The red and blue curves are contours of the streamfunction; see figure 8 caption.

4.2. Class (ii)

Interesting things happen when considering negative $b$ in (4.2). There is no longer any guarantee that the basic flow is class (i), as the sign of the $M^{-1}_0$ profile (4.3) may change. The neutral curve does depend on the choice of $U$, because the phase speed of neutral modes is not always zero. Here we mainly consider the second zonal wind profile $U(y)=\tanh (y)$, fixing $L=16, L_d=1$. When $a$ and $b$ are positive, the neutral curve is obtained as shown in figure 7(c). What happens if the neutral curve is extended to the region where $b$ is negative? For example, consider the situation when $a$ is smaller than 1 and $b$ is negative; the results look like figure 10. The emergence of the unstable region is due to the modes with non-zero phase speed, because clearly $M_0^{-1}$ will be everywhere smaller than 1.

Figure 10. The stability of the basic flow $U=\tanh (y)$ with the parameters $L=16$, $L_d=1$. The black curve is the numerically obtained neutral curve. The red shaded area is the unstable parameter region found by Theorem 3.2.

It is easy to see that the PV extrema are zeros of $U(y)$ and $\{a+(b-a){\rm sech}^2 (y)\}$. For example, consider the case $a=0.5$, $b=-0.5$. The zeros of $\{a+(b-a){\rm sech}^2 (y)\}$ are at $\pm 0.8814$, and thus the PV extrema are $Y_Q=\{-0.8814,0,0.8814\}$. One of these three points must coincide with the zero of the $U(y)-\alpha$, for $M^{-1}_{\alpha }$ to be continuous. Thus, there are neutral modes with $c_r=0$ and $c_r\approx \pm U(0.8814)\approx \pm 0.7071$. If $\alpha =0$ is chosen, the $M^{-1}_{\alpha }$ profile is everywhere less than unity as shown by the red curve in figure 11(a). Hence, the flow is stable for steady modes. However, when $\alpha =0.7071$ is used, $M^{-1}_{\alpha }$ exceeds 1 significantly, as indicated by the blue curve in figure 11(b).

Figure 11. (a) Reciprocal Mach number profile for $L_d=1$, $L=16$, $a=0.5$ and $b=-0.5$. The basic flow $U=\tanh (y)$ is used. (b) The numerically obtained neutral modes at the same parameters. They have the phase speed $c_r\approx 0.7071$.

Travelling-wave-type neutral modes with a phase speed of 0.7071 indeed appear for the parameter choice $a=0.5, b=-0.5$ (figure 11b). For $\alpha =0.7071$, class (ii) conditions are satisfied. This follows from the monotonicity of $U(y)$; see the comments below (3.8). As would be expected from these facts, unstable modes appear when $k$ is reduced from the neutral values (figure 12). The eigenfunctions with the largest growth rates are similar to the neutral mode, with vortices concentrated in regions where $M^{-1}_{\alpha }$ exceeds 1.

Figure 12. The same plots as figure 8 but for $L_d=1$, $L=16$, $a=0.5$ and $b=-0.5$. The basic flow $U=\tanh (y)$ is used. The neutral solutions, indicated by the magenta squares, are identical to those shown in figure 11(b). The red and blue curves are contours of the streamfunction; see figure 8 caption.

The fact that the neutral curve for $b<0$ depends on $U$ suggests that the situation is much more complicated when the basic flow is not of class (i). For example, if the basic flow $U(y)=\sin (2{\rm \pi} y)$ is used, the only reference shift $\alpha$ that would make $M^{-1}_{\alpha }$ continuous is 0, when $a>0, b<0$, and $L$ is large. However, this does not mean that considering only steady modes is sufficient. The reason for this is that when $\{U(y)|y\in \varOmega _- \}\cap \{U(y)\,|\,y\in \varOmega _+ \}$ is non-empty, there may be a singular neutral mode (see the remark below (6.4)). Such singular modes are beyond the scope of this paper and, in fact, there is no need to consider them in the planetary context as long as we assume that Jupiter's atmosphere is of class (i).

5. Implications of the model results to planetary atmospheres

In this section, we summarise implications of our model analysis in geophysics problems, including deep jets on Jupiter and Saturn and the inference of Saturn's rotation period (which is otherwise elusive because its magnetic field is not tilted). We then motivate non-dimensional, idealised cases and show the deep jets associated with neutral stability are reasonable, even with slightly varying reciprocal Rossby Mach number. We end the section with a few comments about weakly unstable scenarios for the neighbourhoods of prominent features such as Jupiter's Great Red Spot and Saturn's Polar Hexagon and Ribbon.

Before delving into the analysis of the model introduced in the last section, we emphasise that our purpose here is not to faithfully reproduce the details of Jupiter's atmospheric phenomena, but rather to elucidate key physical mechanisms. It is understood that the beta-plane $1\frac {3}{4}$ quasi-geostrophic model is derived from a series of simplifications, for example, the actual atmospheres of Jupiter and Saturn have continuous vertical structure. In addition, the value of $L_d$ may vary with latitude and we do not have good estimates in hand outside the vicinity of the latitude range shown in figure 3. The assumption that Jupiter's atmosphere falls under class (i) is inferred from observations, with the understanding that the actual planet is dynamically active, with long periods of stability punctuated by episodic storm outbreaks. Remote-sensing observations are accompanied by noise, and when substituting the figure 1 data into (3.3), it is evident that, regardless of how we choose the zonal wind shift, $M^{-1}_{\alpha }(y)$ cannot be made into a continuous function.

Nevertheless, there is a consensus among multiple independent research groups that the observed $Q'(y)$ and $U(y)$ tend to have the same sign on Jupiter and Saturn (Dowling 1993; Read et al. 2006, 2009b; Marcus & Shetty 2011), and when this property is ideal, the basic flow falls under class (i). As a minimum check to test that the link between the stability theory and the correlation seen in figure 3 is robust with respect to noise, the following numerical experiment was performed. Given the base flows, regardless of their class, the value of $L_d$ that makes the flow configuration neutral can be computed by the eigenvalue problem (2.3). Therefore, if $\kappa _0$ is computed from that $L_d$, we can compare $\kappa _0^{-2} Q'$ and $U$ to check their correlation. We spline interpolated $U$ and $Q_y$ in the latitude range shown in figure 3, corresponding to $L \approx 10\,860\,{\rm km}$. The numerical eigenvalue problem (2.3) yields $\kappa _0^{-1}\approx 1700\,{\rm km}$, which is reasonably close to the value $\kappa _0^{-1}= 1750\,{\rm km}$ used in figure 3. Moreover, the neutral wave has a relatively small phase speed $c_r\approx 6\,{\rm m}\,{\rm s}^{-1}$. As a consequence, plotting $(U - c_r)$ and $\kappa _0^{-2}Q'$ with the computed $\kappa _0^{-1}, c_r$ gives the same level of correlation as in figure 3. This result is robust with respect to the artificially set boundary conditions in the computation.

On Jupiter, the jet wavelength is $L_U \approx 2{\rm \pi} L_d$, as can be seen from figure 1 using the fact that $1^\circ$ latitude is approximately $1200\,{\rm km}$; a similar relationship holds for Saturn. Based on this observation, we choose $L_d =1/2{\rm \pi}$ in the non-dimensional model, employing the profiles (4.2) and (4.4). The black curve in figure 13 represents the numerically obtained neutral curve; as expected, it lies between the stable and unstable regions obtained by Theorems 3.1 and 3.2. The hurdle theorem result now better approximates the neutral curve than figure 7(c). In addition, figure 13 clearly demonstrates that even a slight deviation from the KA-II stability boundary may result in flow instability. As discussed in § 1, the introduction of the Rossby Mach number was motivated by the physical interpretation of KA-II. Our discovery that KA-II is relatively sharp under Jupiter's atmospheric conditions reinforces the validity of that physical mechanism.

Figure 13. The black curve is the neutral curve for $L_d =1/2{\rm \pi}$, $L=16$. The basic flow profiles (4.2) and (4.4) are used. The blue shaded region is stable from Theorem 3.1. The red shaded region is unstable from Theorem 3.2.

Stabilising alternating jets with PV extrema in a gas-giant weather layer, or in an analogous $1\frac {3}{4}$ layer system, via KA-II requires that there be alternating jets in the deep layer. Such deep jets must in turn be stabilised by some physical process other than KA-II. The key there appears to be that the deep jets operate in a quite different geometry: that of a rapidly rotating deep sphere instead of a shallow spherical shell, as investigated by Ingersoll & Pollard (1982). Regardless of the physics behind the synchronisation between the weather and deep jets, it is a fortunate circumstance for geophysicists because in practice, to detect stable PV extrema in a weather layer is to infer deep jets, which is how Jupiter's deep jets were first discovered, decades before the Juno gravity mission confirmed their existence (Dowling & Ingersoll 1989; Dowling 1993, 1995b). Our take is that KA-II provides the weather layer jets with the flexibility to adjust to alternating deep jets, and this represents a meaningful step forward in understanding the overall stability of the atmosphere-interior system.

Assuming that Jupiter's atmosphere favours neutral states, the deep layer profile $U_{deep}(y)$ can be calculated from (2.4). In the previous studies, zonal-wind pairs (weather layer and deep layer) appropriate for the $1\frac {3}{4}$ layer model applied to Jupiter were calculated from Voyager winds and vorticity data and plotted in Dowling & Ingersoll (1989) and Dowling (1995b, 2020). The deep-layer westward jets tend to be similar to the cloud-top westward jets, whereas the deep-layer eastward jets tend to be stronger than the weather-layer eastward jets by about 50 % (Dowling 1995b). For our model where Jupiter's weather-layer profile is idealised as a sinusoid, $U(y) = \sin (2{\rm \pi} y)$, if we take $M^{-1} = 1$ (i.e. $a=b=1$) with $\kappa _0^{-1} \approx L_d$, (2.4) yields the dimensional deep jet profile $U_{deep}(y) = U(y)+\beta L_d^2/U_0$, such that the deep-layer profile is the same sinusoid but with a positive shift. The size of this shift can be estimated for Jupiter assuming $\beta \approx 3.5\times 10^{-12}\,{\rm m}^{-1}\,{\rm s}^{-1}$, $L_d \approx 1750\,{\rm km}$ and a stratospheric wind amplitude, $U_0 \approx 25\,{\rm m}\,{\rm s}^{-1}$, as in figure 3, which yields $\beta L_d^2/U_0 \approx 10.7/25 \approx 0.4$. Alternatively, the $1\frac {3}{4}$ layer model applied to Jupiter's Great Red Spot, which is centred closer to the equator at latitude $-23^\circ$, would typically use $L_d \approx 2000\,{\rm km}$ and a tropospheric wind amplitude, $U_0 \approx 50\,{\rm m}\,{\rm s}^{-1}$, which yields $\beta L_d^2/U_0 \approx 14/50 \approx 0.3$; both are consistent with Voyager results (Dowling 1995a). The non-dimensionalised $U_{deep}(y)$ with $\beta L_d^2/U_0 = 0.3$ is depicted in figure 14(a) by the blue curve. We can also calculate $M^{-1}$ using other points from the neutral curve, for example, the neutral point $(a,b) = (0.6,1.1)$ yields the $U_{deep}(y)$ profile represented by the red curve in figure 14(a). The family of neutral solutions provides helpful information about expected variations in observations.

Figure 14. Idealised, non-dimensional zonal wind profiles relevant to Jupiter. Note that the domain of latitudes shown is an enlarged view near the origin of the full domain, and $L/L_d \gg 1$ is assumed such that $\kappa _0^{-2}\approx L_d^2$. (a) The solid black curve is $U(y) = \sin (2{\rm \pi} y)$. The blue curve is $U_{deep}(y)$ corresponding $M^{-1} = 1$, which amounts to a positive shift by $\beta = 0.3$. The dashed black curve is the profile $M^{-1}(y) = a+(b-a){\rm sech}^2(y)$, with $(a,b) = (0.6,1.1)$, which is a point on the neutral curve shown in figure 13; the red curve is the corresponding $U_{deep}(y)$. (b) Profiles of $U$ (solid) and $\kappa _0^{-2}Q_y$ (dashed) corresponding to the dashed $M^{-1}$ and red $U_{deep}$ profiles in (a); compare with figure 3.

Recall that for class (i), the choice of $\alpha$ is unique. Furthermore, on the neutral curve, the value of $\alpha$ must equal to the phase speed of the neutral modes (to be shown in Theorem 6.2), consistent to the observation by Read et al. (2009a) that long-wavelength Rossby waves in Saturn have the same $10^{\rm h} 34^{\rm m}$ rotation period. Consequently, for those neutral modes, the PV extrema and critical latitudes must coincide. As seen in the previous section, when $L_d$ is small, the middle critical latitude in the model is nearly sonic ($M^{-1}=1$), when the flow is nearly neutral. This situation is consistent with the observations by Read et al. (2006, 2009b) where it is revealed that Jupiter and Saturn each have at least a dozen stable or marginally stable (i.e. supersonic or sonic, $M^{-1} \leqslant 1$) PV extrema, consistent with the fact that the alternating jets on those planets persist on a decadal time scale (Porco 2003).

In addition, as first spotted by Dowling & Ingersoll (1989) and later confirmed by Marcus & Shetty (2011), Jupiter's Great Red Spot straddles a PV extremum. Numerical experiments (e.g. Dowling 1993) suggest that vortices appear at subsonic ($M<1$) latitudes, as theoretically expected for nearly neutral unstable modes. As seen in figure 1, the nature of the correlation of $U$ and $Q_y$ differs between the latitude ranges of $-30^\circ$ to $-20^\circ$ and $-50^\circ$ to $-30^\circ$, and the former range is where features such as the Great Red Spot and Oval BA are observed. As a useful exercise, one can plot $M^{-1}(y)$ profile using the data shown in figure 1 and $L_d$ estimated in figure 3, to demonstrate the presence of a hurdle that aligns with the range of latitudes $-30$ to $-20$ degrees, and that at other latitudes, $M^{-1}(y)$ is characterised by several humps that peak close to unity.

Further to this point, Dowling (2020) pointed out that Saturn has two different, persistently wavy jets sandwiched between straight jets in its northern hemisphere, the Ribbon and the Polar Hexagon, and offered the hypothesis that these are examples of ‘subsonic’ regions sandwiched between ‘sonic’ or ‘supersonic’ regions. Although the temperature fields necessary to empirically determine the stretching vorticity, and hence the PV field, at and below Saturn's cloud tops are difficult to obtain with the data in hand, the morphology of these wavy jets is well established, which suggests that an analysis couched in terms of $M^{-1}(y)$ along the lines outlined here would be instructive.

6. Mathematical proofs

Our strategy to show sufficient conditions for instability (i.e. necessary conditions for stability) is similar to the two-step procedure in Howard (1964) and Lin (2003): first prove the conditions that guarantee the existence of a neutral solution, and then establish the existence of unstable modes in its neighbourhood in parameter space. A regular mode, to be defined in the following, must exist for these processes to take place, and guaranteeing this essentially corresponds to the second class (ii) condition.

6.1. Choice of $\alpha$ for class (i) basic flows

We first note that the following theorem holds.

Theorem 6.1 Suppose the basic flow is class (i). Then there is only one $\alpha \in \mathcal {R}$ that makes $M_{\alpha }^{-1}$ continuous. Moreover, $M_{\alpha }^{-1}\neq 0$ for all $y\in \varOmega$.

This theorem can be shown as follows. The possibility of $M^{-1}_{\alpha }$ vanishing only occurs at $y_l \in Y_Q$. For class (i), $U'(y_l)\neq 0$ for all $y_l \in Y_Q$, because otherwise $M_{\alpha }^{-1}$ becomes singular, in view of the assumption that $Q''(y_l)\neq 0$ for all $y_l \in Y_Q$. L'Hôpital's rule now suggests that $M_{\alpha }^{-1}(y_l)=Q''(y_l)/U'(y_l) \neq 0$. We can also show that the $\alpha \in \mathbb {R}$ that makes $M^{-1}_{\alpha }$ continuous and one-signed is uniquely determined. Suppose there are two possible values of $\alpha$, $\alpha _1$ and $\alpha _2$, say. Since $M_{\alpha _1}^{-1}$ and $M_{\alpha _2}^{-1}$ are continuous functions that do not vanish in $\varOmega$, $(M_{\alpha _1}-M_{\alpha _2})$ is a continuous function. However, this implies that $\frac {\alpha _2-\alpha _1}{Q'}$ is a continuous function, which is not possible unless $\alpha _1=\alpha _2$.

Class (i) is an extension of class $\mathcal {H}$ in Howard (1964) and class $\mathcal {K}^+$ in Lin (2003) to the quasi-geostrophic equation. However, such strong restrictions for basic flows are not necessary to derive the hurdle theorem, as remarked in § 3.

6.2. Classification of neutral modes

Let us consider the neutral solutions, setting $c_i=0$. The key to classify neutral modes is to note that they can potentially become singular at the critical latitudes, the points at which $U(y)$ coincides with the phase speed $c_r$. Mathematically, the set of the critical latitudes, $Y_{U,c_r}=\{y \in \varOmega \,|\, U(y)=c_r\}$, are regular singular points of (2.3) when $c_i=0$. The appropriate tool for analysing the behaviour of solutions around such singularities is Frobenius’ method, from which it can be shown that $\psi$ is continuous but $\psi '$ may be discontinuous at the critical latitudes (Lin 1955). A necessary and sufficient condition for no jumps to exist at a critical latitude is that either $Q'=0$ or $\psi =0$, or both, occur there. Here, following Drazin & Reid (1981), we briefly explain the nature of the critical latitudes without going into the details of Frobenius’ method.

Because of the singularities, neutral solutions must be understood as the vanishing $c_i$ limit in unstable solutions. Let us multiply (2.3) by $\psi ^*$ and take the imaginary part, keeping finite $c_i$:

(6.1)\begin{equation} \psi^*\psi''-\psi \psi^*\,'' +\frac{2ic_i Q'}{(U-c_r)^2+c_i^2}|\psi|^2=0. \end{equation}

Integrate (6.1) across a critical latitude, $y=y_c$ (i.e. $U(y_c)=c_r$), and take the limit $c_i\rightarrow 0+$,

(6.2)\begin{equation} [\psi^*\psi'-\psi \psi^*\,']^{y_c+}_{y_c-}= \lim_{\epsilon\rightarrow 0+}\lim_{c_i\rightarrow 0+} \int^{y_c+\epsilon}_{y_c-\epsilon} \frac{-2ic_i Q'}{(U-c_r)^2+c_i^2}|\psi|^2 \,{{\rm d}y} .\end{equation}

An intuitive way to evaluate the right-hand side is as follows. If $\epsilon >0$ is sufficiently small, we may be able to assume that $Q'|\psi |^2$ is almost a constant, and that $U-c_r$ is approximately $U_c'(y-y_c)$, where $U'_c=U'(y_c)$. Then the right-hand side of (6.2) can be explicitly worked out by using the well-known formula that can be found by differentiating the arctangent function:

(6.3)\begin{equation} {-}2i(Q'|\psi|^2)|_{y=y_c}\lim_{\epsilon\rightarrow 0+}\lim_{c_i\rightarrow 0+} \frac{2}{U_c'}\arctan\left (\frac{\epsilon U_c'}{c_i}\right )={-}2i{\rm \pi} \left. \left (\frac{Q'|\psi|^2}{|U'|}\right )\right |_{y=y_c}. \end{equation}

The above argument is consistent with the viscous problem when the singularity of the inviscid solution is regularised by viscosity (Lin 1955). The regularisation by inertia is also possible (Haberman 1972; Robinson 1974), although such a situation is relevant to nonlinear equilibrium states (e.g. Deguchi & Walton 2018), here we only consider the linear problem.

The above result suggests that if multiple critical latitudes appear as $Y_{U,c_r}=\{y_1, y_2, \ldots, y_N\}$, integrating (6.1) over $\varOmega$ yields

(6.4)\begin{equation} 0=\sum_{j=1}^{N}\left. \left (\frac{Q'|\psi|^2}{|U'|}\right )\right |_{y=y_j}. \end{equation}

If $Q'=0$ or $\psi =0$ happen at all the critical latitudes, there are no jumps at all, and hence the neutral solution $\psi$ is real and $C^2_0$. Here, it is convenient to define terminology to distinguish between neutral-mode types.

1. (a) Pathological mode: empty $Y_{U,c_r}$ or $\psi$ vanishes at all critical latitudes.

2. (b) Regular mode: eigenfunction is not pathological and is $C^2_0$.

3. (c) Singular mode: eigenfunction is not $C^2_0$.

In standard textbooks such as Drazin & Reid (1981) and in previous research, there has been no distinction made between pathological modes and regular modes. We shall show in § 6.3 that when the first class (ii) condition is satisfied and the reciprocal Rossby Mach number surpasses a hurdle, neutral solutions exist with some choices of $k$, and at least one of them must be a regular mode. Moreover, if the second class (ii) condition holds, unstable modes must exist around the (non-pathological, least oscillatory) regular neutral mode (§ 6.4).

If the intersection of sets $\{U(y)\,|\,y\in \varOmega _+ \}$ and $\{U(y)\,|\,y\in \varOmega _- \}$ is non-empty, the summation in (6.4) may cancel and, hence, a singular mode may exist (thus, for class (i) there are no singular modes). For the singular modes we cannot use the standard Sturm–Liouville theory to be used in §§ 6.3 and 6.4, and for the pathological modes we have difficulty in showing neighbourhood instability.

6.3. A sufficient condition for existence of a neutral mode

Let us assume that for some $\alpha \in \mathbb {R}$ the reciprocal Rossby Mach number $M^{-1}_{\alpha }(y)$ becomes continuous. Write $\lambda =-k^2$. Then from (2.3) the neutral modes with the phase speed $c=\alpha$ must satisfy

(6.5)\begin{equation} {-}\psi'' + (L_d^{{-}2}-\kappa_{0}^2M_{\alpha}^{{-}1})\psi=\lambda \psi,\quad y \in ({-}L/2,L/2),\end{equation}

and the Dirichlet boundary conditions. This is a regular Sturm–Liouville problem with eigenvalue $\lambda$. It is well known that all eigenvalues are real, and they can be ordered as $\lambda _0<\lambda _1<\lambda _2<\dots$, where $\lambda _n \rightarrow \infty$ as $n\rightarrow \infty$. By integration by parts, it is easy to check that the associated eigenfunctions, $\psi _0, \psi _1,\psi _2,\ldots$, are real and satisfy the equation

(6.6)\begin{equation} R_{\alpha}(\psi_n)\equiv \frac{\displaystyle\int_{\varOmega} \{|\psi_n'|^2+(L_d^{{-}2}-\kappa_{0}^2M_{\alpha}^{{-}1})|\psi_n|^2\}\,{{\rm d}y}}{ \displaystyle\int_{\varOmega} |\psi_n|^2\,{\rm d} y}=\lambda_n.\end{equation}

Here the Rayleigh quotient, $R_{\alpha }$, depends on $\alpha$. It should also be noted from Sturm's oscillation theorem that the $n$th eigenfunction $\psi _n$ has $n$ zeros in the interval $(-L/2,L/2)$, a useful fact to be used below. In figure 6(a), the red curve is $\psi _0$, and the other curves may be $\psi _1, \psi _2, \psi _3$. Likewise the modes shown in figure 16(b) are the zeroth and first modes.

The minimum eigenvalue $\lambda _0$ can be found by the optimisation problem

(6.7)\begin{equation} \lambda_0=\min_{\phi\in C^2_0}R_{\alpha}(\phi),\end{equation}

where the unique minimiser is the zeroth eigenfunction $\phi =\psi _0$. Note that there is no problem to extend the search space to

(6.8)\begin{equation} H_0^1=\left\{\phi \left | \int_{\varOmega}|\phi|^2\,{\rm d} y<\infty, \int_{\varOmega}|\phi'|^2\,{\rm d} y<\infty, \phi({-}L/2)=\phi(L/2)=0\right. \right \} \end{equation}

by a density argument. Here, following the usual notation in mathematics, $H$ implies that the space is a Hilbert space, the superscript $1$ implies that the square integrability of the first derivative and the subscript $0$ implies that the Dirichlet boundary conditions are satisfied.

A sufficient condition for existence of a neutral mode is then summarised as follows.

Theorem 6.2 Suppose there is $\alpha \in \mathbb {R}$ that makes $M^{-1}_{\alpha }$ continuous. If there exists $g(y)\in H^1_0$ such that $R_{\alpha }(g)\leq 0$, there is a neutral mode with the phase speed $c_r=\alpha$.

If the assumption of the theorem is met, a neutral mode can be found because $\lambda _0=\min _{\phi \in H_0^1}R(\phi ) \leq R(g) \leq 0$ implies that the minimiser of (6.7) is the neutral eigenfunction having the wavenumber $k=k_0\equiv \sqrt {-\lambda _0}\geq 0$. This neutral mode is $\psi _0$ introduced earlier. The neutral curves shown in figures 7(a), 10 and 13 are indeed determined by $\psi _0$ and, crucially, this mode must be a regular mode.

Since $c_r=\alpha$ for $\psi _0$, the critical latitudes (the points $y_j$ at which $U(y_j)=c_r$ is satisfied) must appear on the PV extrema (i.e. $Y_{U,c_r} \subset Y_Q$). However, for class (i) basic flows, the stronger result $Y_Q=Y_{U,c_r}$ can be shown for $\psi _0$, because $Y_Q=Y_{U,\alpha }$ as remarked just below (3.5), and there is only one choice of $\alpha$ from Theorem 6.1.

In the next section, we shall show the following theorem.

Theorem 6.3 Suppose the basic flow is class (ii), and fix $\alpha$ so that the class (ii) conditions are satisfied. If there exists $g(y)\in H^1_0$ such that $R_{\alpha }(g)< 0$, the flow is unstable. Furthermore, if the basic flow is class (i), a necessary and sufficient condition for stability is that no such $g(y)\in H^1_0$ exists.

The second half of the theorem is somewhat similar to that derived by Howard (1964) and Lin (2003) for shear flows, but the first half is entirely new. We also remark that in the case of the Rayleigh equation, it can be proved that the phase speed of the neutral solution is in the range of $U(y)$, denoted $\mathcal {R}$, by Howard's semicircle theorem, and that there are no pathological modes when $k>0$. Moreover, when $k=0$, a regular mode can be found analytically. Those facts are used, for example, in proofs by Howard (1964), Balmforth & Morrison (1999) and Hirota et al. (2014), but they do not in general apply to stability problems more complex than the Rayleigh equation. This is essentially the reason why Tung (1981) struggled to incorporate the effect of non-zero $\beta$ into the theory, but as we will see in the next section, the solution is, in fact, simple.

6.4. Existence of an unstable mode

Here we show Theorem 6.3. The proof is rather technical and readers who are not interested in the mathematical details may skip this subsection without losing the thread of the discussion. We first note that upon writing $\gamma \equiv -L_d^{-2}-k^2=\lambda -L_d^{-2}$, (2.3) becomes

(6.9)\begin{equation} \psi''+\left(\frac{Q'}{U-c}+\gamma \right ) \psi =0.\end{equation}

The corresponding dispersion relation, $F(c,\gamma )=0$, can be formulated by two linearly independent solutions of (6.9). They are analytic functions in the upper or lower half complex $c$-plane, and behave regularly with respect to $\gamma$, and so does $F(c,\gamma )$. A neutral solution is obtained when $c_i$ is brought close to zero in the dispersion relation. From the symmetry of the inviscid problem, neutral solutions always exist as pairs, i.e. the $c_i=0+$ mode and the $c_i=0-$ mode.

Let us suppose that the assumptions of Theorem 6.3 are satisfied. From Theorem 6.2 we know that there is a regular neutral mode $\psi _0$ with wavenumber $k_0=\sqrt {-\lambda _0}>0$ and phase speed $\alpha$. This mode satisfies

(6.10)\begin{equation} \mathcal{L}\psi_0\equiv \psi_0''+\left(\frac{Q'}{U-\alpha}+\gamma_0\right) \psi_0=0. \end{equation}

Here, we define the linear operator $\mathcal {L}$ for later use.

In view of the argument just below (6.9), we can compute the coefficients of the Taylor expansion of $c(\gamma )$ around the neutral mode in the upper half complex plane. Writing $\gamma _0\equiv \lambda _0-L_d^{-2}$, around the neutral mode, the following expansion holds,

(6.11)\begin{equation} c_i(\gamma)=\left.\frac{{\rm d} c_i}{{\rm d}\lambda}\right|_0(\gamma-\gamma_0)+O((\gamma-\gamma_0)^2). \end{equation}

Here and hereafter, the symbol ‘$|_0$’ attached to a quantity means that it is evaluated at the $c_i=0+$ neutral mode. The expression (6.11) is valid as long as $c_i$ does not become negative. Our goal below is to show ${{\rm d} c_i}/{{\rm d}\gamma }|_0\neq 0$; then (6.11) implies that an unstable mode can be found by slightly varying $\gamma$ from the neutral value.

Differentiate (6.9) with respect to $\gamma$,

(6.12)\begin{equation} \psi_{\gamma}''+\left(\frac{Q'}{U-c}+\gamma \right ) \psi_{\gamma} +\left (1+\left (\frac{Q'}{U-c}\right )_{\gamma} \right ) \psi=0, \end{equation}

where the subscript $\gamma$ denotes the differentiation. Evaluate this equation at the neutral parameter,

(6.13)\begin{equation} \mathcal{L} \psi_{\gamma}|_0+\left (1+\left.\left(\frac{Q'}{U-c}\right)_{\gamma}\right|_0\right) \psi_0=0. \end{equation}

Now, combine (6.10) and (6.13) to obtain

(6.14)\begin{equation} 0=\int_{\varOmega}(\psi_0\mathcal{L} \psi_{\gamma}|_0-\psi_{\gamma}|_0\mathcal{L} \psi_0)\,{{\rm d}y}+\int_{\varOmega}\left(1+\left.\left(\frac{Q'}{U-c}\right)_{\gamma}\right |_0 \right ) \psi_0^2\,{\rm d} y. \end{equation}

The first integral vanishes after integration by parts. Therefore

(6.15)\begin{equation} \left. \frac{{\rm d} c}{{\rm d}\gamma}\right |_0={-}\frac{\displaystyle\int_{\varOmega} \psi_0^2\,{\rm d} y}{K}={-}\frac{(K_r-iK_i)\displaystyle\int_{\varOmega} \psi_0^2\,{\rm d} y}{K_r^2+K_i^2}, \end{equation}

where

(6.16)\begin{equation} K=K_r+iK_i=\lim_{c_i \rightarrow 0+}\int_{\varOmega} \frac{Q'\psi_0^2}{(U-\alpha-ic_i)^2}{{\rm d}y}. \end{equation}

The limit can be worked out as

(6.17)\begin{equation} K=\int\kern-10pt_{\varOmega} \frac{Q'\psi_0^2}{(U-\alpha)^2}{{\rm d}y}+{\rm i}{\rm \pi}\kappa_0^2 \sum_{j=1}^N \left.\left(\frac{M_{\alpha}^{{-}1}\psi_0^2}{|U'|} \right) \right|_{y=y_j}, \end{equation}

noting that the integrand of (6.16) becomes singular at $Y_{U,\alpha }=\{y_1, y_2, \ldots, y_N\}$. The dashed integral represents the Cauchy principle value integral.

Here, since $\alpha \in \mathcal {R}$, the set $Y_{U,\alpha }$ is non-empty. The terms in the summation in (6.17) cannot cancel out because all the $M_{\alpha }^{-1}(y_j)$ are one-signed from the second Class (ii) condition. Also $U'(y_j)\neq 0$ for all $j$, as otherwise $M_{\alpha }^{-1}(y)$ cannot be continuous. Moreover, $\psi _0$ does not vanish in $\varOmega$ because it is the least oscillatory eigenmode, as remarked just below (6.6). From (6.17) $K_i$ is non-zero, and therefore we can conclude $({{\rm d} c_i}/{{\rm d}\gamma }) |_0\neq 0$ using (6.15).

We still need to show the latter half of Theorem 6.3 for class (i) basic flows. It is sufficient to consider the case $M_{\alpha }^{-1}(y)>0$ for all $y$, as otherwise the flow must be stable from Theorem 3.1. For any unstable mode $\psi$ we can deduce (A5), which yields $R_{\alpha }(\psi )<0$. Thus, if we suppose there is no $g\in H^1_0$ such that $R_{\alpha }(g)<0$, then there should be no unstable modes.

6.5. Simple integral-based stability conditions

If a function $g\in H^1_0$ that makes the Rayleigh quotient negative is found, the existence of instability is guaranteed by Theorem 6.3. The remaining question for deriving simple necessary conditions for stability is how to choose $g$. A good test function is $g=\varphi _0$, the function introduced in (3.4b). Using (6.6), the condition $R_{\alpha }(\varphi _0)<0$ becomes

(6.18)\begin{equation} 1<\frac{2}{L}\int^{{L}/{2}}_{-{L}/{2}}M_{\alpha}^{{-}1} \cos^2 \left(\frac{{\rm \pi} y}{L} \right) {{\rm d}y}. \end{equation}

Note that this result depends on the boundary conditions. For example, if periodic boundary conditions are used, the eigenvalue problem (3.2) produces (4.1b), and thus the condition (6.18) must be replaced by

(6.19)\begin{equation} 1<\frac{1}{L}\int^{{L}/{2}}_{-{L}/{2}}M_{\alpha}^{{-}1}\,{{\rm d}y}.\end{equation}

6.6. A simple stability condition using a hurdle

The above conditions require us to examine the profile of $M^{-1}_{\alpha }(y)$ across the entire domain, $\varOmega$. However, it turns out we can show the existence of instability by just checking a local part of the basic flow. Consider a comparison problem in a subinterval $[y_1,y_2] \subset \varOmega$ with some $\tilde {M}^{-1}(y)$ profile.

(6.20)\begin{equation} \tilde{\psi}''-(\tilde{k}^2+L_d^{{-}2})\tilde{\psi}+\kappa_{0}^2\tilde{M}^{{-}1} \tilde{\psi}=0,\quad y \in (y_1,y_2). \end{equation}

We impose the Dirichlet condition, $\tilde {\psi }(y_1)=\tilde {\psi }(y_2)=0$. If this problem has a neutral mode, $\tilde {\psi }_0(y)$ say, and $\tilde {M}^{-1}< M_{\alpha }^{-1}$ on $[y_1,y_2]$, the original problem (6.5) must be unstable. The reason is as follows. We set the test function as

(6.21)\begin{equation} g(y) =\begin{cases} \tilde{\psi}_0(y), & \text{if}\ y\in [y_1,y_2], \\ 0, & \text{otherwise}, \end{cases}\end{equation}

using the neutral mode. Then the Rayleigh quotient can be computed as

(6.22)\begin{align} R(g) &= \frac{\displaystyle\int_{y_1}^{y_2} |\tilde{\psi}_0'|^2+ (L_d^{{-}2}-\kappa_{0}^2M_{\alpha}^{{-}1})|\tilde{\psi}_0|^2\,{\rm d} y}{\displaystyle\int_{y_1}^{y_2} |\tilde{\psi}_0|^2\,{\rm d} y} \nonumber\\ &< \frac{\displaystyle\int_{y_1}^{y_2} |\tilde{\psi}_0'|^2+ (L_d^{{-}2}-\kappa_{0}^2\tilde{M}^{{-}1})|\tilde{\psi}_0|^2\,{\rm d} y}{\displaystyle\int_{y_1}^{y_2} |\tilde{\psi}_0|^2\,{\rm d} y}={-}\tilde{k}^2\leq 0, \end{align}

where $\tilde {k}$ is the wavenumber of the neutral mode $\tilde {\psi }_0(y)$. A particularly convenient comparison problem is when $\tilde {M}^{-1}$ is a constant, since the problem can be solved analytically using trigonometric functions. It is easy to check that a neutral mode with a wavenumber $\tilde {k}$ can be found when

(6.23)\begin{equation} \tilde{M}^{{-}1}=\frac{\dfrac{{\rm \pi}^2}{(y_2-y_1)^2}+(\tilde{k}^2+L_d^{{-}2})}{\kappa_{0}^2}. \end{equation}

We want to make $\tilde {M}^{-1}$ as small as possible to produce instability criteria as sharp as possible, so we let $\tilde {k}\rightarrow 0$. Then the condition $\tilde {M}^{-1}< M_{\alpha }^{-1}$ on $[y_1,y_2]$ yields Theorem 3.2. In brief, the hurdle theorem corresponds to selecting the test function $g(y)$ in (6.21) with $\tilde {\psi }_0(y)=\sin ({\rm \pi} \frac {y-y_1}{y_2-y_1})$. It is possible to design test functions that provide greater sensitivity to the conditions, but in this article we have chosen to prioritise simplicity of the criteria.

7. Conclusion and discussion

Rayleigh's inflection-point theorem of 1880 ushered in a century-long series of sufficient conditions for the stability of inviscid, parallel shear flows. Yet, well into the 21st century, sufficient conditions for instability (or, equivalently, necessary conditions for stability) have remained elusive for applied researchers. We have developed a simple methodology that allows for the detection of inviscid instability, inspired by empirical observations of the stability of Jupiter's alternating jets. All that is needed to operate our method is a comparison of $h$ defined in (3.9) with the profile of the reciprocal Rossby Mach number, $M_{\alpha }^{-1}(y)$, which can be easily calculated from the basic flow (see § 3.3). The true stability boundary in the parameter space should be found between the unstable regions detected by Theorem 3.2 and the stable regions deduced by the KA criteria.

Our method relies solely on the fundamental tools of Sturm–Liouville theory, but exhibits a high level of applicability beyond Rayleigh's equation framework. Similar hurdle stability conditions guaranteeing instability are likely to be found for various inviscid shear-flow stability problems in science and engineering, including hypersonic, stratified, magneto-hydrodynamic and/or non-Newtonian flows. To the best of the authors’ knowledge, the possibility of such a general and practical instability theory has not been pointed out previously. The hurdle theory could even yield new insights into classical Rayleigh equation problems as illustrated in Appendix C, where an internally heated vertical channel is analysed.

As noted in § 1, in order for the sufficient condition of instability (Theorem 3.2) to apply, the basic flow must belong to a class with certain favourable properties. We assume there is at least one PV extremum in the domain and identify two basic-flow classes of interest, for which $M^{-1}_{\alpha }(y)$ is continuous in some reference zonal wind shift, $\alpha$. From the perspective of the general inviscid stability problem, one of the novelties of our work is the extension of the applicability of stability conditions through the introduction of new classifications. As illustrated by Theorem 6.3, our sufficient condition of instability apply to class (ii), which covers almost all monotonic flows and a wide range of non-monotonic flows. If the basic flow belongs to this class, it can be asserted that there exists the least oscillatory regular neutral mode (§ 6.3) and that instability occurs when the parameters are slightly altered (§ 6.4). Theorem 6.3 requires a test function, but the stability condition can be simplified to the hurdle form using a kind of comparison principle (§ 6.6). A convenient method to determine whether the basic flow belongs to class (ii) is as follows. (1) Plot $U(y)$ and mark all the points at the latitudes where $Q'=0$, $y=y_j$ say. Then for each point, draw a line at constant $U=U(y_j)$. If all of these lines intersect the $U(y)$ curve at points which are not marked, then the profile is not in class (ii). (2) If the profile passes the first test, then check whether there is an odd number of zeros of $Q'$ in between any successive pair of zeros of $U-\alpha$. If so, this is also not in class (ii). If the profile passes both tests, the profile is in class (ii).

In general, class (ii) admits more than one reference-frame shift, $\alpha$ (see § 4.2). However, for a subclass of it, referred to as class (i), $\alpha$ is uniquely determined (Theorem 6.1). For the latter class, $M^{-1}(y)$, which can now be safely unsubscripted, does not change sign. Class (i) is approximately identified with Jupiter and Saturn. An important feature of class (i) is that on the neutral curve the PV extrema must coincide with the critical latitudes, where the zonal wind speed matches with the phase speed of the neutral mode $\alpha$. Around the neutral parameters, vortices are produced near the critical latitudes, although mathematically we can show that all neutral modes have no singularity there. Furthermore, for class (i), Theorem 3.2 gives the necessary and sufficient condition for the stability when the limit $L_d\rightarrow 0$ is taken.

The stability results for the linearised $1\frac {3}{4}$ layer model are summarised as follows.

1. (a) If $Q'$ has no zeros, the flow is stable (Charney–Stern).

2. (b) Even if $Q'$ has zeros, if the flow is class (i), it is stable if either $M^{-1}\leq 0$ everywhere (KA-I) or $ 0 \leq M^{-1}\leq 1$ everywhere (KA-II).

3. (c) If the flow is class (ii), it is unstable if there are $\alpha, h$ such that $M_{\alpha }^{-1}> h$ in the interval where the hurdle $h$ is defined (Hurdle theorem).

If the flow is not class (ii), then numerical methods are currently the primary recourse to determining stability.

Our results also offer new insights into the theoretical understanding of pattern formation in planetary atmospheres. The KA stability theorems were re-expressed by Dowling (2014) in terms the Rossby Mach number, revealing a key attribute that supersonic PV extrema are stable. However, this left the natural question of how subsonic a PV extremum must be before becoming unstable, which is answered precisely by Theorem 3.2, and illustrated with detailed model analyses in §§ 4 and 5. In the case where $M^{-1}$ is a constant, the sonic condition, $M^{-1}=1$, must give the stability boundary, according to Theorems 3.1 and 3.2. Therefore, the conjecture of Stamp & Dowling (1993) is mathematically confirmed to be correct. In § 4 we extended the Stamp & Dowling (1993) model to study the case where $M^{-1}$ is not constant. The ${\rm sech}^2$ bump profile considered for $M^{-1}(y)$ (see (4.3)) is interesting from the perspective of planetary physics because the behaviour of neutral curves can be analysed analytically to some extent. There exist qualitatively distinct neutral states, oscillatory and localised modes. In the context of Jupiter's atmosphere, only the latter appears, and the KA-II condition becomes almost sharp (figure 13). This implies that in atmospheres sustained over long periods, the $M^{-1}$ profile cannot become too subsonic, aligning with observational facts. However, a closer look in figure 1 reveals that the $Q_y$ curve sits above $U$ in the latitude range of $-$30 to $-$20 degrees. Unsteady dynamics may occur at neutral or subsonic PV extrema there (see Read et al. 2006, 2009b).

From the $M^{-1}$ profile under the assumption of neutrality, in the context of $1$–${3}/{4}$ layer dynamics, we can calculate $U_{deep}(y)$ from (2.4). It is noteworthy that the estimate of $U_{deep}(y)$ has been historically important for probing the deep jets on Jupiter (Dowling 1995b), though additional considerations for Jupiter in combination with Juno gravity inversions of interior circulations, and for Saturn in combination with ring-wave seismology (i.e. kronoseismology), are warranted in the future.

The quasi-correlations of the zonal flow $U$ and the PV gradient $Q'$ seen in figure 1 suggests that the Jupiter's alternate jets are formed to be nearly linearly neutral with respect to inviscid Rossby wave instability. To understand why such a mean flow is achieved, it is necessary to consider the nonlinear evolution of the disturbances and their mutual interaction with the mean flow. It is numerically shown in Dowling (1993) that when Jupiter's observed basic flow is made the initial condition for a local unforced shallow-water model, it will rapidly evolve to become marginally stable. In addition, slowly moving, planetary-scale thermal features have been regularly observed on Jupiter (see Fisher et al. 2016 and references therein). These lines of evidence suggest that a gas giant's alternating jets tend to evolve until the longest Rossby waves are coherent across them, becoming stationary with respect to deep-seated, large-scale pressure anomalies; this has been called a ‘princess and the pea’ phenomenon (Stamp & Dowling 1993).

The spontaneous generation mechanism of zonal jets has long remained a subject of debate among experts, as documented in Read (2024). The central inquiry lies in understanding how the energy provided by the heat from the Sun or the interior of gas giants is converted into the kinetic energy of zonal and equatorial jets. While some numerical successes have emerged over the decade (see Schneider & Liu 2009; Lian & Showman 2010), a comprehensive and rational explanation is still lacking. Decomposing the flow into non-axisymmetric and axisymmetric components marks the first step in observing how Rossby waves transfer angular momentum within zonal mean flows, facilitated by Reynolds stress (Andrews & McIntyre 1976, 1978). The observed strong correlation between zonal flow and Reynolds stress suggests that this is indeed an indispensable mechanism (e.g. Ingersoll et al. 1981). Coupling the mean flow with linear perturbation equations through Reynolds flux terms is often called the mean-field approximation or quasi-linear theory in the modelling community (e.g. O'Gorman & Schneider 2007). It is noteworthy that studies on near-wall turbulence have demonstrated how the interaction between mean vortex fields and inviscid neutral waves rationally elucidates the generation of streaks (Wang, Gibson & Waleffe 2007; Hall & Sherwin 2010). Recent investigations have further pinpointed the presence of neutral waves within large-scale laminar–turbulent patterns in rotating flows (Wang et al. 2022). Thus, we expect that the formation and evolution of Jupiter's alternating jets might be similarly explored by coupling the mean flow equations for the planet's atmosphere and interior with the inviscid stability problem of the mean flow, along the lines of similar theories developed in near-wall turbulence.