Analysis of inviscid shear instability of axisymmetric flows

Abstract

Simple analytical criteria are derived to determine whether axisymmetric base flows in annuli and pipes are stable or unstable. Both axisymmetric and non-axisymmetric inviscid disturbances are considered. Our sufficient condition for stability improves upon the classical result of Batchelor & Gill (1962) J. Fluid Mech. 14(4), 529–551 following the idea of the second Kelvin–Arnol’d stability theorem. A novel sufficient condition for instability is also derived by extending the recently proposed hurdle theorem for parallel flows (Deguchi et al. 2024 J. Fluid Mech. 997, A25). These analytical criteria are applied to annular and pipe model flows and are shown to effectively predict the neutral parameters obtained from eigenvalue computations of the stability problem.

1. Introduction

Axisymmetric flows arise in a variety of applications, including flows through annuli and pipes as well as round jets, and their stability has long been of considerable interest. In most practical situations, the Reynolds number is high, so viscous effects may be neglected. The inviscid approximation greatly simplifies the stability equations and often leads to fruitful mathematical results. An important contribution to the inviscid stability of axisymmetric base flows was made by Batchelor & Gill (1962), who derived a convenient sufficient condition for stability. This result can be viewed as an extension of the well-known Rayleigh–Fjortoft stability condition for parallel flows (Rayleigh 1880; Fjørtoft 1950). Although significant progress has since been made in the inviscid stability analysis of parallel base flows, somewhat surprisingly, applications to axisymmetric base flows remain scarce.

Arnol’d (1966) showed that two distinct sufficient conditions for stability exist through a pseudoenergy-based analysis, which is now known as Arnol’d’s method (see Shepherd 1993 also). Earlier, Lord Kelvin had also suggested the presence of two such conditions (Thomson 1880). Owing to this historical development, those stability conditions are called the Kelvin–Arnol’d first and second shear-stability theorems (KA-I and KA-II). The Rayleigh–Fjørtoft condition corresponds to KA-I, as does the result of Batchelor & Gill (1962). For parallel flows, the KA-II conditions have also been well studied, particularly in the context of instabilities in planetary atmospheres (see Stamp & Dowling 1993; Dowling 2020; Read & Dowling 2026, for example). Therefore, one might expect that the KA-II condition could also be applied to axisymmetric base flows to yield useful stability criteria; however, no such results have been explicitly derived in the literature.

Except for a few special cases, the KA stability conditions generally fail to sharply detect points of neutral stability, as first demonstrated in the pioneering work of Tollmien (1935). This limitation has motivated efforts to identify classes of base flows that admit sharp stability criteria and to simplify those criteria as much as possible. For example, Rosenbluth & Simon (1964) and Balmforth & Morrison (1999) employed the Nyquist method, whereas Barston (1991) and Hirota, Morrison & Hattori (2014) drew inspiration from the properties of certain quadratic forms. Despite being more numerically tractable than the full eigenvalue problem, these approaches still rely on the solution to a Fredholm integral equation or similarly complex criteria.

Recently, Deguchi, Hirota & Dowling (2024) derived a simple sufficient condition for instability. The condition essentially depends on whether the reciprocal Rossby–Mach number, derived from the base flow, surpasses a threshold called the hurdle. Since this can be conveniently assessed using a simple graphical procedure, the result is referred to as the hurdle theorem. Although this condition does not sharply determine the neutral point, plotting the regions in parameter space where the hurdle theorem and the KA criteria hold allows one to roughly estimate the location of neutral points. In practice, this helps restrict the parameter range that needs to be explored in the eigenvalue problem.

The proof by Deguchi et al. (2024) follows the two-step approach that became popular after the seminal work of Howard (1964): (i) first, demonstrate the existence of a neutral solution, and (ii) then show that an unstable mode emerges when this neutral state is perturbed. The essence of this approach can already be found in Tollmien (1935). A mathematically rigorous proof of step (ii) was later provided by Lin (2003), and more recently simplified by Kumar & Ożański (2025). In step (i), the neutral solution is constructed through minimisation of the Rayleigh quotient. The central idea of Deguchi et al. (2024) is to analytically derive a simple condition ensuring the existence of a neutral solution by estimating the minimum value with a suitable trial function.

In this paper, we extend the sufficient stability condition of Batchelor & Gill (1962) and derive a simple sufficient condition for instability for axisymmetric base flows. The former is obtained following the aforementioned KA-II framework, while the latter is derived by generalising the hurdle theorem. For annular flows, the analysis proceeds in much the same way as for parallel flows, although the azimuthal wavenumber of the perturbations must be treated carefully. Pipe flows, however, require additional consideration, and the conditions must be modified from those for the annular case to obtain satisfactory results.

The proposed stability conditions are applied to representative model flows to assess their practical usefulness. More specifically, we consider cases where shear flow and thermal convection coexist. Parallel shear flows, whether driven by a pressure gradient or by moving boundaries, tend to be inviscidly stable, as exemplified by plane Couette flow or plane Poiseuille flow. The first model flow we consider is sliding Couette flow (an annular version of Couette flow, see Gittler 1993; Deguchi & Nagata 2011), coupled with thermal convection due to the temperature difference between the walls. The flow is designed so that, in the absence of the wall sliding effect, the system reduces to natural convection in a vertical annulus, a configuration frequently studied in the natural convection literature (Choi & Korpela 1980; Yao & Rogers 1989; Kang, Yang & Mutabazi 2015; Wang & Chen 2022). The second model is a vertically oriented Hagen–Poiseuille flow with homogeneous internal heating. This model can be regarded as the pipe analogue of the channel model flow used to test the hurdle theorem in Deguchi et al. (2024), and it was also investigated by Senoo, Deguchi & Nagata (2012) and Marensi, He & Willis (2021). The performance of the stability conditions will be assessed by comparing their predictions with the numerically obtained neutral points. Those inviscid stability results are further compared with linearised Navier–Stokes computations.

The structure of the paper is as follows. In the next section, we formulate the stability problem for axisymmetric base flows. Numerical methods for both inviscid and viscous stability problems are also introduced. Section 3 summarises the sufficient conditions for stability and instability derived in this paper. These conditions are applied to model flows in annuli and pipes in § 4, with their derivation presented in § 5. Finally, § 6 presents the conclusions. In the same section, we also discuss the implications of the present results for round jets, which originally motivated the work of Batchelor & Gill (1962).

2. Formulation of the problem

2.1. Inviscid stability problem for axisymmetric base flows

Let $\boldsymbol{u}=u_r\boldsymbol{e}_r+u_{\varphi }\boldsymbol{e}_{\varphi }+u_z\boldsymbol{e}_z$ and $p$ denote the velocity and pressure of an incompressible fluid in cylindrical coordinates $(r,\varphi ,z)$ , with $\boldsymbol{e}_r$ , $\boldsymbol{e}_{\varphi }$ and $\boldsymbol{e}_z$ the associated unit vectors. Periodicity of $2\pi /k$ is imposed in the $z$ direction, with $k$ denoting the axial wavenumber. We assume a steady axisymmetric base flow $\boldsymbol{u}=U(r)\boldsymbol{e}_z$ directed along the axial direction, and assess its stability by examining the evolution of infinitesimal perturbations $\tilde {\boldsymbol{u}}=\tilde {u}_r\boldsymbol{e}_r+\tilde {u}_{\varphi }\boldsymbol{e}_{\varphi }+\tilde {u}_z\boldsymbol{e}_z$ superposed on it. The azimuthal wavenumber of the perturbation is denoted by $n$ . The physical requirement of $2\pi$ periodicity in the $\varphi$ direction implies that $n$ must be an integer.

The stability of the flow can be analysed by solving the linearised Navier–Stokes equations; if viscosity is neglected, the analysis can begin with the linearised Euler equations. Batchelor & Gill (1962) showed that the latter inviscid stability problem boils down to solving the differential equation

(2.1) \begin{align} \left\{\frac {r}{N^2+r^2}(rG)^{\prime} \right \}^{\prime}-k^2G+\frac {Q^{\prime}}{U-c}rG=0,\qquad Q=\frac {-rU^{\prime}}{N^2+ r^2}, \end{align}

with suitable boundary conditions. Here, a prime denotes ordinary differentiation with respect to radius, $r$ . For a fixed wavenumber ratio $N = n/k$ , (2.1) together with suitable boundary conditions constitutes an eigenvalue problem, with the complex phase speed $c=c_r+ic_i$ as the eigenvalue. By symmetry, it suffices to consider $k \gt 0$ and $n \geqslant 0$ . From the eigenfunction $G(r)$ , the perturbation field can be reconstructed as

(2.2) \begin{align} \,[\tilde {u}_r,\tilde {u}_{\varphi },\tilde {u}_{z},\tilde {p}]=[iG(r),H(r),F(r),P(r)]e^{in\varphi +ik(z-ct)}+\text{c.c.}, \end{align}

where

(2.3) \begin{align} \left [ \begin{array}{c} F\\ H \end{array} \right ] &= \frac {1}{k(N^2+r^2)} \left [ \begin{array}{cc} -N & -r\\ r & -N \end{array} \right ] \left [ \begin{array}{c} NU^{\prime}G/(U-c)\\ rG^{\prime}+G \end{array} \right ], \end{align}

(2.4) \begin{align} P &=-(U-c)F-k^{-1}U^{\prime}G, \end{align}

and c.c. denotes the complex conjugate. The existence of an eigenvalue with $c_i \gt 0$ indicates the presence of disturbances that grow exponentially with time $t$ , and hence the flow is unstable. The flow is said to be stable when no such unstable mode exists (the mode with $c_i=0$ is neutrally stable).

The general properties of the inviscid stability problem (2.1) are the main focus of this paper. For concreteness, however, we examine the stability conditions for specific base flows. The model base flows to be introduced in § 4 can be described using the Navier–Stokes equations together with the Boussinesq approximation. We specifically consider the case of an infinitesimally small Prandtl number, for which the stability problem reduces to the linearised Navier–Stokes equations. Note that, since the Prandtl number is the ratio of viscous to thermal diffusion, a small value implies that the diffusion term dominates advection in the temperature equation. Consequently, temperature perturbations decay and the buoyancy term is removed from the momentum equations for perturbations (see Lignières 1999 for example).

Section 4.1 examines a flow between coaxial cylinders, while § 4.2 studies a pipe flow. For the annular geometry, we set the inner and outer cylinder walls at $r=r_i$ and $r=r_o$ . The radius ratio $\eta =r_i/r_o \leqslant 1$ is the only parameter needed to specify the geometry. Under the non-dimensionalisation in which the gap is set to 2

(2.5) \begin{align} r_i=2\eta /(1-\eta ), \qquad r_o=2/(1-\eta ). \end{align}

The no-penetration boundary condition then requires that $G$ vanishes at the walls.

For the pipe flow problem, physical quantities are required to be regular at the centreline. If we assume that the perturbation velocity field is Taylor expandable about the pipe axis in Cartesian coordinates, then $G$ has the following expansion (see Eisen, Heinrichs & Witsch 1991):

(2.6a) \begin{align} G=\mathcal{G}_1r+\mathcal{G}_3r^3+\ldots , \,\,\,\,&\text{if}\,\,\,\, n= 0, \end{align}

(2.6b) \begin{align} G=\mathcal{G}_{n-1}r^{n-1}+\mathcal{G}_{n+1}r^{n+1}+\ldots \,\,\,\,&\text{if}\,\,\,\, n\neq 0. \end{align}

Therefore, the eigenfunction for the pipe problem is expected to satisfy

(2.7a) \begin{align} G^{\prime}=0, \,\,\,\,\text{if}\,\,\,\, n= 1, \end{align}

(2.7b) \begin{align} G=0 \,\,\,\,\text{if}\,\,\,\, n\neq 1, \end{align}

at $r=0$ . Note that (2.7) is not a boundary condition of (2.1). This conclusion follows from the standard singularity theory of ordinary differential equations. As commented in Lessen & Singh (1973), the method of Frobenius expansion yields that one of the linearly independent solutions is singular, while the other, non-singular solution satisfies (2.7).

2.2. Numerical method

The eigenvalue problem (2.1) can be solved using the Chebyshev-collocation method, in which $G(r)$ is expanded by the basis functions $\varPsi _l(r)$

(2.8) \begin{align} G(r)=\sum _{l=0}^L \hat {G}_l \varPsi _l(r). \end{align}

The annular case is straightforward. We use the basis $\varPsi _l(r)=(1-y^2)T_l(y)$ , where $T_l(y)$ are the Chebyshev polynomials in the shifted variable $y=r-r_m$ . To map $y$ into the interval $(-1,1)$ , we set $r_m=(r_o+r_i)/2$ , the mid-gap. Collocation points are chosen as $r_j=r_m+\cos (({j+1})/({L+2})\pi )$ , $j=0,1,\ldots ,L$ .

For the pipe case, we checked the consistency of the numerical codes using the following two methods:

1. (i) Expand $G$ in modified Chebyshev polynomials (see Deguchi & Walton 2013)

   (2.9) \begin{align} \varPsi _l(r)= \left \{ \begin{array}{c} (1-r^2)T_{2l+1}(r),\qquad \text{if $n$ is even},\qquad \\ (1-r^2)T_{2l}(r),\qquad \,\,\, \text{if $n$ is odd}. \end{array} \right . \end{align}
   The governing equation is then evaluated at the collocation points $r_j=\cos ((j+1)\pi /(2L+4)), j=0,1,\ldots ,L$ .

2. (ii) Treat the domain as an annulus with a virtual inner cylinder of small radius $\epsilon \gt 0$ , and impose condition (2.7) there. The code developed in Song & Dong (2023) is used, which employs the fourth-order compact finite-difference scheme (see Malik 1990). We choose $\epsilon =0.001$ .

The results of the inviscid analysis are compared with the calculations of the following linearised Navier–Stokes equations

(2.10a) \begin{align} k(U-c)G &=P^{\prime}-\frac {i}{Re} \left(G^{\prime \prime}+\frac {G^{\prime}}{r}- \left(k^2+\frac {n^2+1}{r^2}\right)G-\frac {2n}{r^2}H \right), \end{align}

(2.10b) \begin{align} k(U-c)H &=-\frac {n}{r}P-\frac {i}{Re}\left(H^{\prime \prime}+\frac {H^{\prime}}{r}-\left(k^2+\frac {n^2+1}{r^2}\right)H-\frac {2n}{r^2}G\right), \end{align}

(2.10c) \begin{align} k(U-c)F+U^{\prime}G &=-kP-\frac {i}{Re}\left(F^{\prime \prime}+\frac {F^{\prime}}{r}-\left(k^2+\frac {n^2}{r^2}\right)F\right), \end{align}

(2.10d) \begin{align} G^{\prime}+\frac {G}{r} &+\frac {n}{r}H+kF =0, \end{align}

where $Re$ is the Reynolds number. To solve the above eigenvalue problem, we employed a numerical code using the Chebyshev-collocation method originally developed by Deguchi & Nagata (2011) for the annular domain. For the pipe flow problem, the basis functions are replaced with those of Deguchi & Walton (2013). The codes have been thoroughly tested and validated against results from other research groups (see, e.g. Deguchi 2017; He et al. 2025).

3. Stability conditions

Throughout the paper, we assume that $U(r)$ possesses $C^1$ smoothness in $\varOmega$ (i.e. the first derivative is continuous in the domain). Here, for the annular case, $\varOmega =(r_i,r_o)$ , and for the pipe case, $\varOmega =(0,1)$ . For the pipe problem, we further assume that the even extension of $U(r)$ is $C^1$ in $(-1,1)$ .

To state the main result of this paper, it is convenient to define

(3.1) \begin{align} W_{\alpha ,N}(r)=\frac {rQ^{\prime}}{U-\alpha }. \end{align}

Batchelor & Gill (1962) also used a similar quantity in their analysis; however, note that our $Q$ is defined differently from theirs. The right-hand side of (3.1) depends on a constant $\alpha$ and the wavenumber ratio $N$ via $Q$ , hence the subscripts.

The main results of this paper are presented in the following theorems, which provide convenient criteria for assessing the inviscid stability problem (2.1).

Theorem 1. The base flow is stable if the KA-I and/or KA-II conditions are satisfied for all $N$ .

1. (i) KA-I: There exists $\alpha \in \mathbb{R}$ such that $W_{\alpha ,N}\leqslant 0$ for all $r \in \varOmega$ .

2. (ii) KA-II: There exists $\beta \in \mathbb{R}$ such that $W_{\beta ,N}=({rQ^{\prime}}/({U-\beta }))\in [0,H(r)]$ for all $r \in \varOmega$ . Here,

   (3.2) \begin{align} H=\left \{ \begin{array}{c} \left(\frac {\pi }{\ln \eta }\right)^2 \frac {1}{N^2+r_o^2}+\frac {1}{N^2}\qquad \text{if} \qquad N\geqslant 1,\\ \left(\frac {\pi }{\ln \eta }\right)^2 \frac {1}{r_o^2} \qquad \quad \qquad \text{if} \qquad N= 0 ,\end{array} \right . \end{align}
   for the annular case,
   (3.3) \begin{align} H=\left \{ \begin{array}{c} (j_{1,1}^2r^2+N^{-2})/(N^2+1)\qquad \text{if} \qquad N\geqslant 1,\\ j_{1,1}^2 \qquad \text{if} \qquad N= 0 ,\end{array} \right. \end{align}
   for the pipe case ( $j_{1,1}\approx 3.8317$ is the first positive zero of the Bessel function of the first kind of order one).

Theorem 2. The base flow is unstable if the following conditions are satisfied for some $N$ :

1. (i) $Q^{\prime}(r)$ has only one zero at $r=r_c$ and it is isolated.

2. (ii) $Q^{\prime \prime}(r_c)\neq 0$ .

3. (iii) $U^{\prime}(r_c)\neq 0$ .

4. (iv) $W_{\alpha ,N}\gt h(r)$ with $\alpha =U(r_c)$ for all $r\in \varOmega$ . Here,

   (3.4) \begin{align} h=\frac {(\pi /\ln \eta )^2}{N^2+ r_i^2}+\frac {1}{N^2} , \end{align}
   for the annular case,
   (3.5) \begin{align} h=Cr^p, \qquad p\geqslant 2, \qquad C=\frac {1}{\rho _p}\left(\rho _N+\frac {1}{6N^2}\right)\!,\nonumber \\ \rho _N=\frac {(3N^2+1)^2}{2}\ln \left(\frac {N^2+1}{N^2}\right)-\frac {3(6N^2+1)}{4},\nonumber \\ \rho _p=\frac {8}{(p+2)(p+4)(p+6)}, \end{align}
   for the pipe case.

Theorem 1 corresponds to KA-I and KA-II for axisymmetric base flows, and guarantees regions in parameter space where the base flow is stable. KA-I is due to Batchelor & Gill (1962), whereas KA-II is a new result. Note that when $Q^{\prime}$ has no zeros for any $N$ , the base flow is stable, since in this case one can choose $|\alpha |$ sufficiently large to fulfil the KA-I condition. This is an extension of Rayleigh’s inflection point theorem.

Theorem 2 serves as a handy means of demonstrating the existence of instability in the parameter space. It extends the hurdle theorem by Deguchi et al. (2024); instability emerges when $W_{\alpha ,N}$ surpasses the hurdle $h$ . Deguchi et al. (2024) defined the ‘reciprocal Rossby–Mach number’ so that the hurdle height equals unity, inspired by studies of Jupiter’s jets. An analogous quantity can be defined by multiplying $W_{\alpha ,N}$ by a constant factor. Readers interested in the interpretation of our stability conditions using that number are referred to Appendix A.4.

The proofs of above theorems are given in § 5.

4. Analysis of model profiles

4.1. Flow through an annulus

Consider the flow of a fluid between co-axial cylinders whose common axis is aligned with gravity (see figure 1 a). The fluid motion is driven by two effects: the differential axial movement of the cylinders and buoyancy resulting from an imposed temperature difference between the cylinders.

Figure 1. The annular model flow used in § 4.1. (a) Schematic of the model flow in the dimensional cylindrical coordinates $(r_*,\varphi ,z_*)$ . Here, the radii $r_{i*}$ and $r_{o*}$ satisfy $r_{o*}-r_{i*}=2d_*$ , $r_{i*}/r_{o*}=\eta$ . (b) The base flow $U(r)$ given in (4.5) for $\chi =-7,-1,1,7$ .

We take half of the cylinder gap $d_*$ , half of the differential speed $u_*$ and half of the imposed temperature difference $\Delta \theta _*$ as the characteristic length, velocity and temperature scales, respectively. The non-dimensional Navier–Stokes equations under the Boussinesq approximation yield the governing equations for the axial base velocity $U(r)$ and the base temperature $\varTheta (r)$

(4.1a) \begin{align} \varTheta^{\prime \prime}+r^{-1}\varTheta^{\prime} &=0, \end{align}

(4.1b) \begin{align} U^{\prime \prime}+r^{-1}U^{\prime} &=\chi \varTheta . \end{align}

The left-hand sides represent diffusivity, while the term proportional to $\chi$ corresponds to buoyancy. The ratio of the Grashof number $Gr$ to the Reynolds number $Re$ defines the parameter $\chi$ as follows:

(4.2) \begin{align} \chi =\frac {Gr}{Re},\qquad Gr=\frac {\gamma _*g_* d_*^3 \Delta \theta _*}{\nu _*^2},\qquad Re=\frac {u_*h_*}{\nu _*}. \end{align}

Here, $g_*$ denotes the gravitational acceleration, $\gamma _*$ the thermal expansion coefficient and $\nu _*$ the kinematic viscosity of the fluid. The boundary conditions are

(4.3a) \begin{align} U=1, \,\,\varTheta =-1\qquad \text{at}\qquad r=r_i, \end{align}

(4.3b) \begin{align} U=-1, \,\,\varTheta =1\qquad \text{at}\qquad r=r_o. \end{align}

The limit $\chi \rightarrow 0$ corresponds to sliding Couette flow (Gittler 1993; Deguchi & Nagata 2011), while for large $|\chi |$ , the system approaches natural convection in a vertical annulus (Choi & Korpela 1980; Yao & Rogers 1989; Kang et al. 2015; Wang & Chen 2022).

The solution of (4.1) with the boundary conditions (4.3) is given by

(4.4) \begin{align} \varTheta &= C_1 \ln r+C_2, \end{align}

(4.5) \begin{align} U &= C_3\ln r+C_4 +\frac {\chi }{4} r^2(C_1-C_2- C_1 \ln r ), \end{align}

where

(4.6) \begin{align} C_1=-\frac {2}{\ln \eta }, \quad C_2=-1- C_1\ln r_i, \quad C_3=\frac {C_5-C_6}{\ln \eta }, \quad C_4=C_5-C_3\ln r_i,\nonumber \\ C_5=1+ \chi \frac {r_i^2}{4}\{ C_1\ln r_i-C_1+C_2 \},\qquad C_6=-1+\chi \frac {r_o^2}{4}\{ C_1 \ln r_o-C_1+C_2 \}. \end{align}

The base flow depends on $\eta$ and $\chi$ . Figure 1(b) shows the profiles at $\eta =0.7$ for various values of $\chi$ .

The function $Q$ defined in (2.1) is readily found. Its derivative, required for computation of $W_{\alpha ,N}$ in (3.1), is given by

(4.7) \begin{align} Q^{\prime}=r\frac {N^2\chi (C_2+C_1\ln r)+2C_3+\chi r^2 \frac {C_1}{2}}{(N^2+r^2)^2}. \end{align}

If $N$ is sufficiently large, $Q^{\prime}$ varies monotonically and therefore has at most one zero. This is convenient for applying Theorem 2.

Figure 2 summarises the inviscid stability in the $\eta$ – $\chi$ parameter plane. The blue region indicates stability guaranteed by Theorem 1, while the red region indicates instability guaranteed by Theorem 2. As expected, the neutral curve obtained by numerically solving the differential equation (the solid black line) lies between the two regions.

Figure 2. Stability diagram of the annular model flow (figure 1 a) in the $\chi$ – $\eta$ plane. All physically possible wavenumbers are covered. The black solid line represents the neutral curve of the inviscid problem (2.1). Stability is guaranteed by Theorem 1 in the blue region, while unstable modes are found by Theorem 2 in the red region.

In the limit $\eta \rightarrow 1$ , the analysis becomes particularly simple, so we begin our discussion from this case. This limit, referred to as the narrow-gap limit, corresponds to $r_i,r_o\rightarrow \infty$ in our non-dimensionalisation where the gap width is set to 2. We introduce the new variable $y=r-r_m \in (-1,1)$ , where $r_m=(r_o+r_i)/2$ is the mid-gap. Taking the limit $N\rightarrow \infty$ and $\eta \rightarrow 1$ while keeping $N\ll r_m$ , we obtain

(4.8) \begin{align} W_{\alpha ,N}=-\frac {1}{U-\alpha }\frac {d^2U}{{\textrm{d}}y^2},\qquad H=h=\left(\frac {\pi }{2}\right)^2, \end{align}

from (3.1), (3.2) and (3.4). Theorem 2 reduces to the hurdle theorem for parallel flows derived in Deguchi et al. (2024). The limiting form of the base flow is

(4.9) \begin{align} U=\frac {\chi }{6}y(y^2-1)-y. \end{align}

The numerator of $W_{\alpha ,N}$ , i.e. $ {d^2U}/{{\textrm{d}}y^2}$ , vanishes only at $y=y_c=0$ . Setting $\alpha =U(y_c)=0$ , we have

(4.10) \begin{align} W_{\alpha ,N}=\frac {6\chi }{6+\chi (1-y^2)}. \end{align}

The results of applying Theorems 1 and 2 for the narrow-gap limit case, shown in figure 3, can be understood as follows. When $\chi \gt 0$ , the maximum of $W_{\alpha ,N}$ is $\chi$ and the minimum is $6\chi /(6+\chi )$ . Therefore, for $\chi \in (0,( {\pi }/{2})^2)$ , where $( {\pi }/{2})^2\approx 2.47$ , the KA-II condition of Theorem 1 is satisfied. Moreover, when $6\chi /(6+\chi )\gt ({\pi }/{2})^2$ , that is, for $\chi \gt {6\pi ^2}/{24-\pi ^2}\approx 4.19$ , Theorem 2 guarantees the existence of an instability. When $\chi \in (-6,0)$ , the KA-I condition in Theorem 1 is satisfied. For $\chi \leqslant -6$ , neither theorem can be applied, since the numerator of (4.10) vanishes at $y=\pm \sqrt {1+6/\chi }$ and $W_{\alpha ,N}$ is not continuous. By numerically solving the Rayleigh equation,

(4.11) \begin{align} \frac {d^2G}{{\textrm{d}}y^2}-k^2G-\frac {\frac {d^2U}{{\textrm{d}}y^2}}{U-c} G=0, \qquad G(-1)=G(1)=0, \end{align}

we find that the true neutral point occurs at $\chi \approx 3.810$ . Note that in the narrow-gap limit, Squire’s theorem implies that the stability is independent of the wavenumber ratio. Hence, there is no need to consider $N$ at all.

Figure 3. Stability diagram of the annular model flow at the narrow-gap limit $\eta \rightarrow 1$ . The eigenvalue problem (4.11) indicates the presence of unstable modes for $\chi \lt -6$ and $\chi \gt 3.81$ . The grey line shows that the profile of $W_{\alpha ,N}$ , given in (4.10), becomes singular when $\chi \leqslant -6$ .

Even for $\eta \lt 1$ , the behaviour of $W_{\alpha ,N}$ remains qualitatively similar to the narrow-gap case, as illustrated in figure 4 for $\eta =0.7$ . However, the results now depend on $N$ , since Squire’s theorem no longer applies. In the figure, we choose $N=10$ . Panel (a) illustrates the case $\chi =7$ , where $W_{\alpha ,N}$ surpasses the hurdle $h$ (defined in (3.4) and shown by the red line). The location $r=r_c$ indicates where $Q^{\prime}$ in (4.7) vanishes. When $\chi$ is reduced to 1, as shown in panel (b), the KA-II condition of Theorem 1 is satisfied; the blue line in the figure indicates $H$ in (3.2). Panel (c) is the case where $\chi =-1$ , for which $W_{\alpha ,N}$ is everywhere negative, and hence the flow is stable for $N=10$ perturbations by KA-I (i.e. the result by Batchelor & Gill 1962). Panel (d) corresponds to $\chi =-7$ , which lies below the lower boundary of the stable region (blue) shown in figure 2. In this case, $W_{\alpha ,N}$ cannot be made continuous for any choice of $\alpha$ , and neither Theorem 1 nor Theorem 2 can be used to determine stability.

Figure 4. Profiles of $W_{\alpha ,N}$ with $N=10$ for the annular model flow at $\eta =0.7$ . The constant $\alpha$ is set equal to $U(r_c)$ . Panels show (a) $\chi =7$ ; (b) $\chi =1$ ; (c) $\chi =-1$ ; (d) $\chi =-7$ . In panel (a), the red line shows $h$ from (3.4). In panel (b), the blue line shows $H$ defined in (3.2).

Figure 4 shows only the case for $N=10$ , but physically, perturbations with all values of $N$ can occur. Therefore, to use Theorem 1 to establish stability, it is not sufficient to examine only a handful of $N$ values. In figure 2, we verified the reliability of the blue stable region using a wide range of $N$ values. Note that the limit $N\rightarrow \infty$ is straightforward to analyse. This means that for sufficiently large $N$ , asymptotic behaviour emerges, so in practice it is sufficient to examine the conditions of the theorem up to moderately large $N$ . Likewise, at each point in the parameter plane, the value of $N$ is optimised to maximise the red unstable region. However, fixing $N=10$ yields almost the same outcome as that in figure 2. This suggests that identifying the unstable region using Theorem 2 is easier than finding the stable region using Theorem 1.

The neutral curve in figure 2 is also obtained by optimising over the wavenumbers $k$ and $n$ : above the black line an unstable mode is found for some $(k,n)$ , while below it no unstable mode exists for any choice of $(k,n)$ .

Figure 5 shows $c_i$ obtained from numerical computations of (2.1) for $\eta =0.7$ and $\chi =7$ . Panel (a) is computed with $N=10$ ; hence at $k=0.1$ , the mode with $n=Nk=1$ exhibits instability, as expected by Theorem 2. The neutral solution at $k=0.502$ is non-physical because $n=5.02$ is not an integer. Physically admissible neutral solutions can be obtained by reducing $N$ . For example, choosing $N=1/0.9475$ yields a neutral point at $k=k_0=0.9475$ , i.e. $n=1$ . Panel (b) shows the corresponding eigenfunction. This function is free of singularities due to the regularity of $W_{\alpha ,N}$ . The red dashed line in panel (a) shows an analytic linear approximation of how $c_i$ behaves near the neutral point. As will be shown in § 5.2 (see (5.23)), we can derive the identity

(4.12) \begin{align} \left . \frac {{\textrm{d}}c_i}{{\textrm{d}}k}\right |_{k=k_0} =-\frac {2k_0K_i}{K_r^2+K_i^2}\int _{\varOmega }rG_0^2{\textrm{d}}r, \end{align}

where

(4.13) \begin{align} K_r=\lim _{\epsilon \rightarrow 0}\left\{\int _{r_i}^{r_c-\epsilon }\frac {r^2Q^{\prime}G_0^2}{(U-\alpha )^2}{\textrm{d}}r+\int ^{r_o}_{r_c+\epsilon }\frac {r^2Q^{\prime}G_0^2}{(U-\alpha )^2}{\textrm{d}}r\right\}, \quad K_i=\pi \left . \left (\frac {rW_{\alpha ,N}G_0^2}{|U^{\prime}|} \right ) \right |_{r=r_c} \end{align}

can be found using the neutral wavenumber $k_0$ , the regular neutral eigenfunction $G_0$ and $\alpha =U(r_c)$ . The right-hand side of (4.12) is evaluated as $-0.886$ , which gives the slope of the red dashed line. As similar results hold for $N\approx 1/0.9475$ , one can conclude that a slight increase of $N$ from its neutral value necessarily results in the emergence of an unstable mode with $n=1$ . The observation here highlights the strategy to be used in § 5.2 for proving Theorem 2.

Figure 5. Inviscid stability result for the annular model flow at $(\eta ,\chi )=(0.7,7)$ . (a) Imaginary part of the phase speed $c_i$ for $N=10$ . The neutral point is at $k=k_0=0.502$ . The dashed red line indicates the result using (4.12). (b) Eigenfunction of the neutral mode found at $N=1/k$ , $k=0.9475$ (i.e. $n=1$ ).

Along the neutral curve in figure 2, the eigenfunction is regular. However, when $W_{\alpha ,N}$ is singular (as in figure 4 d), neutral solutions with singularities at the critical levels (i.e. locations where $U=c$ ) may arise. In the parameter plane shown in figure 2, such singular modes indeed occur when $\chi$ lies slightly below the lower boundary of the stable region. A detailed check in the narrow-gap limit confirmed that the neutral point lies very close to $\chi =-6$ (see figure 3). The computation was performed using 30 000 grid points in order to capture the nearly singular eigenfunctions with small $c_i$ . Likewise, the blue lower stability boundary in figure 2 may lie very close to the neutral curve, but we do not go into further detail. Note that such nearly singular unstable modes cannot be captured by Theorem 2 and are therefore outside the scope of the present study.

Figure 6 confirms the above inviscid stability results through computations of the linearised Navier–Stokes equation, (2.10). The dashed line indicates the neutral curve obtained by optimising over wavenumbers $k$ , $n$ . The most unstable mode is always $n=1$ , which asymptotes to the inviscid result $\chi =3.86$ as $Re$ increases. The stability for $Re\gt 10^4$ is relatively well captured by Theorems 1 and 2.

Figure 6. Comparison between the viscous and inviscid stability analyses for the annular model flow at $\eta =0.7$ . The dashed line represents the neutral curve obtained from the viscous stability problem (2.10), covering all physically possible wavenumbers. The blue, black and red vertical lines correspond to the inviscid stability results shown in figure 2.

4.2. Flow through a pipe

In this section, we consider vertically oriented Hagen–Poiseuille flow with homogeneous internal heating (figure 7 a, see Senoo et al. 2012; Marensi et al. 2021 also). We choose the pipe radius $d_*$ as the length scale, and the centreline velocity of the laminar Hagen–Poiseuille flow $u_*$ as the velocity scale. For the temperature scale, we use the difference $\Delta \theta _*$ between the base temperature at the centreline, $\theta _{0*}+\Delta \theta _*$ , and that at the wall, $\theta _{0*}$ . The velocity and temperature scales are related to the dimensional axial pressure gradient $\Delta p_*/L_*$ and the dimensional internal heat source $q_*$ as follows:

(4.14) \begin{align} u_*=\frac {d_*^2}{4\mu _*}\frac {\Delta p_*}{L_*},\qquad \theta _*=\frac {d_*^2q_*}{4\kappa _*}. \end{align}

Here, $\mu _*$ is the dynamic viscosity and $\kappa _*$ is the thermal diffusivity of the fluid.

Figure 7. The model flow used in § 4.2. (a) Schematic of the model flow in the dimensional cylindrical coordinates $(r_*,\varphi ,z_*)$ . (b) The base flow $U(r)$ given in (4.19) for $\chi =-10,5,7,14$ .

From the Boussinesq-approximated Navier–Stokes equations, the equations for the axial base velocity $U(r)$ and the base temperature $\varTheta (r)$ are obtained as

(4.15) \begin{align} \varTheta^{\prime \prime}+r^{-1}\varTheta^{\prime} &=-4, \end{align}

(4.16) \begin{align} U^{\prime \prime}+r^{-1}U^{\prime} &=\chi \varTheta -4. \end{align}

The solutions of the above equations, satisfying the boundary conditions

(4.17) \begin{align} U=0, \,\,\varTheta =0\qquad \text{at}\qquad r=1, \end{align}

and the centreline regularity, are

(4.18) \begin{align} \varTheta &=(1-r^2), \end{align}

(4.19) \begin{align} U &=(1-r^2)-\frac {\chi }{16}(1-r^2)(3-r^2). \end{align}

The base flow profiles for selected values of $\chi$ are shown in figure 7(b). An inflection point appears when $\chi$ is in a certain range, and by analogy with Rayleigh’s theorem, Senoo et al. (2012) speculated that this could be a possible cause of instability. However, the mathematics of inviscid instability is not that simple.

Substituting the base flow $U$ into the definition of $Q$ (see (2.1)) and differentiating, we obtain

(4.20) \begin{align} Q^{\prime}=r\frac {N^2(4-\chi (1-r^2))+\frac {\chi }{2}r^4}{(N^2+r^2)^2}. \end{align}

It is easy to see that $Q^{\prime}$ vanishes at $r=r_c$ , where

(4.21) \begin{align} r_c=N\sqrt {\sqrt {1+\frac {2}{N^2}\big(1-\frac {4}{\chi }\big)}-1 }. \end{align}

Clearly, the axisymmetric mode is stable, since for $N=0$ , the function $Q^{\prime}$ does not change sign in the domain (the KA-I condition is satisfied, see the remark in § 3). Even for non-axisymmetric modes, the parameter range $\chi \in (0,4)$ remains stable, as no real $r_c$ exists. Determining inviscid stability for other values of $\chi$ is not trivial, but by applying Theorems 1 and 2, we can find the results shown in figure 8.

Figure 8. Stability diagram of the pipe model flow (figure 7 a). The eigenvalue problem (2.1) indicates the presence of unstable modes for $\chi \in (4,7.86)$ . The grey line shows that the profile $W_{\alpha ,N}$ becomes singular when $\chi \in [4,6]$ .

When $\chi$ is negative, either KA-I or KA-II holds. Figure 9(a) shows the profile of $W_{\alpha ,N}$ for $\chi =-10$ and $N=1$ . This is an example in which the KA-II condition defined in (3.3) is satisfied; $W_{\alpha ,N}$ is positive everywhere and lies below $H$ indicated by the blue line.

Figure 9. Profiles of $W_{\alpha ,N}$ with $\alpha =U(r_c)$ and $N=1$ for the pipe model flow. Panels show (a) $\chi =-10$ ; (b) $\chi =5$ ; (c) $\chi =7$ ; (d) $\chi =14$ . In panel (c), the red line shows $h$ from (3.5). In panels (a) and (d), the blue line shows $H$ defined in (3.3).

For $\chi \in [5.72,7.42]$ , there exist $N\gt 0$ and $p\geqslant 2$ such that $W_{\alpha ,N}$ exceeds the hurdle (3.5). Figure 9(c) shows the profile of $W_{\alpha ,N}$ for $\chi =7$ and $N=1$ . In this figure, the hurdle for $p=2$ is indicated by the red line. Thus, by Theorem 2, the flow is unstable. Note that a horizontal hurdle, as used in the annular case, does not work well for the pipe flow since near $r=0$ , $W_{\alpha ,N}$ behaves like $r^2$ . Figure 10(a) shows the eigenvalue analysis of (2.1) for $\chi =7$ and $N=1$ . This calculation uses method (ii) of § 2.2, with the condition $G^{\prime}(\epsilon )=0$ . At $k=1$ , we indeed have an unstable mode with $n=1$ .

Figure 10. Inviscid stability result for the pipe model flow at $\chi =7$ . (a) Imaginary part of the phase speed $c_i$ for $N=1$ . The neutral point is at $k=k_0=1.159$ . The dashed red line indicates the result using (4.12). (b) Eigenfunction of the neutral mode found at $N=1/k$ , $k=1.46$ (i.e. $n=1$ ).

The red dashed line in figure 10(a) corresponds to results similar to those seen in figure 5(a). The slope $-0.0549$ is obtained by the neutral mode using (4.12). By decreasing the value of $N$ from 1 to $1/k$ with $k=1.46$ , a neutral mode with $n=1$ is obtained. This neutral mode is regular as shown in figure 10(b). We confirmed that the methods (i) and (ii) introduced in § 2.2 produce eigenfunctions that are graphically indistinguishable.

By computing neutral points using (2.1) for all physically admissible wavenumbers, it is found that a neutral point occurs at $\chi =7.86$ (see figure 8). This value is reasonably close to the unstable region predicted by Theorem 2. The eigenvalue problem (2.1) does not produce unstable modes for $\chi \gt 7.86$ . This stabilisation can be detected by KA-II when $\chi \gt 13.18$ . For example, figure 9(d) shows the case $\chi =14$ , $N=1$ , where $W_{\alpha ,N}$ lies below the blue line $H$ determined by (3.3).

Another neutral point is expected to exist between the blue stability region and the red instability region in figure 8. The computations of (2.1) indicates that instability sets in when $\chi$ slightly exceeds 4, implying that the stability boundary $\chi =4$ by Theorem 1 sharply predicts the neutral point. The eigenfunction corresponding to this neutral point is singular. Hence, it makes sense that the boundary of the unstable region by Theorem 2 is not sharp. Careful observation of $W_{\alpha ,N}$ shows that it becomes singular for some $N$ when $\chi \in [4,6]$ . An example of such a singular case is shown in figure 9(b), obtained for $\chi =5$ and $N=1$ .

Figure 11 verifies the above inviscid results using viscous computations. The dashed line represents the envelope of the neutral curves obtained from the linearised Navier–Stokes equations (2.10) over a range of wavenumbers. This curve is determined by the $n=1$ mode and asymptotes to the inviscid neutral points at $\chi =4$ and $7.86$ for large Reynolds numbers. Recall that for $\chi \in [4,6]$ , instability may arise from the singular neutral mode, whereas for $\chi \in [5.72,7.42]$ , Theorem 2 guarantees the existence of regular neutral modes. The neutral curve found by the viscous analysis features two humps, representing two unstable modes, which correspond to a singular mode and a regular mode in the inviscid analysis.

Figure 11. Comparison between the viscous and inviscid stability analyses for the pipe model flow. The dashed line represents the neutral curve obtained by (2.10) varying wavenumbers. The blue, black and red vertical lines correspond to the inviscid stability results shown in figure 8.

5. Mathematical proofs

Here, we prove Theorems 1 and 2 stated in § 3. For definiteness, we introduce the function space $H_0^1$ , following the standard notation in functional analysis. For $f \in H_0^1$ , both $\int _{\varOmega } (f^{\prime})^2 r,{\textrm{d}}r$ and $\int _{\varOmega } f^2 r,{\textrm{d}}r$ are finite, and the boundary conditions are satisfied. Specifically, in the annular case, $f(r_i)=f(r_o)=0$ , while in the pipe case, $f(1)=0$ . Recall that the pipe solution also satisfy the centreline regularity condition (2.6). If $G$ satisfies (2.1) in the classical sense, then $G \in H_0^1$ .

5.1. Proof of Theorem 1

Following Batchelor & Gill (1962), we multiply (2.1) by the complex conjugate of $rG$ and integrate over the entire domain $\varOmega$ . The real and imaginary parts of the integral are obtained as

(5.1) \begin{align} \int _{\varOmega } \left\{\frac {r}{N^2+ r^2}|(rG)^{\prime}|^2 +\frac {k^2}{r}|rG|^2-\frac {Q^{\prime}(U-c_r)}{|U-c|^2}|rG|^2 \right\}{\textrm{d}}r & =0, \end{align}

(5.2) \begin{align} c_i\int _{\varOmega } \frac {Q^{\prime}}{|U-c|^2}|rG|^2{\textrm{d}}r &=0, \end{align}

respectively. The boundary terms arising from integration by parts vanish in all cases we consider.

Now, we assume $c_i\neq 0$ and show that the KA-I assumption leads to a contradiction. This can be found in the same manner as Fjortoft’s analysis for parallel flows. By taking a suitable linear combination of the two integrals (5.1) and (5.2), we obtain

(5.3) \begin{align} \int _{\varOmega } \left\{\frac {r}{N^2+r^2}|(rG)^{\prime}|^2 +\frac {k^2}{r}|rG|^2-W_{\alpha ,N}\frac {(U-\alpha )^2}{|U-c|^2}\frac {|rG|^2}{r} \right\}{\textrm{d}}r=0, \end{align}

where $\alpha \in \mathbb{R}$ is arbitrary. If $W_{\alpha ,N}=rQ^{\prime}/(U-\alpha )\leqslant 0$ for all $r \in \varOmega$ , the above integral cannot be satisfied for non-trivial $G$ . The analysis here is identical to that in Batchelor & Gill (1962) but included for completeness.

KA-II stability is usually derived by defining a Hamiltonian. However, in our case, an equivalent condition can be obtained in a more straightforward manner. We first set $\alpha =2c_r-\beta$ to rewrite (5.3) in the form

(5.4) \begin{align} \int _{\varOmega } \left\{\frac {r}{N^2+ r^2}|(rG)^{\prime}|^2 +\frac {k^2}{r}|rG|^2-W_{\beta ,N}Z\frac {|rG|^2}{r} \right\}{\textrm{d}}r=0, \end{align}

where $Z= ({(U-c_r)^2-(c_r-\beta )^2})/({(U-c_r)^2+c_i^2})\lt 1$ . The choice of $\beta \in \mathbb{R}$ is arbitrary; by assumption, it can be adjusted so that $W_{\beta ,N}\in [0,H(r))$ for $r\in \varOmega$ . From this point onward, the analysis must be treated separately for annular and pipe flows.

For the annular problem, there is a positive constant $\kappa _N^2$ depending on $N$ such that

(5.5) \begin{align} \int _{r_i}^{r_o} \frac {r}{N^2+ r^2}|(rf)^{\prime}|^2 {\textrm{d}}r\geqslant \kappa _N^2 \int _{r_i}^{r_o} r^{-1}|rf|^2{\textrm{d}}r, \end{align}

for all $f(r) \in H_0^1$ ; this is a generalisation of Poincaré’s inequality. Using the non-negativity of $W_{\beta ,N}$ and (5.5) in (5.4), we obtain

(5.6) \begin{align} 0\gt \int _{r_i}^{r_o} \left (\kappa _N^2 +k^2-W_{\beta ,N} \right ) \frac {|rG|^2}{r} {\textrm{d}}r. \end{align}

Computation in Appendix A.1 shows that

(5.7) \begin{align} \kappa _N^2=\frac {(\pi /\ln \eta )^2}{N^2+r_o^2}, \end{align}

and with this expression, (3.2) becomes

(5.8) \begin{align} H=\left \{ \begin{array}{c} \kappa _N^2+\frac {1}{N^2}\qquad \text{if} \qquad N\geqslant 1,\\ \kappa _0^2 \qquad \qquad \, \text{if} \qquad N= 0. \end{array} \right . \end{align}

If $n=0$ (i.e. $N=0$ ) and $W_{\beta ,0}\lt \kappa _0^2$ , clearly the inequality (5.6) cannot be satisfied by non-trivial $G$ . When $n\neq 0$ (i.e. $N\gt 0$ ), we can use $k^2N^2=n^2\geqslant 1$ . Thus, in this case, when $W_{\beta ,N}\lt \kappa _N^2+({N^{-2}})$ the inequality (5.6) cannot be satisfied.

For the pipe problem, different inequalities must be applied separately to the axisymmetric case ( $n=0$ ) and the non-axisymmetric cases ( $n\geqslant 1$ ). If $n=0$ , we use the fact that

(5.9) \begin{align} \int _{0}^{1} \frac {1}{r}|(rf)^{\prime}|^2 {\textrm{d}}r\geqslant j_{1,1}^2 \int _{0}^{1} \frac {1}{r}|rf|^2{\textrm{d}}r, \end{align}

for all $f\in H_0^1$ satisfying the regularity condition (see Appendix A.2). Then (5.9) and the integral (5.4) imply that

(5.10) \begin{align} 0\gt \int _0^1 \left (j_{1,1}^2 +k^2-W_{\beta ,0} \right ) \frac {|rG|^2}{r} {\textrm{d}}r. \end{align}

Clearly this does not happen when $W_{\beta ,0}\leqslant j_{1,1}^2$ .

If $n\geqslant 1$ , we can use the following inequality, which holds for all regular $f\in H_0^1$ (see Appendix A.3):

(5.11) \begin{align} \int ^1_0 \left \{r|(rf)^{\prime}|^2 +\frac {n^2}{r}|rf|^2 \right\}{\textrm{d}}r\geqslant j_{1,1}^2 \int ^1_0 r|rf|^2{\textrm{d}}r. \end{align}

The integral (5.4) then yields the inequality

(5.12) \begin{align} 0\gt & \int ^1_0 \left \{\frac {r}{N^2+ 1}|(rG)^{\prime}|^2 +\frac {n^2}{rN^2}|rG|^2-W_{\beta ,N}\frac {|rG|^2}{r} \right\}{\textrm{d}}r\nonumber \\ & = \int ^1_0 \left \{\frac {r}{N^2+ 1}|(rG)^{\prime}|^2 +\frac {n^2}{r(N^2+1)}|rG|^2 \right\}{\textrm{d}}r\nonumber \\ & +\int ^1_0 \left \{ \frac {n^2}{rN^2(N^2+1)}|rG|^2-W_{\beta ,N}\frac {|rG|^2}{r} \right\}{\textrm{d}}r\nonumber \\ & \geqslant \int ^1_0 \left \{ \frac {j_{1,1}^2r^2}{N^2+ 1} + \frac {1}{N^2(N^2+1)}-W_{\beta ,N} \right\}\frac {|rG|^2}{r} {\textrm{d}}r, \end{align}

which is impossible when $W_{\beta ,N}\leqslant (j_{1,1}^2 r^2+N^{-2})/(N^2+1)$ for all $r \in (0,1)$ .

5.2. Proof of Theorem 2: step (i)

To carry out step (i) described in § 1, we first note that a neutral solution satisfies the integral equation

(5.13) \begin{align} \frac {\int _{\varOmega } \left \{ \frac {r}{N^2+ r^2}|(rG)^{\prime}|^2 -W_{c,N}\frac {1}{r}|rG|^2 \right\} {\textrm{d}}r}{\int _{\varOmega } \frac {1}{r}|rG|^2 {\textrm{d}}r}=-k^2 , \end{align}

with $c \in \mathbb{R}$ . This observation motivates us to define a functional $R_{\alpha ,N}: H_0^1 \rightarrow \mathbb{R}$ depending on $\alpha , N\in \mathbb{R}$ by

(5.14) \begin{align} R_{\alpha ,N}(\phi )=\frac {\int _{\varOmega } \left\{ \frac {r}{N^2+ r^2}|(r\phi )^{\prime}|^2 -W_{\alpha ,N}\frac {1}{r}|r\phi |^2 \right\} {\textrm{d}}r}{\int _{\varOmega } \frac {1}{r}|r\phi |^2 {\textrm{d}}r}. \end{align}

Then, (5.13) can be compactly written as $R_{c,N}(G)=-k^2$ . A function that satisfies this equation with real $k$ and $c$ may correspond to a neutral solution, but at this stage, it is not obvious whether such a solution exists.

The key mathematical fact we use is that when $W_{\alpha ,N}$ is continuous, the minimum of the Rayleigh quotient (5.14)

(5.15) \begin{align} \lambda _0=\min _{\phi \in H_0^1} R_{\alpha ,N}(\phi ), \end{align}

exists, thanks to the completeness of $H_0^1$ . The minimiser $\phi _0$ satisfies the Euler–Lagrange equation

(5.16) \begin{align} \left \{\frac {r}{N^2+r^2}(r\phi )^{\prime} \right\}^{\prime}+(\lambda +W_{\alpha ,N})\phi =0, \end{align}

with $\lambda =\lambda _0$ and $\phi =\phi _0$ . This equation has the same form as (2.1). Hence, by setting $G=\phi _0$ , $c=\alpha$ , and $k=\sqrt {-\lambda _0}$ , a neutral solution is constructed, although physically $\lambda _0$ must be negative at least. The remaining task is to determine the condition under which $\lambda _0$ becomes sufficiently small. For this purpose, we will introduce appropriate trial functions in § 5.4.

Here, we remark that, for the annular domain, (5.16) together with the boundary conditions forms a regular Sturm–Liouville problem, which, as is well known, admits countably many real eigenvalues that can be ordered as $\lambda _0 \lt \lambda _1 \lt \ldots \lt \lambda _m \lt \ldots$ , with $\lambda _m \to \infty$ as $m \to \infty$ . The minimum of the Rayleigh quotient corresponds to the smallest eigenvalue $\lambda _0$ . The associated minimiser $\phi _0$ is real and has no zeros in $\varOmega$ . This latter property follows from Sturm’s oscillation theorem and plays an important role in step (ii).

For the pipe case, however, (5.16) is a singular Sturm–Liouville problem; readers interested in the details of regular and singular Sturm–Liouville problems are referred to the classical textbooks by Titchmarsh (1962), Courant & Hilbert (1989) and Birkhoff & Rota (1989) or the more modern treatment by Zettl (2005). Singular cases can be further classified; our pipe problem belongs to the limit-point non-oscillatory case, where the eigenfunction behaves relatively well near $r=0$ .

That said, there is no need to distinguish between the annular and pipe cases in the proof, except when constructing the trial function that determines $h$ in the theorem in § 5.4. Recall that in method (ii) of § 2.2, we considered an annular problem with a thin virtual inner cylinder of radius $\epsilon \gt 0$ . For this augmented pipe problem, the discussions in §§ 5.2 and 5.3 remain applicable. We then expect that the neutral or unstable solutions constructed there will converge to those of the pipe problem as $\epsilon \rightarrow 0$ . We have already demonstrated numerically that this expectation holds in § 4.2; for a mathematical justification, see the textbook by Titchmarsh (1962), for example.

5.3. Proof of Theorem 2: step (ii)

For both the annular and pipe cases, the assumptions of Theorem 2 are necessary to guarantee the existence of a neutral solution using specific trial functions $\psi$ . The proof is deferred to § 5.4, while step (ii) is presented first.

Throughout this section, we assume that the domain is an annulus. Once the existence of a neutral mode is established, we can then show that unstable modes appear when $k$ is slightly decreased. Let us fix $N$ and vary $\lambda = -k^2$ in (2.1). Both $G$ and $c$ depend on $\lambda$ , and their evaluation at the neutral parameter $\lambda =\lambda _0$ will be denoted by the subscript $0$ . Evaluation of (2.1) at this parameter yields $\mathcal{L} G_0 = 0$ , where $\mathcal{L}$ denotes the operator on the left-hand side of (5.16) with $\lambda =\lambda _0$ . By the properties of regular Sturm–Liouville problems, $G_0$ is real valued.

The eigenvalue $c$ of (2.1) can be expanded in a Taylor series. Thus if the first-order coefficient of the expansion

(5.17) \begin{align} c_i(\lambda )=\left . \frac {{\textrm{d}}c_i}{{\textrm{d}}\lambda }\right |_0(\lambda -\lambda _0)+O((\lambda -\lambda _0)^2) \end{align}

is non-zero, this implies the emergence of an unstable mode. The eigenvalue $c$ behaves well in the complex plane except along the real axis since the two linearly independent solutions themselves remain well behaved; see Lin (2003) for a mathematically rigorous discussion.

To proceed, we first differentiate (2.1) with respect to $\lambda$ to obtain

(5.18) \begin{align} \left(\frac {r}{N^2+r^2}(rG_{\lambda })^{\prime}\right)^{\prime}+\left(\lambda +\frac {rQ^{\prime}}{U-c}\right)G_{\lambda }+\left(1 +\left(\frac {rQ^{\prime}}{U-c}\right)_{\lambda }\right)G=0. \end{align}

Here, the subscript $\lambda$ denotes partial differentiation. Evaluating this equation at $\lambda =\lambda _0$ ,

(5.19) \begin{align} \mathcal{L}G_{\lambda }|_0+\left . \left (1 +\left (\frac {rQ^{\prime}}{U-c}\right )_{\lambda }\right |_0 \right )G_0=0. \end{align}

Multiplying (5.19) by $rG_0$ , subtracting $rG_{\lambda }|_0\mathcal{L}G_0=0$ , and integrating over the domain, we obtain

(5.20) \begin{align} 0=\int _{\varOmega }(G_0\mathcal{L}G_{\lambda }|_0-G_{\lambda }|_0\mathcal{L}G_0)r{\textrm{d}}r+\int _{\varOmega }\left . \left (1 +\left (\frac {rQ^{\prime}}{U-c}\right )_{\lambda }\right |_0 \right )G_0^2r{\textrm{d}}r. \end{align}

The first integral vanishes after integration by parts, so we get

(5.21) \begin{align} 0=\int _{\varOmega }G_0^2r{\textrm{d}}r+\left . K\frac {{\textrm{d}}c}{{\textrm{d}}\lambda } \right |_0, \end{align}

where

(5.22) \begin{align} K & =K_r+\textit{iK}_i=\lim _{c_i\rightarrow 0^+}\int _{\varOmega }\frac {rQ^{\prime}G_0^2}{(U-\alpha -ic_i)^2}r{\textrm{d}}r\nonumber \\ & = {\int \kern-12pt-}_{\varOmega } \frac {rQ^{\prime}G_0^2}{(U-\alpha )^2}r{\textrm{d}}r+i\pi \left . \left (\frac {rW_{\alpha ,N}G_0^2}{|U^{\prime}|} \right ) \right |_{r=r_c}. \end{align}

The dashed integral is the Cauchy principal integral. A heuristic derivation of (5.22) is well known in the shear flow instability community (see Drazin & Reid 1981, for example). A rigorous mathematical justification, along the lines of Kumar & Ożański (2025), can be straightforwardly done using the Plemelj–Sochocki formula.

Combining (5.21) and (5.22), we obtain

(5.23) \begin{align} \left . \frac {{\textrm{d}}c_i}{{\textrm{d}}k}\right |_0 = -2\sqrt {-\lambda _0}\left . \frac {{\textrm{d}}c_i}{{\textrm{d}}\lambda }\right |_0 =-\frac {2\sqrt {-\lambda _0}K_i}{K_r^2+K_i^2}\int _{\varOmega }G_0^2r{\textrm{d}}r. \end{align}

A slight manipulation yields (4.12). As remarked earlier, $G_0$ has no zeros in the domain, so $\int _{\varOmega }G_0^2r{\textrm{d}}r\gt 0$ . Furthermore, from the assumptions, $W_{\alpha ,N}|_{r=r_c}= {r_cQ^{\prime \prime}(r_c)}/{U^{\prime}(r_c)}$ is non-zero. In fact, if $W_{\alpha ,N}$ surpasses the hurdle, then $W_{\alpha ,N}|_{r=r_c}$ is strictly positive, and so is $({{\textrm{d}}c_i}/{{\textrm{d}}\lambda }) |_0$ . Therefore, if $k$ is slightly decreased from its neutral value, $c_i$ must increase, as observed in figures 5(a) and 10(a).

5.4. Proof of Theorem 2: trial functions

To guarantee the existence of a physically relevant neutral solution, let us first consider how small $\lambda _0$ in (5.15) needs to be. The discussion differs between axisymmetric and non-axisymmetric perturbations.

The axisymmetric case is very similar to the corresponding analysis for parallel flows; a neutral mode exists if there is a trial function $\psi \in H^1_0$ such that

(5.24) \begin{align} R_{\alpha ,0}(\psi )\lt 0. \end{align}

If this condition holds, we have $\lambda _0 \leqslant R_{\alpha ,0}(\psi )\lt 0$ , allowing us to set $k=\sqrt {-\lambda _0}$ .

For the non-axisymmetric case, a stronger condition is required: a neutral mode exists if there is a trial function $\psi \in H^1_0$ such that

(5.25) \begin{align} R_{\alpha ,N}(\psi )\lt -N^{-2} \end{align}

for some $N$ , say $N=N_+$ . This condition guarantees $\lambda _0\lt -N^{-2}$ at $N=N_+$ , so setting $k=\sqrt {-\lambda _0}$ yields a neutral solution with $N^{-1}\lt k$ . Note that this neutral solution may not be the one we want, since physically $n=Nk$ must be an integer. Nevertheless, condition (5.25) is sufficient to ensure the existence of a physically relevant neutral solution. The reason is as follows. We first note that the eigenvalue $\lambda _0$ of the regular Sturm–Liouville problem depends continuously on $N$ . Thus, if for some $N=N_-$ we have $\lambda _0\gt -N^2$ , then there is an $N \in (N_-,N_+)$ such that $\lambda _0=-N^{-2}$ from the mean value theorem. We can indeed find such $N_-$ at sufficiently small $N$ , since $-N^{-2}\rightarrow -\infty$ as $N\rightarrow 0$ , while $\lambda _0$ remains finite. By taking $k=\sqrt {-\lambda _0}$ at the value of $N$ where $\lambda _0=-N^{-2}$ is satisfied, we obtain a neutral mode with $n=1$ .

We now turn our attention to showing that the conditions in the theorem imply the existence of a neutral mode. We first prove the annular case. By assumption there is a constant $h\in \mathbb{R}$ such that $W_{\alpha ,N}\gt h$ for all $r \in \varOmega$ . Then using (5.14) we can show that

(5.26) \begin{align} R_{\alpha ,N}(\psi )\lt \frac {\int _{\varOmega } \frac {r}{N^2+ r^2}|(r\psi )^{\prime}|^2{\textrm{d}}r}{\int _{\varOmega } \frac {1}{r}|r\psi |^2 {\textrm{d}}r}-h\leqslant \frac {\frac {1}{N^2+ r_i^2}\int _{\varOmega } r|(r\psi )^{\prime}|^2{\textrm{d}}r}{\int _{\varOmega } \frac {1}{r}|r\psi |^2 {\textrm{d}}r}-h. \end{align}

The right-hand side becomes $( {(\pi /\ln \eta )^2}/({N^2+ r_i^2}))-h$ upon substituting the trial function $\psi (r)=r^{-1}\sin (\pi ( {\ln (r_i/r)}/{\ln \eta }) )$ . This constant equals $-N^2$ when $h$ is set as (3.4), in which case condition (5.25) is satisfied.

The proof of the pipe case is similar, but employs a different trial function. The assumption $W_{\alpha ,N}\gt Cr^{p}$ for all $r\in \varOmega$ implies

(5.27) \begin{align} R_{\alpha ,N}(\psi )\lt \frac {\int _{\varOmega }\left \{ \frac {r}{N^2+ r^2}|(r\psi )^{\prime}|^2 -Cr^p\frac {1}{r}|r\psi |^2 \right\} {\textrm{d}}r}{\int _{\varOmega } \frac {1}{r}|r\psi |^2 {\textrm{d}}r}. \end{align}

With the choice of trial function $\psi =(1-r^2)$ , which is suitable for $n=1$ modes, the right-hand side of (5.27) evaluates to $6(\rho _N-h\rho _p)$ . Substituting $C$ in (3.5), this constant becomes $-N^{-2}$ , which shows that (5.25) is satisfied.

6. Conclusion and discussion

We studied the inviscid stability problem (2.1) for axisymmetric base flows in annuli and pipes. The main results are the two theorems presented in § 3, which provide simple analytical tools to estimate the behaviour of the eigenvalue $c$ . By using those theorems complementarily, we can estimate the location of the neutral point in the parameter space, as demonstrated through the analysis of model flows.

Theorem 1 provides a sufficient condition for stability. The KA-II condition further extends the range of parameters for which stability can be guaranteed, beyond the classical KA-I result of Batchelor & Gill (1962). As shown in § 5.1, the KA-II condition follows from elementary algebraic manipulations together with the application of Poincaré-type inequality. In practice, the condition can be assessed by checking whether the function $W_{\alpha ,N}$ defined in (3.1) remains below a threshold $H$ (see figures 4 b and 9 a). The form of $H$ differs between annular and pipe flows and also depends on whether the disturbance is axisymmetric or non-axisymmetric.

Theorem 2 gives a sufficient condition for instability. For annular flows, a constant hurdle $h$ can be introduced, as in Deguchi et al. (2024) for parallel flows. Surpassing the hurdle $h$ by $W_{\alpha ,N}$ implies instability (figure 4 a). However, for pipe flows, the constant hurdle is ineffective as $W_{\alpha ,N}$ behaves like $r^2$ near the origin. Accordingly, we constructed the instability condition so that $h$ is modified to take a power-law dependence on $r$ (figure 9 c). As shown in the proof in § 5.2, the specific form of $h$ depends on the choice of trial function. While various forms of $h$ are possible, our choice offers an optimal balance between simplicity and effectiveness.

In step (i) of the proof of Theorem 2 (§ 5.2), Sturm–Liouville theory plays a key role in establishing that a neutral solution exists when the minimum of the Rayleigh quotient (5.15) is sufficiently small. In step (ii) (§ 5.3), a perturbation analysis is then used to estimate the change in the growth rate caused by a slight variation in the wavenumber and to identify the onset of instability (the red dashed lines in figures 5 a and 10 a). In § 5.4, test functions are selected to satisfy (5.24) for the axisymmetric case and (5.25) for the non-axisymmetric case.

The theoretical results have been applied to the annular model flow (§ 4.1) and the pipe model flow (§ 4.2). In both cases, the stable and unstable regions in parameter space predicted by Theorems 1 and 2, respectively, are consistent with the exact stability analyses obtained from the eigenvalue problem (2.1). Furthermore, these results agree with stability analyses based on the linearised Navier–Stokes equations in the limit of sufficiently large Reynolds numbers. The theorems require only a simple graphical analysis, which provides a rough estimate of stability in parameter space and thereby reduces the number of eigenvalue computations needed. They are not limited to model flows and could be applied in a variety of situations. For example, they can be used for local temporal stability analyses of slowly spatially developing flows, such as diverging pipe flows (Sahu & Govindarajan 2005).

It should be noted that the results of the theorems are valid only for stability analyses based on the inviscid approximation. As exemplified by Tollmien–Schlichting waves, viscosity can contribute to instability, so even within the stable regions predicted by Theorem 1, the full stability problem without approximation may exhibit instability. One approach to obtain an analytic estimate of instabilities that persist in the presence of viscosity is to use a Wentzel–Kramers–Brillouin type approximation (Kirillov, Stefani & Fukumoto 2014; Kirillov & Mutabazi 2025). However, this method is valid only for short wavelengths, where shear flows are typically stable, and hence does not capture the instabilities predicted by Theorem 2.

Finally, we discuss the stability of axisymmetric jets, which are not covered by our theorems but were the main subject of the study by Batchelor & Gill (1962). Unbounded domains like jets correspond to a singular Sturm–Liouville problem, which is likely to be more difficult than the pipe case and would require additional steps in the proof. Even ignoring this mathematical technicality, there are more fundamental difficulties. First, Poincaré-type inequalities cannot be applied in unbounded domains. Therefore, in Theorem 1, it is difficult to improve the stability results. Also, the form of $h$ given in Theorem 2 is not optimal for detecting unstable modes, as $W_{\alpha ,N}$ tends to zero in the far field. This is illustrated, for example, from figure 3 of Batchelor & Gill (1962), where $N^2$ corresponds to $W_{\alpha ,N}$ in the Schlichting jet

(6.1) \begin{align} U(r)=\frac {1}{(1+r^2)^2}. \end{align}

That said, there is a convenient way to design suitable trial functions for the jet flows, as shown in Appendix B. For the Schlichting jet, the integral in condition (5.25) can be calculated analytically, showing that instability occurs, which is consistent with the Navier–Stokes computations of Lessen & Singh (1973). However, in general, numerical integration is required to use (5.25). Whether simple conditions like those we derived for annular and pipe flows can be obtained for jets remains a topic for future work. If successful, this could facilitate the stability analysis of more practically relevant flows, such as plumes (Chakravarthy, Lesshafft & Huerre 2015).