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 Reference Vallis2017). 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 (Reference Ingersoll and Cuong1981). 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 Reference Majda and Wang2005). 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 Reference Flierl, Morrison and VilasurSwaminathan2019).

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 Reference Dowling1995b; Read et al. Reference Read, Gierasch, Conrath, Simon-Miller, Fouchet and Yamazaki2006, Reference Read, Conrath, Fletcher, Gierasch, Simon-Miller and Zuchowski2009b,Reference Read, Dowling and Schuberta). 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 Reference Numata, Ball and Dewar2007; Zhu, Zhou & Dodin Reference Zhu, Zhou and Dodin2018). 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.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 (Reference Dowling2014, Reference Dowling2020). 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.

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 Reference Dowling2014). 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 Reference Lindzen and Tung1978). 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 Reference Dowling2020; Afanasyev & Dowling Reference Afanasyev and Dowling2022). As mentioned in § 1, observations support that $M^{-1}$ is around unity in Jupiter's and Saturn's atmosphere (Dowling Reference Dowling1993; Read et al. Reference Read, Gierasch, Conrath, Simon-Miller, Fouchet and Yamazaki2006, Reference Read, Conrath, Fletcher, Gierasch, Simon-Miller and Zuchowski2009b,Reference Read, Dowling and Schuberta). 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. Reference Hammel1995) 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 Reference Du, Dowling and Bradley2015). 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.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 (Reference Stamp and Dowling1993). 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 (Reference Deguchi2021)). 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).

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

We first note that the following theorem holds.

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 (Reference Howard1964) and class $\mathcal {K}^+$ in Lin (Reference Lin2003) 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 Reference Lin1955). 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 (Reference Drazin and Reid1981), 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 Reference Lin1955). The regularisation by inertia is also possible (Haberman Reference Haberman1972; Robinson Reference Robinson1974), although such a situation is relevant to nonlinear equilibrium states (e.g. Deguchi & Walton Reference Deguchi and Walton2018), 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 (Reference Drazin and Reid1981) 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.

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.

The second half of the theorem is somewhat similar to that derived by Howard (Reference Howard1964) and Lin (Reference Lin2003) 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 (Reference Howard1964), Balmforth & Morrison (Reference Balmforth and Morrison1999) and Hirota et al. (Reference Hirota, Morrison and Hattori2014), but they do not in general apply to stability problems more complex than the Rayleigh equation. This is essentially the reason why Tung (Reference Tung1981) 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.