Instabilities in visco-thermodiffusive swirling flows

Abstract

An analytical theory is presented for linear, local, short-wavelength instabilities in swirling flows, in which axial shear, differential rotation, radial thermal stratification, viscosity and thermal diffusivity are all taken into account. A geometrical optics approach is applied to the Navier–Stokes equations, coupled with the energy equation, leading to a set of amplitude transport equations. From these, a dispersion relation is derived, capturing two distinct types of instability: a stationary centrifugal instability and an oscillatory, visco-diffusive McIntyre instability. Instability regions corresponding to different axial or azimuthal wavenumbers are found to possess envelopes in the plane of physical parameters, which are explicitly determined using the discriminants of polynomials. As these envelopes are shown to bound the union of instability regions associated with particular wavenumbers, it is concluded that the envelopes correspond to curves of critical values of physical parameters, thereby providing compact, closed-form criteria for the onset of instability. The derived analytical criteria are validated for swirling flows modelled by a cylindrical, differentially rotating annulus with axial flow induced by either a sliding inner cylinder, an axial pressure gradient or a radial temperature gradient combined with vertical gravity. These criteria unify and extend, to viscous and thermodiffusive differentially heated swirling flows, the Rayleigh criterion for centrifugally driven instabilities, the Ludwieg–Eckhoff–Leibovich–Stewartson criterion for isothermal swirling flows and the Goldreich–Schubert–Fricke criterion for non-isothermal azimuthal flows. Additionally, they predict oscillatory modes in swirling, differentially heated, visco-diffusive flows, thereby generalising the McIntyre instability criterion to these systems.

1. Introduction

The stability of swirling flows and their transition to turbulence represent a significant scientific challenge due to their occurrence in diverse natural and industrial settings.

These flows appear in natural phenomena such as tropical cyclones and tornadoes (Emanuel 1984, 2018), rotating convection (Lappa 2012; Oruba et al. 2017; Horn & Aurnou 2018; Oruba et al. 2018; Castano, Navarro & Herrero 2021; Soward et al. 2022), swirling vortex rings (Lifschitz & Hameiri 1993; Hattori & Hijiya 2010), meandering rivers (Park & Ahn 2019) and astrophysical flows, including magnetic tornadoes in the solar atmosphere (Knobloch & Spruit 1982; Busse & Pesch 2006; Wedemeyer-Böhm et al. 2012; Lopez et al. 2013; Tziotziou et al. 2023).

In engineering applications, swirling flows develop behind aircraft wingtips as trailing vortices (Eckhoff & Storesletten 1978; Leibovich & Stewartson 1983; Eckhoff 1984; Lucca-Negro & O’Doherty 2001; Di Pierro & Abid 2010; Billant & Gallaire 2013; Feys & Maslowe 2016) and play a crucial role in industrial processes such as combustion (Candel et al. 2014), isotope separation via centrifugation (Lamarsh 2017), cooling and lubrication of rotating machinery (Kreith 1969; Lee & Minkowycz 1989; Fénot et al. 2011; Seibold et al. 2022), wastewater purification (Ollis et al. 1991), everyday piping systems and physiological flows (Ault et al. 2016; Baratchi et al. 2020), oil drilling operations (MacAndrew et al. 1993), metal solidification (Vivès 1988) and crystal growth (Dhanaraj et al. 2010).

Hydrodynamic modelling describes swirling flows as a result of the combined effects of rotation and shear in two orthogonal directions (e.g. azimuthal and axial in cylindrical geometry). Swirling flows can occur either in confined configurations common in engineering applications (e.g. heat transfer systems, motor turbines, etc.) or in open geometries typical of natural phenomena such as atmospheric flows.

The circular Couette flow (CCF), in which two coaxial cylinders can rotate independently, is the simplest system for experimental and theoretical studies of swirling flows (Hollerbach et al. 2023; Lueptow et al. 2023). Axial flow in this system can be induced by various mechanisms, including the sliding motion of the inner cylinder, as in sliding Couette flow (Deguchi & Nagata 2011); an external pressure gradient, as in annular Poiseuille flow (Heaton 2008b ; Masuda et al. 2008); or a radial temperature gradient combined with axial gravity, as in baroclinic convection (Bahloul et al. 2000; Lepiller et al. 2007; Wang & Chen 2022). These mechanisms yield the following different variants of CCF with axial flow that replicate swirling flows.

1. (i) Spiral Couette flow (SCF): where axial flow is driven by inner-cylinder sliding (Ludwieg 1960, 1964; Ng & Turner 1982; Ali & Weidman 1993; Meseguer & Marques 2000).

2. (ii) Spiral Poiseuille flow (SPF): where axial flow is induced by an external pressure gradient (Takeuchi & Jankowski 1981; Wereley & Lueptow 1999; Meseguer & Marques 2002, Cotrell & Pearlstein 2004; Cotrell & McFadden 2005; Meseguer & Marques 2005; Brockmann et al. 2023).

3. (iii) Baroclinic Couette flow (BCF): where a radial temperature gradient and axial gravity cause baroclinic convection (Snyder & Karlsson 1964; Ali & Weidman 1990; Lepiller et al. 2008; Yoshikawa et al. 2013; Guillerm et al. 2015; Kang et al. 2015, 2023).

Note that in short Couette–Taylor cells and rotating cylindrical tanks, Ekman pumping induced by stationary or rotating endwalls can generate a secondary axial flow in the sidewall Stewartson boundary layer, which in turn gives rise to a spiralling motion within this region (Hart & Kittelman 1996; Lopez & Marques 2010); see also Gelfgat et al. (1996) and Vishnu & Sameen (2020). In the present study, however, we focus on Couette–Taylor cells with an infinite aspect ratio.

The stability of CCF has been extensively studied experimentally, theoretically and numerically (Hollerbach et al. 2023; Lueptow et al. 2023), as it provides a well-controlled system for investigating the onset of turbulence. Rayleigh’s inviscid stability criterion states that a flow with curved streamlines becomes unstable if angular momentum decreases outward. However, this criterion must be supplemented by linear stability analysis to account for viscosity and determine the instability threshold as a function of the radii ratio and the relative rotation velocities of the cylinders (Chandrasekhar 1961). Additionally, axial flows in a stationary cylindrical annulus are susceptible to wall-driven shear instabilities, with thresholds that depend on the radii ratio (Bahloul et al. 2000; Lepiller et al. 2007; Masuda et al. 2008; Heaton 2008b ; Deguchi & Nagata 2011; Wang & Chen 2022).

The presence of a radial temperature gradient and thermal diffusion introduces new destabilisation mechanisms in differentially rotating azimuthal flows, such as CCF. These include the visco-thermodiffusive extension of Rayleigh’s centrifugal instability, known as the Goldreich–Schubert–Fricke (GSF) instability (Acheson & Gibbons 1978; Maeder et al. 2013; Kirillov & Mutabazi 2017; Dymott et al. 2023), and the visco-thermodiffusive oscillatory McIntyre instability (McIntyre 1970; Labarbe & Kirillov 2021). Both instabilities can destabilise Rayleigh-stable flows, including quasi-Keplerian ones, depending on the Prandtl number and the direction of heating (Kirillov & Mutabazi 2017; Meyer et al. 2021).

A radial temperature gradient naturally arises in many applications of swirling flows. For example, recent simulations (Oruba et al. 2017) of rotating convection between two differentially heated horizontal plates in a shallow cylindrical domain revealed the formation of a vertical swirling base flow with a naturally developed hot core and an outward radial temperature gradient, consistent with observations of tropical cyclones (Emanuel 2018). In such non-isothermal visco-thermodiffusive swirling flows, all the previously mentioned linear instabilities may be present, along with algebraically growing and spatio-temporal instabilities (Eckhoff & Storesletten 1978; Gallaire & Chomaz 2003a , b ; Heaton & Peake 2006; Heaton 2008a ; Martinand et al. 2009) and nonlinear effects (Szeri & Holmes 1988), potentially leading to complex new interactions.

For isothermal, inviscid, incompressible swirling flows, the following three general analytical criteria for stationary centrifugal instability have been developed in the literature.

1. (i) The Ludwieg criterion for SCF in the narrow-gap limit (Ludwieg 1960, 1964).

2. (ii) The Leibovich-Stewartson criterion for swirling jets (Leibovich & Stewartson 1983).

3. (iii) The Eckhoff criterion, derived using a local geometrical optics approach (Eckhoff 1984).

The equivalence between the local and global approaches leading to these instability conditions – commonly referred to as the Ludwieg–Eckhoff–Leibovich–Stewartson (LELS) criterion – was formally established by Leblanc & Le Duc (2005). Emanuel (1984) rederived the LELS criterion by identifying a plane that locally contains all the velocity shear relative to a rotating coordinate system and performing computations in the resulting helical coordinate system. In a later experimental study, Hart & Kittelman (1996) used physical arguments to justify a result similar to that of Emanuel (1984), in the context of centrifugal instability of the helical secondary flow within the Stewartson sidewall boundary layer in a rotating cylinder with a more rapidly corotating upper lid.

However, the LELS criterion applies only to isothermal, non-stratified and inviscid swirling flows. In particular, by this reason it does not capture, e.g. the folded neutral stability surface observed in SCF and SPF through numerical linear stability analysis by Meseguer & Marques (2000, 2002, 2005), nor does it account for non-isothermal (Cotrell & McFadden 2005; Yoshikawa et al. 2013) or stratification (Yavneh et al. 2001; Di Pierro & Abid 2010) effects. Thus, similar to the case of rotating plane shear flows (Oxley & Kerswell 2024), there is a strong need for a unified theory of instabilities in swirling flows that incorporates axial shear, differential rotation, thermal stratification, viscosity and thermal diffusion and aligns with numerical and experimental results.

The present work is a comprehensive theoretical investigation of both isothermal and non-isothermal viscous and thermodiffusive swirling flows; it presents general criteria of instability in analytical form. It provides the complete development of our recent brief communication (Kirillov & Mutabazi 2024), which focused solely on centrifugal instability. More importantly, in the present work we extend this analytical approach to identify the conditions for oscillatory visco-thermodiffusive McIntyre instability in non-isothermal swirling flows.

The methodology of this study is based on the local geometrical optics approach, initially introduced in hydrodynamics for inviscid flows (Eckhoff & Storesletten 1978; Friedlander & Vishik 1991; Lifschitz & Hameiri 1991, 1993; Mathur et al. 2014; Ionescu-Kruse 2017) and later extended to visco-diffusive flows (Kirillov & Stefani 2013; Kirillov et al. 2014; Kirillov 2017, 2021, 2025; Kirillov & Mutabazi 2017, 2024; Labarbe & Kirillov 2021) and to multiple-diffusive flows (Singh & Mathur 2019, 2021; Vidal et al. 2019). We apply this approach to helical stationary solutions of the Navier–Stokes equations coupled with the energy equation in the Boussinesq–Oberbeck approximation, which – depending on boundary conditions – represent SCF, SPF or BCF.

We establish that the neutral stability curves in the plane spanned by the azimuthal and axial Reynolds numbers – corresponding to both stationary (for isothermal and non-isothermal flows) and oscillatory (for non-isothermal flows only) instabilities – form families that possess envelopes, regardless of whether they are parameterised by the azimuthal or axial wavenumber. Since the axial wavenumber can take any real value, these envelopes define the boundaries of the unions of individual stationary and oscillatory instability domains for specific values of wavenumbers.

We find that the envelopes, and thus, the instability domains, undergo splitting during the transition from Rayleigh-unstable to Rayleigh-stable flows. Notably, this splitting occurs at different Rossby number values for stationary and oscillatory instabilities of non-isothermal flows. This discrepancy offers a predictive tool for determining whether stationary or oscillatory instability will dominate in a visco-thermodiffusive swirling flow.

For Rayleigh-stable flows, as the azimuthal Reynolds number approaches infinity, an asymptotic line to the envelope of stationary instability domains corresponds to the inviscid LELS criterion for isothermal flows and provides an analytical expression for the unified LELS-GSF criterion in non-isothermal flows. In the isothermal case, we derive a compact closed-form expression for the envelope, extending the inviscid LELS criterion to viscous swirling flows and broadening its applicability across a wide range of azimuthal Reynolds numbers, from infinity to small but finite values. For non-isothermal flows, we establish explicit equations that define the asymptotic lines of the instability domains, offering a complete characterisation of their boundaries.

The paper is organised as follows. In § 2 we present the dimensionless nonlinear equations governing swirling flows in the Boussinesq–Oberbeck approximation, identify three main helical base states (detailed in Appendix A), and derive the linearised equations about these base states. Section 3 develops the geometrical optics asymptotic solution to the linearised equations, leading to the amplitude transport equations for the localised wavepacket evolving along the helical streamlines and the eikonal equation for its wavevector. In § 4 we derive the dispersion relation and demonstrate its connection to previous works in Appendix B. Section 5 introduces the new instability criteria, applying them to viscous isothermal swirling flows in § 5.1 and to visco-thermodiffusive non-isothermal swirling flows in § 5.2. Appendix C establishes the connection between our results and the classical inviscid LELS criterion. Appendix D relates our results to the criterion of Di Pierro & Abid (2010) for centrifugal instability in an inviscid swirling flow with radial density stratification but without mass diffusivity. Finally, § 6 presents the conclusions.

Figure 1. The helical base state as a superposition of the CCF $v_B(R)$ and an axial annular flow $w_B(R)$ in a differentially heated infinitely long cylindrical annulus.

2. Mathematical setting

2.1. Nonlinear governing equations

We consider an incompressible Newtonian fluid with constant reference density $\rho$ as well as constant thermal expansion coefficient $\alpha$ , kinematic viscosity $\nu$ and thermal diffusivity $\kappa$ . This fluid is confined within an infinitely long cylindrical annulus with a gap width $d = R_2 - R_1$ , where $R_1$ is the radius of the inner cylinder at temperature $T_1$ , rotating with angular velocity $\varOmega _1$ , and $R_2$ is the radius of the outer cylinder at temperature $T_2 = T_1 - \Delta T$ , rotating with angular velocity $\varOmega _2$ ; see figure 1. We denote the ratios of radii and angular velocities as

(2.1) \begin{equation} \eta =\frac {R_1}{R_2}, \quad \mu =\frac {\varOmega _2}{\varOmega _1}. \end{equation}

The system is subjected to a uniform gravity field with acceleration $g$ along the $Z$ axis of the cylindrical coordinates $(R, \varphi , Z)$ , which aligns with the common rotation axis of the cylinders (figure 1).

We choose $V_0 = \varOmega _1 R_1$ as the velocity scale, $d$ as the length scale, ${d}/{V_0}$ as the time scale and $\rho V_0^{2}$ as the pressure scale. Under the Boussinesq–Oberbeck approximation, in which fluid properties are assumed constant except for density, which varies linearly with temperature in the driving forces, the dimensionless governing equations can be written as

(2.2a) \begin{align} {\boldsymbol{\nabla}} {\boldsymbol{\cdot}} {\boldsymbol{u}}&=0, \end{align}

(2.2b) \begin{align} \frac {{\textrm{d}} {\boldsymbol{u}}}{{\textrm{d}} t}+{\boldsymbol{\nabla}} p-\frac {1}{\textit{Re}}{\nabla }^{2}{\boldsymbol{u}}&+ \left (\gamma \frac {v^{2}}{r}\boldsymbol{e}_r-Ri\boldsymbol{e}_z\right )\theta =0, \end{align}

(2.2c) \begin{align} \frac {{\textrm{d}}\theta }{{\textrm{d}} t}&-\frac {1}{\textit{Re}\textit{Pr}}{\boldsymbol{\nabla}}^{2}\theta =0, \end{align}

where $ \textrm{d}/{{\textrm{d}}t} = {\partial }/{\partial t} + {\boldsymbol{u}} {\boldsymbol{\cdot}} {\boldsymbol{\nabla}}$ , $ r = {R}/{d}$ and $ z = {Z}/{d}$ . Here, $ p$ is the pressure, $ {\boldsymbol{u}} = (u, v, w)$ is the velocity field and $\theta = (T-T_2)/{\Delta T}$ represents the temperature deviation from the reference temperature $ T_2$ . The parameter $ \gamma = \alpha \Delta T \gt 0$ corresponds to outward heating ( $ T_1 \gt T_2$ ), while $ \gamma \lt 0$ corresponds to inward heating ( $ T_1 \lt T_2$ ). The velocity field satisfies the no-slip boundary conditions at the surfaces of the cylinders.

The last two terms on the left-hand side of (2.2b ) represent the centrifugal buoyancy ( $ \gamma ({v^{2}}/{r}) \boldsymbol{e}_r$ ) and Archimedean buoyancy ( $ -Ri \boldsymbol{e}_z$ ), respectively.

The dimensionless control parameters in (2.2) are defined as

(2.3) \begin{eqnarray} \quad \textit{Re}=\frac {V_0 d}{\nu }, \quad \textit{Pr}=\frac {\nu }{\kappa }, \quad S=\frac {V_0}{W_0}, \quad Ri=\frac {W_T}{W_0}\frac {1}{S \textit{Re}}, \end{eqnarray}

where $ \textit{Re}$ is the Reynolds number associated with the rotation of the inner cylinder, $\textit{Pr}$ is the Prandtl number, $W_0$ is the characteristic axial velocity of the flow, $ S$ is the swirl parameter (Ali & Weidman 1990, 1993) and $ Ri$ is the Richardson number. The Richardson number is defined using the characteristic thermal velocity (Choi & Korpela 1980):

(2.4) \begin{equation} W_T=\frac {\gamma g d^{2}}{\nu }. \end{equation}

Using $ W_T$ , we can also introduce the Grashof number ( $ Gr$ ), which characterises the strength of the baroclinic flow:

(2.5) \begin{equation} Gr=\frac {W_T d}{\nu }. \end{equation}

It is worth noting that if $ W_0 = W_T$ then $ S = {\textit{Re}}/{Gr}$ , and consequently, the Richardson and Grashof numbers are related as

(2.6) \begin{equation} Ri = \frac {1}{S\textit{Re}} = \frac {Gr}{\textit{Re}^{2}}. \end{equation}

2.2. Helical base state

In an infinitely long cylindrical annulus, the base flow is a steady helical flow, invariant in both the azimuthal and axial directions. The temperature and velocity field depend only on the radial coordinate $r$ , while the pressure varies with both $r$ and $z$ (Kirillov & Mutabazi 2024):

(2.7a,b,c) \begin{align} {\boldsymbol{u}}_B(r)=(0,V(r),S^{-1}W(r)), \quad \theta _B(r)=\varTheta (r), \quad p_B(r,z)=p_1(r)+z p_2. \end{align}

Explicit expressions for the radial distributions of temperature $\varTheta (r)$ , azimuthal velocity $V(r)$ , axial velocity $W(r)$ and the characteristic axial velocity $W_0$ for various particular helical base flows – including isothermal and non-isothermal cases, such as SCF, SPF and BCF – are derived in Appendix A. Table 1 presents the base state profiles and their values at the dimensionless geometric mean radius $r_g$  (A4).

Table 1. Expressions of base flow in cylindrical annulus and their values evaluated at the geometric mean radius $r_{g} = \sqrt {r_{1} r_{2}}$ , where $r_1= {\eta }/(1-\eta)$ and $r_2={1}/(1-\eta)$ , for CCF, BCF, SCF and SPF.

2.3. Linearisation

To test the stability of the base state (2.7), we introduce small three-dimensional perturbations $ ({\boldsymbol{u}}^{\prime}, {p^{\prime}},\theta ^{\prime} )$ and linearise the nonlinear equations (2.2) about the base state and write the linearised equations in the matrix form (Kirillov 2017, 2021)

(2.8) \begin{equation} {\left (\! \begin{array}{cc} \partial _t+ {\mathcal{U}}+\frac {2V}{r}\gamma \varTheta \boldsymbol{ e}_r\boldsymbol{ e}_{\varphi }^{T}+{{\boldsymbol{u}}}_B{\boldsymbol{\cdot}} {\boldsymbol{\nabla}} -\frac {1}{\textit{Re}}{\boldsymbol{\nabla}}^2 &\! \gamma \frac {V^2}{r}\boldsymbol{ e}_r-Ri\boldsymbol{ e}_z \\ ({\boldsymbol{\nabla}} \varTheta )^T &\! \partial _t+{{\boldsymbol{u}}}_B{\boldsymbol{\cdot}} {\boldsymbol{\nabla}} -\frac {1}{\textit{Re}\textit{Pr}}{\boldsymbol{\nabla}}^2\\ \end{array}\! \right )}{\left (\! \begin{array}{c} {{\boldsymbol{u}}}^{\prime} \\ \theta ^{\prime} \\ \end{array}\! \right )}=-{\left (\! \begin{array}{c} {\boldsymbol{\nabla}} p^{\prime} \\0 \\ \end{array}\! \right )}, \end{equation}

where ${\boldsymbol{u}}^{\prime}$ is subject to the incompressibility constraint

(2.9) \begin{equation} {\boldsymbol{\nabla}} {\boldsymbol{\cdot}} {\boldsymbol{u}}^{\prime}=0. \end{equation}

The gradients of the base state are given by (Kirillov & Mutabazi 2024)

(2.10) \begin{equation} {\mathcal{U}}={\boldsymbol{\nabla}} {{\boldsymbol{u}}}_B=\left ( \begin{array}{ccc} 0 & -\varOmega & 0 \\ (1+2Ro)\varOmega & 0 & 0 \\ \frac {1}{S} DW & 0 & 0 \\ \end{array} \right )\!, \quad {\boldsymbol{\nabla}} \varTheta = \left (\begin{array}{c} D\varTheta \\ 0\\ 0\\ \end{array}\right )\!,\end{equation}

with $D=\textrm{d}/{{\textrm{d}}r}$ , $\varOmega ={V}/{r}$ and the Rossby number defined as (Kirillov & Stefani 2013)

(2.11) \begin{equation} Ro=\frac {rD\varOmega }{2\varOmega }. \end{equation}

The perturbations must vanish at the boundaries of the cylindrical annulus; that is, ${\boldsymbol{u}}^{\prime} = 0$ and $\theta ^{\prime} = 0$ at both $r = r_1$ and $r = r_2$ . However, these boundary conditions are not applied in the analysis that follows, as we restrict our attention to local instabilities.

3. Geometrical optics equations

Using a small parameter $0\lt \epsilon \ll 1$ , we represent the perturbations as asymptotic expansions (Eckhoff & Storesletten 1978; Friedlander & Vishik 1991; Lifschitz & Hameiri 1991):

(3.1) \begin{align} {\boldsymbol{u}}^{\prime } &=\big({\boldsymbol{u}}^{(0)}(\boldsymbol{x}, t){+}\epsilon {\boldsymbol{u}}^{(1)}(\boldsymbol{x}, t)\big){e}^{\frac {{\textrm{i}} \varPhi (x, t)}{ \epsilon }}+\epsilon {\boldsymbol{u}}^{(r)}(\boldsymbol{x}, t,\epsilon ){\, +\, \textrm{c.c}.},\nonumber \\ \theta ^{\prime } &=\big(\theta ^{(0)}(\boldsymbol{x}, t){+}\epsilon \theta ^{(1)}(\boldsymbol{x}, t)\big){e}^{\frac {{\textrm{i}} \varPhi (x, t)}{ \epsilon }}+\epsilon \theta ^{(r)}(\boldsymbol{x}, t,\epsilon ){\, +\,\textrm{c.c}.},\nonumber \\ p^{\prime } &=\big(p^{(0)}(\boldsymbol{x}, t){+}\epsilon p^{(1)}(\boldsymbol{x}, t)\big){e}^{\frac {{\textrm{i}} \varPhi (x, t)}{ \epsilon }}+\epsilon p^{(r)}(\boldsymbol{x}, t,\epsilon ){\, +\, \textrm{c.c}.} \end{align}

Here, ${\textrm{i}} = \sqrt {-1}$ and $\varPhi$ is generally a complex-valued scalar function representing the phase of the wave or the eikonal. The quantities ${\boldsymbol{u}}^{(i)}$ , $\theta ^{(i)}$ and $p^{(i)}$ , for $i = 0, 1, 2, \ldots$ , are complex amplitudes, and the complex remainder terms ${\boldsymbol{u}}^{(r)}$ , $\theta ^{(r)}$ and $p^{(r)}$ are assumed to be uniformly bounded in $\epsilon$ over any fixed time interval. The notation c.c. denotes complex conjugate terms, included to ensure that the expansions are real-valued, consistent with the real-valued nature of the primed quantities on the left-hand side.

Maslov (1986) observed that high-frequency oscillations $\exp ({\textrm{i}}\epsilon ^{-1}\varPhi (\boldsymbol{ x}, t) )$ quickly die out because of viscosity unless one assumes a quadratic dependency of viscosity and diffusivity on the small parameter $\epsilon$ : $\nu =\epsilon ^2 \widehat {\nu }$ and $\kappa =\epsilon ^2\widehat {\kappa }$ . Hence, following Maslov (1986), we have $\textit{Re} = \epsilon ^{-2}\widehat {\textit{Re}}$ ; see also Lagnado et al. (1984), Craik & Criminale (1986), Leblanc (2003).

With the assumptions made, substitution of the asymptotic series into the incompressibility condition (2.9) and collection of terms of the order $\epsilon ^{-1}$ and $\epsilon ^{0}$ yield

(3.2) \begin{align} \epsilon ^{-1}:\quad {{\boldsymbol{u}}}^{(0)} {\boldsymbol{\cdot}} {\boldsymbol{\nabla}} \varPhi &= 0, \end{align}

(3.3) \begin{align} \epsilon ^{0}:\quad {\boldsymbol{\nabla}} {\boldsymbol{\cdot}} {{\boldsymbol{u}}}^{(0)} + {\textrm{i}} {{\boldsymbol{u}}}^{(1)} {\boldsymbol{\cdot}} {\boldsymbol{\nabla}} \varPhi &= 0. \end{align}

A similar procedure applied to the linearised Navier–Stokes and energy equations (2.8) yields the two systems of equations

(3.4a) \begin{align} \epsilon ^{-1}:\quad &{\left (\begin{array}{cc} \dfrac {\partial \varPhi }{\partial t}+{{\boldsymbol{u}}}_B{\boldsymbol{\cdot}} {\boldsymbol{\nabla}} \varPhi & 0\\ 0 & \dfrac {\partial \varPhi }{\partial t}+{{\boldsymbol{u}}}_B{\boldsymbol{\cdot}} {\boldsymbol{\nabla}} \varPhi \\ \end{array}\right )}{\left (\begin{array}{c} {{\boldsymbol{u}}}^{(0)}\\ \theta ^{(0)} \end{array}\right )}=- p^{(0)}{\left (\begin{array}{c} {\boldsymbol{\nabla}} \varPhi \\ 0\\ \end{array} \right )}, \end{align}

(3.4b) \begin{align} \epsilon ^{0}:\quad &{\textrm{i}} {\left (\begin{array}{cc} \dfrac {\partial \varPhi }{\partial t}+{{\boldsymbol{u}}}_B{\boldsymbol{\cdot}} {\boldsymbol{\nabla}} \varPhi & 0\\ 0 & \dfrac {\partial \varPhi }{\partial t}+{{\boldsymbol{u}}}_B{\boldsymbol{\cdot}} {\boldsymbol{\nabla}} \varPhi \\ \end{array} \right )}{\left (\begin{array}{c} {{\boldsymbol{u}}}^{(1)}\\ \theta ^{(1)} \end{array}\right )}=-{\textrm{i}} p^{(1)} {\left (\begin{array}{c} {\boldsymbol{\nabla}} \varPhi \\ 0\\ \end{array} \right )}\nonumber \\ &-{\left (\begin{array}{cc} \dfrac {\textrm{d}}{{\textrm{d}} t}+{\mathcal{U}} +\dfrac {1}{\widehat {\textit{Re}}}\left ({\boldsymbol{\nabla}} \varPhi \right )^2 & 0\\ 0 & \dfrac {\textrm{d}}{{\textrm{d}} t} +\dfrac {1}{\vphantom{\dfrac{1}{2}}\widehat {\textit{Re}}\textit{Pr}}\left ({\boldsymbol{\nabla}} \varPhi \right )^2\\ \end{array}\right )}\left (\begin{array}{c} {{\boldsymbol{u}}}^{(0)}\\ \theta ^{(0)} \end{array}\right )\nonumber \\ &-{\left ( \begin{array}{cc} \dfrac {2V}{r}\gamma \varTheta \boldsymbol{ e}_r{\textbf{e}}_{\varphi }^{T} & \gamma \dfrac {V^{2}}{r}\boldsymbol{ e}_r-{Ri}\boldsymbol{ e}_z\\ \left ({\boldsymbol{\nabla}}\varTheta \right )^{T} & 0\\ \end{array}\right )}{\left (\begin{array}{c} {{\boldsymbol{u}}}^{(0)}\\ \theta ^{(0)} \end{array}\right )}- {\left (\begin{array}{c} {\boldsymbol{\nabla}} p^{(0)}\\ 0\\ \end{array} \right )}, \end{align}

where $ \textrm{d}/{{\textrm{d}}t} = {\partial }/{\partial t} + {\boldsymbol{u}}_B {\boldsymbol{\cdot}} {\boldsymbol{\nabla}}$ .

Taking the dot product of the first of the equations in the system (3.4a ) with ${\boldsymbol{\nabla}} \varPhi$ under the constraint (3.2), we find that, for ${\boldsymbol{\nabla}} \varPhi \ne 0$ ,

(3.5) \begin{equation} p^{(0)}=0. \end{equation}

Under the condition (3.5) the system (3.4a ) has a non-trivial solution if the determinant of the $4\times 4$ matrix on its left-hand side is vanishing. This gives us a four-fold characteristic root corresponding to the Hamilton–Jacobi equation

(3.6) \begin{equation} \frac {\partial \varPhi }{\partial t}+{{\boldsymbol{u}}}_B{\boldsymbol{\cdot}} {\boldsymbol{\nabla}} \varPhi =0, \end{equation}

with the initial data $\varPhi (\boldsymbol{ x},0)=\varPhi _0(\boldsymbol{ x})$ . Taking the gradient of (3.6) yields the eikonal equation

(3.7) \begin{equation} \frac {\textrm{d}}{{\textrm{d}}t}{\boldsymbol{\nabla}} \varPhi =-{\boldsymbol{\nabla}} {{\boldsymbol{u}}}_B {\boldsymbol{\cdot}} {\boldsymbol{\nabla}} \varPhi =-{\mathcal{U}}^{T} {\boldsymbol{\nabla}} \varPhi \end{equation}

with the initial condition ${\boldsymbol{\nabla}}\varPhi (\boldsymbol{ x},0)={\boldsymbol{\nabla}}\varPhi _0(\boldsymbol{ x})$ , where ${\mathcal{U}}^{T}$ denotes the transposed $3 \times 3$ matrix ${\mathcal{U}}$ defined by (2.10).

Relations (3.6) and (3.7) allow us to reduce the system (3.4b ) to

(3.8a) \begin{align} \left (\frac {\textrm{d}}{{\textrm{d}}t}+{\mathcal{U}}+\frac {1}{\widehat {\textit{Re}}}\left ({\boldsymbol{\nabla}} \varPhi \right )^2 {\mathcal{I}}+\frac {2V}{r}\gamma \varTheta \boldsymbol{ e}_r\boldsymbol{ e}_{\varphi }^{T}\right ){{\boldsymbol{u}}}^{(0)}+\left (\gamma \frac {V^2}{r}\boldsymbol{ e}_r-Ri\boldsymbol{ e}_z\right )\theta ^{(0)}&=-{\textrm{i}}{\boldsymbol{\nabla}} \varPhi p^{(1)},\end{align}

(3.8b) \begin{align} \left ({\boldsymbol{\nabla}}\varTheta \right )^T{{\boldsymbol{u}}}^{(0)}+\left (\frac {\textrm{d}}{{\textrm{d}}t} +\frac {1}{\widehat {\textit{Re}}\textit{Pr}}\left ({\boldsymbol{\nabla}} \varPhi \right )^2\right )\theta ^{(0)}&=0,\end{align}

where ${\mathcal{I}}$ is the $3\times 3$ identity matrix, $\boldsymbol{ e}_r\boldsymbol{ e}_{\varphi }^{T}$ is a $3\times 3$ matrix and $ ({\boldsymbol{\nabla}}\varTheta )^T{{\boldsymbol{u}}}^{(0)}={\boldsymbol{\nabla}}\varTheta {\boldsymbol{\cdot}} {{\boldsymbol{u}}}^{(0)}$ .

Multiplying (3.8a ) with the vector ${\boldsymbol{\nabla}} \varPhi$ from the left and taking into account the relation (3.2), we isolate the pressure term

(3.9) \begin{equation} p^{(1)}={\textrm{i}}\frac {{\boldsymbol{\nabla}} \varPhi }{({\boldsymbol{\nabla}} \varPhi )^2}{\boldsymbol{\cdot}}\left [ \left (\frac {\textrm{d}}{{\textrm{d}}t}+{\mathcal{U}}+\frac {2V}{r}\gamma \varTheta \boldsymbol{ e}_r\boldsymbol{ e}_{\varphi }^{T}\right ){{\boldsymbol{u}}}^{(0)}+\left (\gamma \frac {V^2}{r}\boldsymbol{ e}_r-Ri\boldsymbol{ e}_z\right )\theta ^{(0)} \right ]\!. \end{equation}

Taking into account the identity

(3.10) \begin{equation} \frac {\textrm{d}}{{\textrm{d}}t}\big({\boldsymbol{\nabla}} \varPhi {\boldsymbol{\cdot}} {{\boldsymbol{u}}}^{(0)}\big)=\frac {{\textrm{d}} {\boldsymbol{\nabla}} \varPhi }{{\textrm{d}}t}{\boldsymbol{\cdot}} {{\boldsymbol{u}}}^{(0)}+{\boldsymbol{\nabla}} \varPhi {\boldsymbol{\cdot}} \frac {{\textrm{d}} {{\boldsymbol{u}}}^{(0)}}{{\textrm{d}}t}=0 \end{equation}

we modify (3.9) in the following way:

(3.11) \begin{align} p^{(1)}&={\textrm{i}}\frac {{\boldsymbol{\nabla}} \varPhi }{({\boldsymbol{\nabla}} \varPhi )^2}{\boldsymbol{\cdot}}\left [ \left ({\mathcal{U}}+\frac {2V}{r} \gamma \varTheta \boldsymbol{ e}_r\boldsymbol{ e}_{\varphi }^{T}\right ){{\boldsymbol{u}}}^{(0)}+\left (\gamma \frac {V^2}{r}\boldsymbol{ e}_r-Ri\boldsymbol{ e}_z\right )\theta ^{(0)} \right ]\nonumber \\[4pt] &\quad -{\textrm{i}} \frac {1}{({\boldsymbol{\nabla}} \varPhi )^2}\frac {{\textrm{d}} {\boldsymbol{\nabla}} \varPhi }{{\textrm{d}} t}{\boldsymbol{\cdot}} {{\boldsymbol{u}}}^{(0)}. \end{align}

Now using the eikonal equation (3.7) we transform the last term in (3.11) to obtain

(3.12) \begin{align} p^{(1)}&={\textrm{i}}\frac {2V}{r}\gamma \varTheta \frac {({\boldsymbol{\nabla}} \varPhi )^{T}\boldsymbol{ e}_r{\textbf{e}}_{\varphi }^{T}}{({\boldsymbol{\nabla}} \varPhi )^{2}}{{\boldsymbol{u}}}^{(0)}+{\textrm{i}}\left (\gamma \frac {V^2}{r}\frac {({\boldsymbol{\nabla}} \varPhi )^{T}\boldsymbol{ e}_r}{({\boldsymbol{\nabla}} \varPhi )^{2}}-Ri\frac {({\boldsymbol{\nabla}} \varPhi )^{T}\boldsymbol{ e}_z}{({\boldsymbol{\nabla}} \varPhi )^{2}}\right )\theta ^{(0)}\nonumber \\[4pt] &\quad+2{\textrm{i}} \frac {({\boldsymbol{\nabla}} \varPhi )^{T}{\mathcal{U}}}{({\boldsymbol{\nabla}} \varPhi )^{2}} {{\boldsymbol{u}}}^{(0)}. \end{align}

Substituting (3.12) into (3.8a ) we finally arrive at the transport equations for the leading-order amplitudes ${{\boldsymbol{u}}}^{(0)}$ and $\theta ^{(0)}$ of the localised wavepacket moving along the streamlines of the base flow written in the stationary frame (Kirillov & Mutabazi 2024):

(3.13) \begin{align} \frac {{\textrm{d}}{\boldsymbol u}^{(0)}}{{\textrm{d}}t}+\frac {|{\boldsymbol k}|^2}{{\widehat {\textit{Re}}}}{\boldsymbol u}^{(0)} &=-\left (\mathcal{I}-\frac {{\boldsymbol k}{\boldsymbol k}^T}{|{\boldsymbol k}|^2}\right )\left (\gamma \frac {V^2}{r}{\boldsymbol e}_r-Ri{\boldsymbol e}_z\right )\theta ^{(0)}\nonumber \\[4pt]&\quad -\left (\mathcal{I}-2 \frac {{\boldsymbol k}{\boldsymbol k}^T}{|{\boldsymbol k}|^2}\right ) {\mathcal{U}}{\boldsymbol u}^{(0)}-2\gamma \varTheta \varOmega \left ({\mathcal{I}}-\frac {{\boldsymbol k}{\boldsymbol k}^T}{|{\boldsymbol k}|^2}\right ){\boldsymbol e}_r{\boldsymbol e}_{\varphi }^{T}{\boldsymbol u}^{(0)},\nonumber \\ \frac {{\textrm{d}} \theta ^{(0)}}{{\textrm{d}}t}+\frac {|{\boldsymbol k}|^2}{{\widehat {\textit{Re}}}\textit{Pr}}\theta ^{(0)}&=-\left ({\boldsymbol{\nabla}}\varTheta \right )^T{\boldsymbol u}^{(0)}. \end{align}

Here ${\boldsymbol k}{\boldsymbol k}^T$ is a $3\times 3$ matrix and we denote the wavenumber of the perturbations as $\boldsymbol{ k}={\boldsymbol{\nabla}} \varPhi$ . In this notation the eikonal equation (3.6) determines evolution of the wavevector in the stationary frame

(3.14) \begin{equation} \frac {{\textrm{d}}\boldsymbol{ k}}{{\textrm{d}} t}=-{\mathcal{U}}^T\boldsymbol{ k} \end{equation}

under the constraint (following from (3.2))

(3.15) \begin{equation} \boldsymbol{ k}{\boldsymbol{\cdot}} {{\boldsymbol{u}}}^{(0)}=0. \end{equation}

Note that the equations (3.13) differ from the amplitude equations derived earlier in Kirillov & Mutabazi 2017 by the Archimedean buoyancy term $-Ri\boldsymbol{e}_z$ and the term $S^{-1}DW$ in the matrix ${\mathcal{U}}$ .

4. Dispersion relation

The derivative $ {\partial {\boldsymbol k}}/{\partial t}$ of the wavevector ${\boldsymbol k}=(k_r,k_{\varphi },k_z)$ in the frame of the wavepacket, rotating about the vertical axis with the angular velocity $\varOmega$ is related to the derivative of this vector in the stationary frame as (Eckhardt & Yao 1995)

(4.1) \begin{equation} \frac {{\textrm{d}}{\boldsymbol k}}{{\textrm{d}}t}=\frac {\partial {\boldsymbol k}}{\partial t}+{\mathcal{J}}{\boldsymbol k}, \quad {\mathcal{J}}=\left ( \begin{array}{rrr} 0 & -\varOmega & 0 \\ \varOmega & 0 & 0 \\ 0 & 0 & 0 \\ \end{array} \right )\!\!. \end{equation}

Taking into account (4.1) in (3.14), we get

(4.2) \begin{equation} \frac {\partial {\boldsymbol k}}{\partial t}=\left ( \begin{array}{ccc} 0 & -2\varOmega \textit{Ro} & -S^{-1}DW \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{array} \right )\left ( \begin{array}{c} k_r \\ k_{\varphi } \\ k_z \\ \end{array} \right )\!\!. \end{equation}

Under the condition (Eckhoff & Storesletten 1978; Leibovich & Stewartson 1983; Emanuel 1984)

(4.3a,b) \begin{align} k_\varphi =-\overline {DW} k_z, \quad {\textrm{with}} \quad \overline {DW}=\frac {DW}{2\varOmega RoS}, \end{align}

the components of the wavevector $\boldsymbol{k}$ are time-independent in the rotating frame: $k_r = \textrm{const.}$ , $k_{\varphi }=\textrm{const}.$ , $k_z= \textrm{const.}$ Note that, since the vorticity of the helical base flow (2.7a ) has two components, $\omega _B = (0, -S^{-1}DW, 2\varOmega (1 + \textit{Ro}))$ , we can conclude from (4.3) that the azimuthal wavenumber of the perturbations is linked to their axial wavenumber through the azimuthal vorticity component ( $-S^{-1}DW$ ) of the base flow. These perturbations exhibit helical symmetry (Emanuel 1984), remaining invariant along circular helices with a pitch of $2\pi r \overline {DW}$ , and correspond to the most unstable and exponentially growing three-dimensional perturbations, as demonstrated in the works of Eckhoff & Storesletten (1978), Leibovich & Stewartson (1983), Eckhoff (1984), Billant & Gallaire (2013).

Note that perturbations violating the constraint (4.3) were shown already by Eckhoff & Storesletten (1978) to grow algebraically and to be weaker than those satisfying the constraint. For this reason, the former have been a less popular subject of investigation in the context of swirling flows; however, see the works of Heaton & Peake (2006), Di Pierro & Abid (2010) for more details on this subject. These weaker instabilities also lie outside the scope of the present paper.

Under the condition (4.3) the amplitude equations (3.13) are autonomous in the rotating frame. Following the standard procedure (see, e.g. Eckhardt & Yao (1995) and Friedlander & Vishik 1995) we write explicitly the material derivative $\textrm{d}/{{\textrm{d}}t} = {\partial }/{\partial t} + {\boldsymbol{u}}_B {\boldsymbol{\cdot}} {\boldsymbol{\nabla}}$ on the left-hand side of equations (3.13) and explicitly compute their right-hand side, exploiting the relation

(4.4) \begin{equation} |\boldsymbol{k}|^2=k_r^2+k_z^2\big[1+\overline {DW}^2\big]. \end{equation}

Since ${\boldsymbol{u}}_B(r)=(0,V(r),S^{-1}W(r))$ , the left-hand side of the resulting linear equations with constant coefficients does not contain the derivative with respect to $r$ . This allows us to seek for their solution in the modal form ${\boldsymbol{u}}^{(0)}, \theta ^{(0)}\sim e^{s t+{\textrm{i}} m\varphi +{\textrm{i}}k_z z}$ , where $s=\sigma + {\textrm{i}} \omega$ , $\sigma ,\omega \in {\mathbb{R}}$ , is the complex growth rate, and $m=k_\varphi r$ and $k_z$ are the integer azimuthal and real axial wavenumbers. This yields the system of linear algebraic equations

(4.5a) \begin{align} \lambda u_r^{(0)} &= - \frac {DW}{S \textit{Ro} } \frac {k_r k_z}{|\boldsymbol{ k}|^2} u_r^{(0)} -\frac {|\boldsymbol{ k}|^2}{\widehat {\textit{Re}}} u_r^{(0)} + 2\varOmega (1{-}\gamma \varTheta )\left (1 - \frac {k_r^2}{|\boldsymbol{ k}|^2}\right )u_{\varphi }^{(0)}\nonumber \\ &\quad - \left ( \gamma r \varOmega ^2\left (1 - \frac {k_r^2}{|\boldsymbol{ k}|^2}\right )+Ri \frac {k_rk_z}{|\boldsymbol{ k}|^2}\right )\theta ^{(0)}, \end{align}

(4.5b) \begin{align} \lambda u_{\varphi }^{(0)} &= - 2 \varOmega \left (\textit{Ro} + \frac {k_r^2 + k_z^2}{|\boldsymbol{ k}|^2}\right )u_r^{(0)} -\frac {|\boldsymbol{ k}|^2}{\widehat {\textit{Re}}} u_{\varphi }^{(0)} + (1-\gamma \varTheta )\frac {DW}{S \textit{Ro} }\frac {k_rk_z}{|\boldsymbol{ k}|^2}u_{\varphi }^{(0)}\nonumber \\ &\quad - \left ( \gamma r\varOmega ^2-Ri \frac {k_z}{k_r}\right )\frac {1}{2\varOmega }\frac {DW}{S \textit{Ro} }\frac {k_r k_z}{|\boldsymbol{ k}|^2}\theta ^{(0)}, \end{align}

(4.5c) \begin{align} \lambda u_z^{(0)} &= -\frac {|\boldsymbol{ k}|^2}{\widehat {\textit{Re}}} u_z^{(0)} - \frac {DW}{S Ro} \left (\textit{Ro} + \frac {k_z^2}{|\boldsymbol{ k}|^2}\right ) u_r^{(0)} - 2 \varOmega (1-\gamma \varTheta )\frac {k_rk_z}{|\boldsymbol{ k}|^2}u_{\varphi }^{(0)} \nonumber \\ &\quad + \gamma r\varOmega ^2\frac {k_r k_z}{|\boldsymbol{ k}|^2}\theta ^{(0)}+ Ri\left (1 - \frac {k_z^2}{|\boldsymbol{ k}|^2}\right )\theta ^{(0)}, \end{align}

(4.5d) \begin{align} \lambda \theta ^{(0)} &=-D\varTheta u_r^{(0)} -\frac {|\boldsymbol{ k}|^2}{\textit{Pr} \widehat {\textit{Re}}}\theta ^{(0)}, \end{align}

where

(4.6) \begin{equation} \lambda =\sigma +{\textrm{i}}\left (\omega +m\varOmega +\frac {k_zW}{S} \right )\!. \end{equation}

The generalised eigenvalue $\lambda$ includes the Doppler frequency $\omega _D=m\varOmega + {k_z W}/{S}$ of the fluid particle moving along the helix around the cylindrical annulus at the radius $r\in [r_1,r_2]$ . One can see that multiplying (4.5a ) with $ {k_r}/{k_z}$ , subtracting the result from (4.5b ) multiplied with $\overline {DW}$ , and taking into account (3.15) yields (4.5c ), which is therefore decoupled from the first two. The remaining equations (4.5a ), (4.5b ) and (4.5d ) result in the matrix eigenvalue problem for the $3 \times 3$ matrix ${\mathcal{H}}$ with the eigenvalue parameter $\lambda$ given by (4.6), where (Kirillov & Mutabazi 2024)

(4.7) \begin{equation} {\mathcal{H}}=\left (\! \begin{array}{ccc} \!-\dfrac { DW}{ SRo}\dfrac {k_r k_z}{|\boldsymbol{ k}|^2}-\dfrac {|\boldsymbol{ k}|^2}{\vphantom{frac{1}{2}}\widehat {\textit{Re}}} &\! 2\varOmega (1-\gamma \varTheta ) \!\left (1-\dfrac {k_r^2}{|\boldsymbol{ k}|^2}\! \right ) & -\left (r \gamma \varOmega ^2 \!\left (1-\dfrac {k_r^2}{|\boldsymbol{ k}|^2}\! \right ) \!+\! Ri\dfrac {k_r k_z}{|\boldsymbol{ k}|^2}\! \right ) \\ \!-2\varOmega \!\left (Ro+\dfrac {k_r^2+k_z^2}{|\boldsymbol{ k}|^2}\! \right ) &\! (1-\gamma \varTheta )\dfrac { DW}{ S Ro}\dfrac {k_r k_z}{|\boldsymbol{ k}|^2} \!-\! \dfrac {|\boldsymbol{ k}|^2}{\vphantom{frac{1}{2}}\widehat {\textit{Re}}}& \!-\left (r \gamma \varOmega ^2 - Ri\dfrac {k_z}{k_r}\right )\dfrac {1 }{2\varOmega } \dfrac { DW}{S Ro}\dfrac {k_r k_z}{|\boldsymbol{ k}|^2} \\ \!-D\varTheta & 0 &\! \dfrac {|\boldsymbol{k}|^2}{\vphantom{frac{1}{2}}\widehat {\textit{Re}}} \dfrac {\textit{Pr}-1}{\textit{Pr}}-\dfrac {|\boldsymbol{k}|^2}{\vphantom{frac{1}{2}}\widehat {\textit{Re}}} \\ \end{array}\! \right )\!. \end{equation}

We find the dispersion relation by computing the characteristic polynomial of ${\mathcal{H}}$ , i.e.

(4.8) \begin{equation} p(\lambda )=-\det ( {\mathcal{H}}-\lambda {\mathcal{I}})=a_3\lambda ^3+a_2\lambda ^2+a_1\lambda +a_0, \end{equation}

with the real coefficients

(4.9a) \begin{align} a_3&=1, \end{align}

(4.9b) \begin{align} a_2&=\frac {k_r k_z}{|\boldsymbol{ k}|^2} \frac {DW}{SRo}\gamma \varTheta +\frac {2\textit{Pr}+1}{\textit{Pr}}\frac {|\boldsymbol{ k}|^2}{\widehat {\textit{Re}}}, \end{align}

(4.9c) \begin{align} a_1&= \varOmega ^2\left (1 - \frac {k_r^2}{|\boldsymbol{ k}|^2}\right )\left \{4\left [\textit{Ro} + \frac {k_r^2 + k_z^2}{|\boldsymbol{ k}|^2}\right ](1-\gamma \varTheta ) - \gamma r D\varTheta \right \}+\frac {\textit{Pr}+2}{\textit{Pr}}\frac {|\boldsymbol{ k}|^4}{\widehat {\textit{Re}}^2} \nonumber \\ &\quad - \frac {(DW)^2}{S^2Ro^2}\frac {k_r^2k_z^2}{|\boldsymbol{ k}|^4}(1-\gamma \varTheta ) - Ri \frac {k_r k_z}{|\boldsymbol{ k}|^2}D\varTheta +\frac {\textit{Pr}+1}{\textit{Pr}}\frac {DW}{SRo}\frac {k_r k_z}{|\boldsymbol{ k}|^2}\frac {|\boldsymbol{ k}|^2}{ \widehat {\textit{Re}}}\gamma \varTheta , \end{align}

(4.9d) \begin{align} a_0&= \frac {D\varTheta DW}{SRo} \frac {k_z^2}{|\boldsymbol{ k}|^2} Ri (1-\gamma \varTheta )-\frac {k_r k_z}{|\boldsymbol{ k}|^2}\frac {|\boldsymbol{ k}|^2}{\widehat {\textit{Re}}}RiD\varTheta \nonumber \\ &\quad +\frac {|\boldsymbol{ k}|^2\varOmega ^2}{\widehat {\textit{Re}} \textit{Pr}}\left (1 - \frac {k_r^2}{|\boldsymbol{ k}|^2}\right )\left \{4\left [\textit{Ro} + \frac {k_r^2 + k_z^2}{|\boldsymbol{ k}|^2}\right ](1-\gamma \varTheta ) - \gamma r D\varTheta \textit{Pr} \right \}\nonumber \\ &\quad +\frac {|\boldsymbol{ k}|^2}{ \widehat {\textit{Re}} \textit{Pr}}\left ( \frac {|\boldsymbol{ k}|^2}{ \widehat {\textit{Re}}}-\frac {k_rk_z}{|\boldsymbol{ k}|^2}\frac {DW}{SRo}(1-\gamma \varTheta ) \right )\left (\frac {|\boldsymbol{ k}|^2}{ \widehat {\textit{Re}}} + \frac {k_rk_z}{|\boldsymbol{ k}|^2}\frac {DW}{SRo}\right )\!. \end{align}

In view of $Ri=({W_T}/{W_0})( {1}/{S\widehat {\textit{Re}}})$ the coefficients (4.9) and the matrix (4.7) reduce to those derived by Kirillov & Mutabazi (2017) in the limit $S \rightarrow \infty$ , as shown in Appendix B.

5. Local linear stability analysis of swirling flows

In the following we derive general stability conditions of the base flow by applying the Liénard–Chipart stability criterion (Kirillov & Mutabazi 2017; Kirillov 2021) to the characteristic polynomial (4.8) with the coefficients (4.9). This criterion guarantees $\lambda$  (4.6) to have only a negative real part:

(5.1a,b,c) \begin{align} a_0 \gt 0, \quad a_2 \gt 0,\quad a_1a_2-a_0\gt 0. \end{align}

We first examine isothermal flows ( $\gamma =0$ , $Ri=0$ ), such as spiral Couette and spiral Poiseuille flows, where the coefficients (4.9) simplify enough to allow for a fully analytical treatment. Next, we explore the most general scenario, involving both radial heating and Archimedean buoyancy, which includes BCF. For the sake of simplicity of notation in § 5, we omit the hat over the Reynolds numbers.

5.1. Viscous isothermal swirling flows

5.1.1. Viscous extension of the LELS centrifugal instability criterion

Assuming that $\gamma =0$ and $Ri=0$ in (4.9) automatically satisfies the condition $a_2\gt 0$ . This leaves only two inequalities, (5.1a ) and (5.1c ), to determine stability. The condition $a_0=0$ following from the inequality (5.1a ) defines the neutral stability curve as

(5.2) \begin{eqnarray} 4\varOmega ^2\left (1 - \frac {k_r^2}{|\boldsymbol{ k}|^2}\right )\left [\textit{Ro} + \frac {k_r^2 + k_z^2}{|\boldsymbol{ k}|^2}\right ]-\frac {k_r^2k_z^2}{|\boldsymbol{ k}|^4}\frac {DW^2}{S^2Ro^2} +\frac {|\boldsymbol{ k}|^4}{ \textit{Re}^2}=0, \end{eqnarray}

where $|\boldsymbol{ k}|$ is given by (4.4). The swirl parameter

(5.3) \begin{equation} S=\frac {\textit{Re}}{\textit{Re}_z} \end{equation}

is now defined via the axial Reynolds number $\textit{Re}_z = {W_0 d}/{\nu }$ , based on the characteristic axial flow velocity $W_0$ for isothermal flows. The details of its definition for SCF and SPF base flows are given in Appendix A.

Taking into account relations (4.3) and (4.4), we can rewrite (5.2) as

(5.4) \begin{equation} \frac {4\varOmega ^2k_z^2\big(\overline {DW}^2\textit{Ro} + \textit{Ro} + 1\big)\textit{Re}^2 + \big(\overline {DW}^2k_z^2 + k_r^2 + k_z^2\big)^3}{4\varOmega ^2\textit{Re}^2\big(\overline {DW}^2k_z^2 + k_r^2 + k_z^2\big)}=0, \end{equation}

which, given the positivity of its denominator, leads to an equation

(5.5) \begin{equation} q(k_z) = 0, \end{equation}

where

(5.6) \begin{align} q(k_z) &= \big(\overline {DW}^2 + 1\big)^3k_z^6 + 3k_r^2\big(\overline {DW}^2 + 1\big)^2k_z^4 \nonumber \\ &\quad + \big[3\big(\overline {DW}^2 + 1\big)k_r^4 + 4\varOmega ^2 \textit{Re}^2(\overline {DW}^2\textit{Ro} + \textit{Ro} + 1\big)\big]k_z^2 + k_r^6 \end{align}

is a real polynomial of degree six in $ k_z$ . Using the explicit expression (4.3b ) for $ \overline {DW}$ and (5.3) for the swirl parameter $ S$ , we can interpret (5.5) as defining a family of neutral stability curves in the $ (\textit{Re}_z, \textit{Re})$ plane.

Masuda et al. (2008) observed, using global numerical linear stability analysis, that the neutral stability curves of three-dimensional perturbations in plane Poiseuille flow with streamwise system rotation share a common envelope in the plane defined by the axial Reynolds number and the rotation number. This envelope exhibits both vertical and horizontal asymptotes at large values of these parameters. They further hypothesised that a similar envelope is likely to exist for SPF between concentric cylinders with $\eta =0.5$ (Meseguer & Marques 2002). Since (5.5) describes the family of neutral stability curves corresponding to the onset of centrifugal instability in isothermal swirling flows – including both SCF and SPF – we aim to test this hypothesis and search for an envelope of the family defined by (5.5).

To determine the envelope, we can compute the discriminant of the polynomial (5.6), using a standard tool in any modern computer algebra system:

(5.7) \begin{align} {\textrm{Disc}}_{k_z}(q) &= \left [27\big(\overline {DW}^2 + 1\big)k_r^4 + 16\varOmega ^2 \textit{Re}^2\big(\overline {DW}^2 \textit{Ro} + \textit{Ro} + 1\big)\right ]^2 \nonumber \\ &\quad \times \big(\overline {DW}^2\textit{Ro} + \textit{Ro} + 1\big)^4(\overline {DW}^2 + 1)^9(-16384)\textit{Re}^8\varOmega ^8k_r^6. \end{align}

Since the envelope corresponds to a portion of the discriminant set of the polynomial, given by ${\textrm{Disc}}_{k_z}(q) = 0$ (Hartman & Wintner 1953; Bruce & Giblin 1992), we focus only on the first factor in (5.7). Introducing the Rayleigh discriminant

(5.8) \begin{equation} N_{\varOmega }^2 = 4(1 + \textit{Ro})\varOmega ^2, \end{equation}

we express $ Ro$ in terms of $ N_{\varOmega }^2$ and set the first factor in (5.7) to zero, obtaining a compact analytical expression that defines the envelope of the neutral stability curves (5.2) for isothermal viscous swirling flows (Kirillov & Mutabazi 2024):

(5.9) \begin{equation} E(\textit{Re}_z, \textit{Re}) = \frac {N_{\varOmega }^2}{\varOmega ^2} - \frac {4\overline {DW}^2}{1 + \overline {DW}^2} + \frac {27}{4\varOmega ^2} \frac {k_r^4}{\textit{Re}^2} = 0. \end{equation}

Alternatively, the envelope (5.9) can be derived using the classical approach (Hartman & Wintner 1953; Bruce & Giblin 1992) by first differentiating (5.5) with respect to $ k_z$ , resulting in a biquadratic equation in $k_z$ . Solving this equation explicitly and substituting the solution back into (5.5) leads directly to the envelope equation (5.9).

According to Hull (2020) (see also Bonnans & Shapiro 2000, Milgrom & Segal 2002), the determination of an envelope can be formulated as a parametric optimisation problem. To illustrate this third approach, we present an example of such a problem that leads to the envelope (5.9).

First, we note that (5.4) can be rewritten as

(5.10) \begin{equation} \textit{Re} = \frac {\big(k_r^2 +\big(1+\overline {DW}^2\big)k_z^2\big)^{3/2}}{2\varOmega k_z\big(-1-Ro\big(1+\overline {DW}^2\big)\big)^{1/2}}. \end{equation}

Since the right-hand side of (5.10) contains $ \textit{Re}$ through the swirl parameter (5.3), which itself is embedded in $ \overline {DW}$  (4.3b ), this equation implicitly defines a neutral stability curve in the $ (\textit{Re}_z, \textit{Re})$ plane. Alternatively, for a fixed $ S$ , we can treat (5.10) as defining $ \textit{Re}$ as a function of $ k_z$ and seek to minimise it with respect to the axial wavenumber.

At the minimiser, given by

(5.11) \begin{equation} {k_{z,c} = k_r \big[2 \big(1+ \overline {DW}^2\big)\big]^{-1/2}, } \end{equation}

the minimal Reynolds number is

(5.12) \begin{equation} {\textit{Re}_c = \frac {3\sqrt {3}k_r^2}{2\varOmega } \left ( \frac {4\overline {DW}^2}{1+\overline {DW}^2} - \frac {N_{\varOmega }^2}{\varOmega ^2} \right )^{-1/2}= \frac {3\sqrt {3}k_r^2}{4\varOmega } \left (\frac {\overline {DW}^2}{1+\overline {DW}^2} - (1+Ro) \right )^{-1/2}. } \end{equation}

Through straightforward algebraic manipulations, the condition $ \textit{Re} = \textit{Re}_c$ leads directly to the envelope equation (5.9), in which we, again, should take into account (5.3) and (4.3b ).

The azimuthal component, $-S^{-1}DW$ , of vorticity $\omega _B$ due to the axial flow destabilises the azimuthal flow in the annulus: it decreases the threshold of the Rayleigh-unstable flows ( $1+Ro\lt 0$ ) whereas the Rayleigh-stable flows are destabilised by the axial flow, if

(5.13) \begin{equation} 0\lt 1+\textit{Ro} \lt \frac {\overline {DW}^2}{1+\overline {DW}^2} \end{equation}

or

(5.14) \begin{equation} \frac {N_\varOmega ^2}{\varOmega ^2}-\frac {4{\overline {DW}}^2}{1+{\overline {DW}}^2} \lt 0. \end{equation}

In (5.14) we recover the inviscid LELS instability criterion of isothermal swirling flow at $\textit{Re}\rightarrow \infty$ ; see Appendix C.

In the absence of the axial flow ( $\overline {DW}=0$ ), the inequality $E(\textit{Re}_z, \textit{Re})\lt 0$ exactly reproduces the result of Eckhardt & Yao (1995) for the centrifugal instability of the viscous Couette–Taylor flow:

(5.15) \begin{equation} \frac {N_{\varOmega }^2}{\varOmega ^2} + \frac {27}{4\varOmega ^2} \frac {k_r^4}{\textit{Re}^2} \lt 0. \end{equation}

5.1.2. Growth rate of the centrifugal LELS instability of isothermal viscous swirling flows

The dispersion relation (4.8) factorises for isothermal viscous swirling flows as

(5.16) \begin{equation} \left (\lambda ^2 + 2\frac {|\boldsymbol{k}|^2}{\textit{Re}}\lambda + \frac {|\boldsymbol{k}|^4}{\textit{Re}^2} + 4\varOmega ^2\frac {k_z^2}{|\boldsymbol{k}|^2}\big(1+Ro\big(1+\overline {DW}^2\big)\big) \right )\left (\lambda +\frac {|\boldsymbol{k}|^2}{\textit{Re}\textit{Pr}}\right )\!, \end{equation}

where $|\boldsymbol{k}|^2$ is defined in (4.4). The simple root associated with the second factor is always negative and corresponds to damped modes of the perturbation. The other two roots,

(5.17) \begin{equation} \lambda =-\frac {|\boldsymbol{k}|^2}{\textit{Re}}\pm 2\varOmega \frac {k_z}{|\boldsymbol{k}|}\sqrt {-1-Ro\big(1+\overline {DW}^2\big)}, \end{equation}

determined by the quadratic polynomial in the first factor, can exhibit positive real parts if and only if $\textit{Re}$ exceeds the marginal value specified on the right-hand side of (5.10), which was derived from the condition $a_0 = 0$ . In this case, one mode is damped and the other is amplified with a frequency $m\varOmega +({k_zW}/{S})$ , according to (4.6).

This analysis confirms that $a_0 \lt 0$ is the sole instability condition following from the three inequalities of the criterion (5.1) for viscous isothermal swirling flows. The envelope (5.9) of the individual instability domains specified by the inequality $a_0 \lt 0$ yields a viscous extension of the inviscid LELS instability criterion (5.14) for isothermal swirling flows.

The growth rate of the centrifugal instability can be estimated analytically at any position within the annulus using (5.17), particularly at the geometric mean radius, for which the expressions for $DW_g$ are provided in table 1.

Specifically, for CCF ( $W = 0$ ), we obtain

(5.18) \begin{equation} \lambda = -\frac {|\boldsymbol{k}|^2}{\textit{Re}} \pm 2\varOmega \frac {k_z}{|\boldsymbol{k}|}\sqrt {-1 - Ro}. \end{equation}

Hence, the marginal stability condition for Rayleigh-unstable CCF $(1 + \textit{Ro} \lt 0)$ is given by

(5.19) \begin{equation} \textit{Re}_m^{CCF} = \pm \frac {|\boldsymbol{k}|^3}{k_z}\frac {1}{2\varOmega \sqrt {-1 - Ro}}. \end{equation}

For Rayleigh-stable CCF $(1 + \textit{Ro} \gt 0)$ , the epicyclic frequency is given by $N_{\varOmega } = 2\varOmega \sqrt {1 + \textit{Ro}}$ .

5.1.3. Spiral Couette flow

To illustrate the inviscid LELS criterion (5.14) and its viscous extension, $E(\textit{Re}_z, \textit{Re})\lt 0$ , with $E$ given by (5.9), we consider the enclosed SCF as a base flow. Then, $W_0=W_1$ , i.e. the sliding speed of the inner cylinder, the dimensionless axial velocity $W(r)$ is given by (A17) and (A18), while $V(r)$ is defined by (A2); see table 1.

• Parameterisation by the axial wavenumber $k_z$ : the neutral stability curves (5.5) for the enclosed SCF with $W(r)$ given by (A17) and (A18), $\eta = 0.4$ and $k_r = 2\sqrt {2}$ are shown in the $(\textit{Re}_z,\textit{Re})$ plane in figure 2 for different values of $k_z$ and $\mu$ as black solid curves bounding the greenish centrifugal LELS instability domains. All computations for the SCF are performed at the geometric mean radius (A4).

For Rayleigh-unstable flows ( $N_{\varOmega }^2\lt 0$ or $\textit{Ro} \lt -1$ ), their envelope (5.9) is a single curve, shown as a red thick curve in figure 2(a) for $\mu =0$ . The envelope has a maximum $\textit{Re}=\textit{Re}_0$ when $\textit{Re}_z=0$ and a horizontal asymptote $\textit{Re}=\textit{Re}_{\infty }$ as $|\textit{Re}_z|\rightarrow \infty$ , where

(5.20a,b) \begin{align} \textit{Re}_0=\frac {3\sqrt {3}k_r^2}{4\varOmega \sqrt {-Ro-1}}, \quad \textit{Re}_\infty =\frac {3\sqrt 3 k_r^2}{4\varOmega \sqrt {-Ro}}. \end{align}

The asymptotes $\textit{Re}=\textit{Re}_{\infty }$ are shown in figure 2 by the black dot-dashed lines.

Figure 2. Isothermal enclosed SCF with $\eta = 0.4$ and $k_r = 2\sqrt {2}$ . (a,c) Greenish centrifugal LELS instability domains in the $(\textit{Re}_z, \textit{Re})$ plane, parameterised by $k_z$ , and their envelope (5.9) (thick red curves) for (a) $\mu = 0$ and $\textit{Re}_{\infty }\approx 18.78$ , and (c) $\mu = 0.5$ and $\textit{Re}_{\infty }\approx 17.71$ . (b) Variation of $k_z$ from $0$ to the terminal value $k_z^{\textit{term}} = (\sqrt{2}/2)k_r = 2$ according to (5.22) as $\textit{Re}$ increases from $\textit{Re}_{\infty }\approx 18.78$ (dot-dashed line) to $\textit{Re}_0\approx 29.70$ (dotted line) for $\mu =0$ . (d) Variation of $k_z$ from $0$ to the terminal value $k_z^{\textit{term}} = ({\sqrt {-2Ro}}/{2}) k_r \approx 1.217$ according to (5.22) as $\textit{Re}$ increases from $\textit{Re}_{\infty }\approx 17.71$ (dot-dashed line) to infinity for $\mu =0.5$ .

As $\textit{Ro} \rightarrow -1$ , $\textit{Re}_0 \rightarrow \infty$ , and for Rayleigh-stable flows $(-1 \lt \textit{Ro} \lt 0)$ , the envelope splits into two curves, each with vertical tangents at $\textit{Re}_z = \pm \textit{Re}_z^{\textit{min}}$ , where

(5.21) \begin{equation} \textit{Re}_z^{\textit{min}}= \frac {3}{2}\frac {\sqrt {3}k_r^2}{DW}\big(1 + \sqrt {Ro+1}\big) \end{equation}

is the minimal critical axial Reynolds number destabilising Rayleigh-stable isothermal azimuthal flows. The vertical solid lines in figure 2(c) show $\pm \textit{Re}_z^{\textit{min}}\approx \pm 27.43$ .

Additionally, in figure 2(c) the classical inviscid LELS criterion (5.14) is represented by the oblique black solid lines touching the upper parts of the envelope (5.9) as $\textit{Re} \rightarrow \infty$ . While (5.14) does not apply to Rayleigh-unstable flows, its viscous extension, $E(\textit{Re}_z, \textit{Re})\lt 0$ , defines the stability boundary in this case, as given by the envelope (5.9); see figure 2(a).

We can interpret these results as indicating that the axial flow induced by the sliding inner cylinder introduces a destabilising mechanism into the Couette–Taylor flow. The value $\textit{Re}_{\infty }$ corresponds to the minimum rotation rate below which the sliding of the inner cylinder is no longer able to destabilise the SCF.

Expressing $\overline {DW}$ from the envelope equation (5.9), substituting the result into (5.11) and subsequently utilising the explicit expressions provided in (5.8) and (5.20) for isothermal swirling flows, we derive the following analytical expression, which describes the dependence of the critical axial wavenumber $k_z$ of the perturbations on the critical Reynolds number $\textit{Re}$ along the envelope (5.9):

(5.22) \begin{equation} |k_z|(\textit{Re})=k_r\frac {\sqrt {2}}{2}\sqrt {\frac {1-\left (\frac {\textit{Re}_{\infty }}{\textit{Re}}\right )^2}{1-\left (\frac {\textit{Re}_{\infty }}{\textit{Re}_0}\right )^2}}. \end{equation}

In particular, it implies that $|k_z|\rightarrow 0$ as $\textit{Re} \rightarrow \textit{Re}_{\infty }$ and reaches its terminal value at $\textit{Re} = \textit{Re}_0$ ,

(5.23) \begin{equation} k_z^{\textit{term}}=\frac {\sqrt {2}}{2} k_r, \end{equation}

for Rayleigh-unstable flows (figure 2 b), and

(5.24) \begin{equation} k_z^{\textit{term}}=\frac {\sqrt {-2Ro}}{2} k_r, \end{equation}

as $\textit{Re} \rightarrow \infty$ for Rayleigh-stable flows (figure 2 d), agreeing with numerical studies by Ludwieg (1964), Ng & Turner (1982), Ali & Weidman (1993) and Meseguer & Marques (2000).

Additionally, figure 2(b,d) shows that the spectrum of axial wavenumbers for perturbations differs between Rayleigh-unstable and Rayleigh-stable flows. In Rayleigh-unstable flows (figure 2 b), the destabilising centrifugal force broadens the wavenumber range, whereas in Rayleigh-stable flows (figure 2 d), it is limited to lower $ k_z$ values. Long axial wavelength perturbations are favoured as $ \textit{Re}$ approaches $ \textit{Re}_{\infty }$ (and $\textit{Re}_z \rightarrow \infty$ ).

From (5.9), we determine that at the points of the envelope,

(5.25) \begin{equation} { \textit{Re} = \textit{Re}_{\infty }\left (1 + \frac {27k_r^4}{8DW^2\textit{Re}_z^2}+ \cdots \!\right ) } \end{equation}

as $ \textit{Re}_z \to \infty$ , where $ \textit{Re}_{\infty }$ is defined by (5.20b ).

Substituting (5.25) into (5.22), we find the relation between axial and radial wavenumbers on the stability boundary,

(5.26) \begin{equation} k_z = \frac {3}{4} \frac {\sqrt {-6Ro}}{DW} \frac {k_r^3}{\textit{Re}_z} + O \left(\textit{Re}_z^{-3} \right)\!, \end{equation}

as $ \textit{Re}_z \to \infty$ , which confirms the numerical result from Ali & Weidman (1993) that $ k_z \sim {1}/{\textit{Re}_z}$ in this limit; cf. figure 11 of their work. Unlike Ali & Weidman (1993) however, we have explicitly derived the coefficient of $ \textit{Re}_z^{-1}$ , showing its dependence on the flow properties. Indeed, substituting $ Ro$ evaluated at the mean geometric radius for $ \mu = 0$ from (A3) to (5.26), we obtain

(5.27) \begin{equation} k_z = \frac {3\sqrt {6}}{4DW\sqrt {1 - \eta }} \frac {k_r^3}{\textit{Re}_z} + O\big(\textit{Re}_z^{-3}\big). \end{equation}

Thus, we confirm that the coefficient of $\textit{Re}_z^{-1}$ increases as $\eta \to 1$ , in agreement with the numerical findings of Ali & Weidman (1993); cf. equation 9 and table II of their work.

• Parameterisation by azimuthal wavenumber $m$ : utilising the relationship (4.3) between the azimuthal and axial wavenumbers, with $m = r k_{\varphi }$ , we can interpret (5.5) as defining a family of curves parameterised by the azimuthal wavenumber $m \in {\mathbb{R}}$ . Although in this parametrisation the individual neutral stability curves (5.5) in the $(\textit{Re}_z, \textit{Re})$ plane look different (cf. figures 2 and 3), they have the same envelope (5.9), as verified either through direct calculation of the envelope or by minimising the critical Reynolds number in (5.10) with respect to $m$ . Specifically, considering the relation (5.3), we find that the neutral stability curves (5.5), parameterised by $m$ , exhibit horizontal asymptotes in the $(\textit{Re}_z,\textit{Re})$ plane as $\textit{Re}_z \to \infty$ :

  (5.28) \begin{equation} \textit{Re}_{\infty }(m) = \frac {\big(k_r^2 r^2 + m^2 \big)^{3/2}}{2|m|\varOmega r^2 \sqrt {-Ro}}. \end{equation}

Figure 3. Isothermal enclosed SCF with $\eta = 0.4$ and $k_r = 2\sqrt {2}$ . (a,c) Greenish centrifugal LELS instability domains in the $(\textit{Re}_z, \textit{Re})$ plane, parameterised by $m$ , and their envelope (5.9) (thick red curves) for (a) $\mu = 0$ and $\textit{Re}_{\infty } \approx 18.78$ , and (c) $\mu = 0.5$ and $\textit{Re}_{\infty } \approx 17.71$ . (b) The asymptotic value $\textit{Re}_\infty (m)$  (5.28) for individual neutral stability curves when $\mu = 0.5$ , showing a minimum equal to $\textit{Re}_{\infty } \approx 17.71$ (green circle), given by (5.20b ), at $m = m_{\textit{min}} \approx 2.108$ , as given by (5.29) or (5.30). (d) Variation of $m_{\textit{min}}$ with $\eta$ for $\mu = 0.5$ , according to (5.30), where the green circle indicates the minimal value $m \approx 2.108$ at $\eta = 0.4$ .

In figure 3(a,c) the greenish domains of centrifugal instability bounded by the solid black neutral stability curves (5.5), parameterised by $m$ , are shown in the $(\textit{Re}_z, \textit{Re})$ plane for the isothermal enclosed SCF with $W(r)$ given by (A17) and (A18), $\eta = 0.4$ and $k_r = 2\sqrt {2}$ for various values of $m$ and $\mu$ . The neutral stability curves are flattening as $\textit{Re}_z \to \infty$ , thus confirming the asymptotes (5.28). The oblique black solid lines in figure 3(c) indicate the inviscid LELS criterion (5.14). Vertical solid lines in figure 3(c) show $\pm \textit{Re}_z^{\textit{min}}\approx \pm 27.43$ given by (5.21). All computations for the SCF are performed at the mean geometric radius (A4).

By differentiating (5.28) with respect to $m$ , we determine the minimiser of $\textit{Re}_{\infty }(m)$ , i.e.

(5.29) \begin{equation} m_{\textit{min}} = \pm \frac {\sqrt {2}}{2} rk_r, \end{equation}

which thus defines the terminal azimuthal wavenumber as the floor of $m_{{min}}$ , $m_{\textit{term}}=\pm \left \lfloor |m_{\textit{min}}|\right \rfloor \in {\mathbb{Z}}$ , previously known only through numerical computations (Ali & Weidman 1993) (cf. figure 5 and table II of their work) and Meseguer & Marques (2000) (cf. figure 2 of their work). Substituting (5.29) into (5.28) yields the minimal value of $\textit{Re}_{\infty }(m)$ , exactly matching the limiting value $\textit{Re}_{\infty }$ in (5.20b ); see figure 3(a,b).

At the mean geometric radius (A4), (5.29) takes the form

(5.30) \begin{equation} m_{\textit{min}} = \pm \frac {k_r}{2}\frac {\sqrt {2\eta }}{1 - \eta }, \end{equation}

allowing us to find the terminal azimuthal wavenumber $m_{\textit{term}}=\pm \left \lfloor |m_{\textit{min}}|\right \rfloor$ at different ratios of the cylinder radii, $\eta$ ; see figure 3(d). Comparing figure 3(d) with figure 13 from Ali & Weidman (1993), one can see full qualitative and very good quantitative agreement between our analytical formula (5.30) and their earlier numerical results.

For example, using the parameters from figure 3, namely $\eta = 0.4$ and $k_r = 2\sqrt {2}$ , in (5.30), we find that $m_{\textit{min}} \approx 2.108$ . This yields a terminal axial wavenumber of $m_{\textit{term}} = 2$ . In figure 3 the neutral stability curve for azimuthal wavenumber $m = 2$ is nearly indistinguishable from the envelope at large $\textit{Re}_z$ , unlike the curves for $m = 1, 3$ and $5$ . This is because the terminal Reynolds number $\textit{Re}_{\infty }(2) \approx 17.74$ , as given by (5.28), is very close to the horizontal asymptote of the envelope $\textit{Re}_{\infty } \approx 17.71$ for the Rayleigh-stable case with $\mu = 0.5$ (figure 3 c), and $\textit{Re}_{\infty } \approx 18.78$ for the Rayleigh-unstable case with $\mu = 0$ (figure 3 a), as given by (5.20b ). The horizontal asymptotes $\textit{Re}=\textit{Re}_{\infty }$ are shown in figure 3(a,c) by black dot-dashed lines.

Explicit expression for the neutral stability surface in $(\textit{Re}_z, \textit{Re}_2, \textit{Re})$ space: note that the envelope equation (5.9) is universal for all isothermal swirling flows. For example, for flows between differentially rotating cylinders, the following relationship between the ratio of the angular velocities of the cylinders and the ratio of their radii, defined in (2.1), holds, i.e.

(5.31a,b,c) \begin{align} \mu = \eta \frac {\textit{Re}_2}{\textit{Re}}, \quad \textit{Re} = \frac {R_1 \varOmega _1 d}{\nu }, \quad \textit{Re}_2 = \frac {R_2 \varOmega _2 d}{\nu }, \end{align}

where $\textit{Re}$ and $\textit{Re}_2$ are the Reynolds numbers of the inner and outer cylinders, respectively (Meseguer & Marques 2000, 2002). This reformulates the azimuthal velocity (A2) and the expressions (A3) evaluated at the mean geometric radius (A4), which incorporate the parameter $\mu$ , in terms of $\textit{Re}_2$ . As a result, (5.9) yields the explicit analytical expression for the neutral stability surface within the $(\textit{Re}_z, \textit{Re}_2, \textit{Re})$ space:

(5.32) \begin{equation} \left [\frac {(1 + \eta )DW_g\textit{Re}_z}{2(\textit{Re}-\eta \textit{Re}_2)}\right ]^2 +\frac {27k_r^4(1 + \eta )^{2} - 16(1 - \eta )(\textit{Re} + \textit{Re}_2)(\eta \textit{Re} - \textit{Re}_2)} {27k_r^4(1 + \eta )^{2} - 16(1 - \eta )(\textit{Re} + \textit{Re}_2)(\textit{Re}-\eta \textit{Re}_2)}=0. \end{equation}

To the best of our knowledge, such a closed-form equation has not previously been reported in the literature on isothermal swirling flows, where numerical results are prevailing.

Note that the ‘individuality’ of an isothermal swirling flow enters (5.32) only by means of the radial derivative, $DW_g$ , of its axial velocity evaluated at the mean geometric radius (A4). The expressions for $DW_g$ derived in Appendix A are presented in table 1 for the SPF (A29), as well as for both open (A22) and enclosed (A21) SCFs.

In the case of the open SCF the expression for $DW_g$ given by (A22) is especially simple: $DW_g=(1-\eta)/(\sqrt{\eta}\ln{\eta})$ . The surface (5.32) and its cross-sections are shown in figure 4 for the open SCF with $k_r = \pi$ and $\eta = 0.5$ to facilitate comparison with the numerical results of Meseguer & Marques (2000). All computations for the open SCF are performed at the mean geometric radius (A4).

Figure 4. Open SCF with $ DW_g$ given by (A22), with $ k_r = \pi$ and $ \eta = 0.5$ . (a) Folded surface of the envelope (5.32) in $(\textit{Re}_z, \textit{Re}_2, \textit{Re})$ space. (b) For $ \textit{Re}_2 = 200$ , the thick solid red line shows the envelope (5.32) in the $(\textit{Re}_z, \textit{Re})$ plane. Shaded areas indicate individual instability domains defined by (5.5) and (4.3), corresponding to $ m = 0.07, 1, 2, 3$ . (c) Cross-sections of the instability domain (shaded regions) defined by the envelope surface (5.32) in the $(\textit{Re}_2, \textit{Re})$ plane, covering a range of $ \textit{Re}_z$ values from 0 to 55. (d) Shaded cross-sections of the surface (5.32) indicating instability domains at various $ \textit{Re}$ values in the $(\textit{Re}_z, \textit{Re}_2)$ plane that, for $\textit{Re}\gt 51.28$ , touch the thick red curve (5.40) that features a cusp point at $(43.79, 48.21)$ shown by the open circle.

At $\textit{Re}_z = 0$ , the surface yields the neutral stability curve of the Couette–Taylor flow obtained by the same local stability analysis by Eckhardt & Yao (1995), i.e.

(5.33) \begin{equation} (\textit{Re} + \textit{Re}_2)(\eta \textit{Re} - \textit{Re}_2)=\frac {27}{16}\frac {(1 + \eta )^{2}}{1-\eta }k_r^4, \end{equation}

which yields, at large $\textit{Re}_2$ ,

(5.34) \begin{equation} \textit{Re}=\frac {1}{\eta }\textit{Re}_2 + \frac {27}{16}\frac {1+\eta }{1-\eta }\frac {k_r^4}{\textit{Re}_2} + O\big(\textit{Re}_2^{-3}\big),\end{equation}

meaning that (5.33) has the inviscid Rayleigh line

(5.35) \begin{equation} \textit{Re}=\eta ^{-1}\textit{Re}_2 \end{equation}

as its asymptote shown by the thick dashed straight line in figure 4(c).

In another limit, $\textit{Re}_z \rightarrow \infty$ , the surface (5.32) produces a different neutral stability curve,

(5.36) \begin{equation} (\textit{Re} + \textit{Re}_2)(\textit{Re}-\eta \textit{Re}_2)=\frac {27}{16}\frac {(1 + \eta )^{2}}{1-\eta }k_r^4, \end{equation}

which yields, at large $\textit{Re}_2$ ,

(5.37) \begin{equation} \textit{Re}=\eta \textit{Re}_2 + \frac {27}{16}\frac {1 + \eta }{1-\eta }\frac {k_r^4}{\textit{Re}_2} + O\big(\textit{Re}_2^{-3}\big), \end{equation}

providing the line of solid-body rotation

(5.38) \begin{equation} \textit{Re}=\eta \textit{Re}_2 \end{equation}

as the asymptote to the curve (5.36), shown by the thick solid straight line in figure 4(c).

Cross-sections of the instability domains (shaded regions), defined by the envelope surface (5.32) in the $(\textit{Re}_2, \textit{Re})$ plane and shown in figure 4(c), span a range of $\textit{Re}_z$ values from 0 to 55. They demonstrate that the axial velocity component destabilises the Rayleigh-stable base flow and elucidate, in detail, the transition of the asymptotes to the stability boundaries from the Rayleigh line to the solid-body rotation line, as illustrated in figure 4(c).

A similar phenomenon occurs in the standard magnetorotational instability of magnetohydrodynamics, known as the Velikhov–Chandrasekhar paradox (Balbus & Hawley 1991). In this case, the stability boundary of magnetised Couette–Taylor flow drops below the Rayleigh line and approaches the solid-body rotation line when an axial magnetic field is applied, provided that the magnetic Prandtl number (the ratio of kinematic viscosity to magnetic diffusivity) differs from unity (Willis & Barenghi 2002; Kirillov et al. 2011; Kirillov & Stefani 2011).

Folds, pleat and cusp of the envelope surface for the open SCF: the transition between the Rayleigh line and the solid-body rotation line arises from the two folds of the surface (5.32), visible in figure 4(a), which converge at the pleat located at the point $(43.79, 48.21, 51.38)$ in the $(\textit{Re}_z, \textit{Re}_2, \textit{Re})$ space. These folds are projected onto a thick red curve in the $(\textit{Re}_z, \textit{Re}_2)$ plane, which exhibits a cusp at $(43.79, 48.21)$ , highlighted in figure 4(d) by the open circle. The open circle at $(48.21, 51.38)$ in figure 4(c) marks the projection of the pleat onto the $(\textit{Re}_2, \textit{Re})$ plane. The neutral stability surface with a pleat for the open SCF was first computed numerically by Meseguer & Marques (2000); our model now provides an explicit equation (5.32), allowing for the analytical computation of this surface’s projections, pleat and cusp points.

Indeed, differentiating (5.32) with respect to $ \textit{Re}$ and setting the result to zero yields the projection of the surface folds onto the $(\textit{Re}_2, \textit{Re})$ plane:

(5.39) \begin{eqnarray} & 729(1 + \eta )^4k_r^8 + 432(1 - \eta )(1 + \eta )^{2}(\textit{Re} + \textit{Re}_2)\big(\textit{Re}_2\eta ^{2} - 2\textit{Re}\eta + \textit{Re}_2 \big)k_r^4& \nonumber \\ &- 128(1 - \eta )^{2}(\textit{Re} + \textit{Re}_2)^2(\textit{Re}-\textit{Re}_2\eta ) \big(\textit{Re}_2\eta ^{2} - 2\textit{Re}\eta + \textit{Re}_2 \big)=0.& \end{eqnarray}

Interpreting (5.32) as a family of curves in the $(\textit{Re}_2, \textit{Re}_z)$ plane parameterised by $ \textit{Re}$ , we compute its discriminant as a polynomial in $ \textit{Re}$ to obtain the folds and pleat projection onto the $(\textit{Re}_2, \textit{Re}_z)$ plane for the open SCF flow (A22):

(5.40) \begin{eqnarray} &\def\lumina\hspace*{-1.5pc} 9\left (k_r^4 + \frac {2\textit{Re}_2^2}{27}\right )\textit{Re}_z^8 + 126(\ln 2)^2\textit{Re}_2^2\left (k_r^4 + \frac {5\textit{Re}_2^2}{63}\right )\textit{Re}_z^6&\nonumber \\ &\def\lumina\hspace*{-1.5pc} +\,2(\ln 2)^4\left(27k_r^4-4\textit{Re}_2^2\right)\left (243k_r^8 - \textit{Re}_2^2\left(27k_r^4 + 4\textit{Re}_2^2\right)\right )\textit{Re}_z^4&\nonumber \\ &\def\lumina\hspace*{-1.2pc} -\, 118098(\ln 2)^6\textit{Re}_2^2\left (k_r^4 + \frac {4\textit{Re}_2^2}{27}\right )^2\left (k_r^4 + \frac {4\textit{Re}_2^2}{243}\right )\textit{Re}_z^2 \nonumber\\ &\def\lumina\hspace*{-1.5pc}+\, 4782969(\ln 2)^8\left (k_r^4 + \frac {4\textit{Re}_2^2}{27}\right )^4 k_r^4= 0.& \end{eqnarray}

By differentiating (5.39) with respect to $ \textit{Re}$ , solving the resulting cubic equation and substituting the roots back into (5.39), we find $ \textit{Re}_2$ at the cusp (and pleat) point for $ \eta = 0.5$ , i.e.

(5.41) \begin{equation} \textit{Re}_2=\frac {3\sqrt {(48\sqrt {3} + 207)^{2/3} - (48\sqrt {3} + 207)^{1/3} + 33}}{2(48\sqrt {3} + 207)^{1/6}}k_r^2\approx 4.88k_r^2, \end{equation}

which, for $ k_r = \pi$ , gives $ \textit{Re}_2 \approx 48.21$ . Then, from (5.40), $ \textit{Re}_z \approx 43.79$ , and from (5.39), $ \textit{Re} \approx 51.38$ , locating the cusp and pleat points for the open SCF with $ k_r = \pi$ and $ \eta = 0.5$ as shown in figure 4.

In addition, figure 4(b) shows, in the $(\textit{Re}_z, \textit{Re})$ plane, the greenish individual instability domains defined by (5.5) and (4.3) for azimuthal wavenumbers $m = 0.07, 1, 2, 3$ . The boundaries of these domains are tangent to the envelope (5.32) shown as a thick red line, but intersect in such a way that point $A$ at $(133.67, 299.53)$ lies at the intersection of the neutral stability curves for $m = 1$ and $m = 2$ , while point $B$ at $(122.59, 174.45)$ lies at the intersection of those for $m = 2$ and $m = 3$ . This implies that the $m = 1$ mode is critical before point $A$ ; $m = 2$ is critical between points $A$ and $B$ ; and $m = 3$ becomes critical after point $B$ , as illustrated in figure 4(b). The azimuthal mode $m = 3$ is terminal, as indicated by (5.30) and consistent with (Meseguer & Marques 2000). The instability domain corresponding to $m = 0.07$ is also shown; as $m \rightarrow 0$ , this domain reduces to the semi-infinite interval $[\textit{Re}_0, \infty )$ , where $\textit{Re}_0$ is given by (5.20a ).

Figure 5. The BCF with $W(r)$ given by (A12) and (A13), and $\eta =0.8$ , $\textit{Pr}=5.5$ , $\gamma =0.0004$ and $k_r = 4\sqrt {2}$ the neutral stability curves $a_0=0$ with $a_0$ given by (4.9d ) in the $(Gr, \textit{Re})$ plane, parameterised (a–c) by $k_z$ or (d–f) by $m$ for the three different $\mu$ : (a,d) Rayleigh-unstable flow, $\mu = 0$ ; (b,e) modified Rayleigh line, $\mu = \mu _R \approx 0.63935$ (from (5.44)); and (c,f) Rayleigh-stable flow, $\mu = 0.8$ . The thick red curves show their envelope that has a horizontal asymptote $\textit{Re}=\textit{Re}_{\infty }$ for $m\gt 0$ with (a,d) $\textit{Re}_{\infty }\approx 24.66$ , (b,e) $\textit{Re}_{\infty }\approx 55.13$ and (c,f) $\textit{Re}_{\infty }\approx 83.99$ ; cf. Figure 7. In (a,d) the two envelope branches intersect at $\textit{Re}=\textit{Re}_0$ , where $\textit{Re}_0\approx 187.14$ is given by (5.42). The oblique black solid lines in (c,f) indicate the new unified visco-thermodiffusive extension (5.46) of the LELS criterion. All computations for the BCF are performed at the mean geometric radius (A4).

5.2. Visco-diffusive non-isothermal swirling flows – BCF

5.2.1. Stationary centrifugal instability as a combination of GSF and LELS instabilities

In general, $a_0=0$ with $a_0$ given by (4.9d ), defines a family of marginal stability curves for the BCF in the $(Gr, \textit{Re})$ plane, parameterised by the axial wavenumber $ k_z$ (or equivalently $ m$ via (4.3)), as shown in figure 5.

Envelope of the family of centrifugal instability domains: the individual marginal stability curves of the BCF have an envelope determined by computing the discriminant of the expression (4.9d ), treated as a polynomial in $k_z$ or $m$ (Hartman & Wintner 1953; Bruce & Giblin 1992; Hull 2020).

Modified Rayleigh line: the envelope, independent of the chosen parameterisation, separates the domain of unstable modes from the stability zone. It consists of two distinct symmetric curves that can intersect at $\textit{Gr} = 0$ and $\textit{Re}=\textit{Re}_0$ , as shown in figure 5(a), where

(5.42) \begin{equation} \textit{Re}_0 = \frac {3\sqrt {3}k_r^2}{2\varOmega }\frac {1}{\sqrt {4(Ro_R-Ro)(1-\gamma \varTheta )}} \end{equation}

and $\textit{Ro}_R$ is the modified Rayleigh line for non-isothermal flows:

(5.43) \begin{equation} Ro_R = -1+\frac {rD\varTheta \gamma \textit{Pr} }{4\left (1{-}\gamma \varTheta \right )}. \end{equation}

In the isothermal case $(\varTheta \equiv 0)$ , $\textit{Ro}_R=-1$ , and the expression (5.42) reduces to (5.20a ). Moreover, rewriting (5.42) in terms of the Taylor number $Ta = ( {2\textit{Re} \varOmega }/{3\sqrt {3}k_r^2})$ exactly reproduces the result obtained in (Kirillov & Mutabazi 2017) for CCF with a radial temperature gradient.

At the geometric mean radius (A4) the modified Rayleigh line (5.43) takes a convenient form for Couette–Taylor applications:

(5.44) \begin{equation} \mu _R=\eta ^{2}\frac {2(\gamma - 2)\ln \eta + \left(1-\frac {1}{\eta }\right)\gamma \textit{Pr}}{2(\gamma - 2)\ln \eta + (1 - \eta )\gamma \textit{Pr}}. \end{equation}

This relation shows that due to the temperature gradient, $\mu _R$ can substantially deviate from the isothermal Rayleigh line $\mu ={\eta }^2$ .

Rayleigh-unstable flow and intersection of envelope branches at $Gr=0$ : in the Rayleigh-unstable regime ( $\textit{Ro} \lt Ro_R$ or $\mu \lt \mu _R$ ), the envelope in the $(Gr,\textit{Re})$ plane can be approximated by its tangent lines computed at the intersection point $(Gr,\textit{Re})=(0,\textit{Re}_0)$ , i.e.

(5.45) \begin{equation} \frac {\textit{Re}}{\textit{Re}_0} = 1 \pm \sqrt {2} \frac {2D\varTheta \textit{Pr} \textit{Ro} - 3\gamma \varTheta DW k_r^2}{27 \textit{Ro} k_r^4} Gr, \end{equation}

where $\textit{Re}_0$ is given by (5.42). This result demonstrates that the temperature gradient plays a pivotal role in the existence of two distinct branches of the envelope for non-isothermal flows. These branches merge into a single curve in the case of isothermal Rayleigh-unstable flows, as established in the previous section; see figures 2 and 3.

In the absence of vertical gravity, as in the case of SPF with a radial temperature gradient (SPFRT) considered by Cotrell & McFadden (2005), the term $2D\varTheta \textit{Pr} Ro$ is absent in (5.45), and the splitting is primarily governed by the parameter $\gamma$ , which is often very small in practical situations (Kirillov 2025). As shown in Kirillov (2025), when $|\gamma | \ll \textit{Pr}$ , the envelope splitting becomes negligible, allowing for a zeroth-order approximation in $\gamma$ , which yields a smooth curve closely matching the exact envelope.

The condition $ a_0 \lt 0$ defines the instability regions for each $ k_z$ (or $ m$ ). The union of these regions for $ k_r = 4\sqrt {2}$ and $ k_z = \pm 0.2$ , $ \pm 0.8$ , $ \pm 4$ and $ \pm 7.3$ is represented as the green shaded area in figures 5(a) and 5(b). The boundaries of the individual instability domains, corresponding to the terminal values of the axial wavenumber, touch the envelope at $ \textit{Gr} = 0$ and $ \textit{Re} = \textit{Re}_0$ . These terminal values are $ k_z^{\textit{term}} = 4$ (dashed) and $ -k_z^{\textit{term}}=-4$ (solid), as determined by (5.23), which follows directly from (4.9d ) evaluated at $ \textit{Gr} = 0$ and $ \textit{Re} = \textit{Re}_0$ , where $\textit{Re}_0$ is defined by (5.42).

For $|k_z|\gt k_z^{\textit{term}}$ , the instability regions lie entirely within the green shaded area, as shown by dashed lines for $ k_z = 7.3$ and by solid lines for $ k_z = -7.3$ in figure 5(a). In figure 5(c) the dashed and solid curves mark the boundaries of the individual instability domains (their union is shown as a shaded green area) for $ k_z = 1$ and $ 3.5$ , and $ k_z = -1$ and $ -3.5$ , respectively.

While intersecting for $\mu \lt \mu _R$ , the envelope branches have vertical asymptotes at the modified Rayleigh line (5.44), delimiting the zone where no instability modes can be obtained (figures 5 band 5 e).

Rayleigh-stable flow and unified LELS-GSF centrifugal instability criterion: in the Rayleigh-stable regime ( $\textit{Ro} \gt Ro_R$ or $\mu \gt \mu _R$ ), the asymptotes to the envelope are inclined, placing the instability domains within the half-planes (Kirillov & Mutabazi 2024)

(5.46) \begin{equation} \frac {N_{\varOmega }^2}{\varOmega ^2}(1-\gamma \varTheta )+\textit{Pr}\frac {N^2}{\varOmega ^2} \lt \frac {\left ( \overline {DW}(2-\gamma \varTheta )-\frac {\textit{Pr} D\varTheta }{2k_r^2\varOmega S}\right )^2}{1 + \overline {DW}^2}, \end{equation}

defined by the asymptotes to the external branch, shown by oblique solid lines in figure 5(c, f):

(5.47) \begin{equation} \textit{Re} = \pm \frac {DW}{2\varOmega Ro} \sqrt {\frac {\varOmega ^2\left (2-\gamma \varTheta -\frac {\textit{Pr} D\varTheta Ro}{k_r^2DW}\right )^2-N_{\varOmega }^2(1-\gamma \varTheta )-\textit{Pr}N^2} {N_{\varOmega }^2(1-\gamma \varTheta )+\textit{Pr}N^2}}Gr. \end{equation}

Here $N_{\varOmega }^2$ is the Rayleigh discriminant (5.8), $S={\textit{Re}}/{Gr}$ and

(5.48) \begin{equation} N^2= -\gamma \varOmega ^2 r D\varTheta \end{equation}

is the square of the centrifugal Brunt–Väisälä frequency (Kirillov & Mutabazi 2017).

The inequality (5.46) yields the inviscid LELS criterion (5.14) in the isothermal case ( $\gamma =0$ , $D\varTheta =0$ ) (Leibovich & Stewartson 1983), and the GSF criterion

(5.49) \begin{equation} N_{\varOmega }^2(1-\gamma \varTheta )+\textit{Pr} N^2 \lt 0 \end{equation}

in the limit of azimuthal flow $S\rightarrow \infty$ (Maeder et al. 2013; Kirillov & Mutabazi 2017; Dymott et al. 2023). Thus, (5.46) presents a novel unified LELS-GSF centrifugal instability criterion for viscous and thermodiffusive swirling flows, incorporating the effects of Prandtl number ( $\textit{Pr}$ ) and temperature variation ( $D\varTheta$ ).

Note that a similar criterion for an inviscid isothermal incompressible swirling flow with radially varying density follows from the result by Eckhoff & Storesletten (1978) for compressible fluids, as shown by Di Pierro & Abid (2010); see Appendix D for details.

For completeness, we also present the equation for an asymptote to the internal branch of the envelope in the Rayleigh-stable case:

(5.50) \begin{equation} {\frac {N_{\varOmega }^2}{\varOmega ^2}(1-\gamma \varTheta ) + \textit{Pr} \frac {N^2}{\varOmega ^2} = 4(1-\gamma \varTheta )\left (\overline {DW} - \frac {\textit{Pr} D\varTheta }{2k_r^2\varOmega S}\right )\frac {\overline {DW}}{1+\overline {DW}^2}. } \end{equation}

The envelope in figure 5 closely matches the critical states curve from both numerical linear stability analysis and experiments (Lepiller et al. 2008; Yoshikawa et al. 2013; Guillerm et al. 2015; Kang et al. 2015, 2023), for example, cf. our figure 5(a,d) and figure 1(a) in (Yoshikawa et al. 2013). The envelope, unlike individual neutral stability curves, has a horizontal asymptote at $\textit{Re} = \textit{Re}_{\infty }$ as $|Gr| \rightarrow \infty$ . This explains the seemingly smooth and nearly $Gr$ -independent stability boundary observed in Lepiller et al. 2008, Yoshikawa et al. 2013, Guillerm et al. 2015 and Kang et al. 2015, 2023 for $\mu =0$ , although this boundary is actually piecewise smooth, with each neutral stability curve touching the common envelope, which flattens at large $|Gr|$ where shear instability dominates. Despite the Rayleigh–Fjörtoft shear instability mechanism due to an inflection point in the axial velocity profile (Drazin 2002) leading to axisymmetric perturbations (Bahloul et al. 2000; Lepiller et al. 2007), rotation ensures that the critical modes of the BCF are three-dimensional with $k_z\ne 0$ and $m\ne 0$ (Dubrulle et al. 2005).

Selection of unstable modes through envelope branches: figure 5, in conjunction with relation (4.3) demonstrates how the envelope selects the critical modes. The lower left branch of the envelope corresponds to the left spiral modes ( $k_z \gt 0$ , $m \gt 0$ ), while the lower right branch corresponds to the right spiral modes ( $k_z \lt 0$ , $m \gt 0$ ). This is consistent with numerical findings by Ali & Weidman (1990). The modes with $m\lt 0$ touch only the upper parts of the envelope and define the instability domains lying entirely within the green shaded area in figure 5, corresponding to the union of the instability regions with different values of $m$ .

Parameterisation  by $m$ and   determination  of  horizontal  asymptotes  to  envelope branches: when parameterised by $m$ , the neutral stability curves $a_0=0$ with $a_0$ given by (4.9d ), corresponding to different values of $m=k_{\varphi } r$ have horizontal (not depending on $Gr$ ) asymptotes in the $(Gr,\textit{Re})$ plane; see figure 5. The exact asymptotic values of the Reynolds number $\textit{Re}_{\infty }(m)$ for a particular $m$ are determined by the quadratic equation, following from $a_0=0$ in the limit $|Gr| \rightarrow \infty$ :

(5.51) \begin{align} q(m,\textit{Re}_{\infty })&:= 4\varOmega ^2m^2r^4(1 \!-\!\gamma \varTheta )\nonumber\\&\quad\times \big[\big(k_r^2r^2 + m^2\big)(1+ \textit{Ro}_R-Ro)DW + r^2 \textit{Ro} D\varTheta \textit{Pr}\big]\textit{Re}_{\infty }^2 \nonumber \\ &\quad -\, 2\varOmega m k_r r^3 \big(k_r^2 r^2 + m^2\big)\big[r^2 \textit{Ro} D\varTheta \textit{Pr} - \gamma \varTheta DW \big(k_r^2r^2 + m^2\big)\big]\textit{Re}_{\infty } \nonumber \\ &\quad -\, DW\big(k_r^2r^2 + m^2\big)^4=0. \end{align}

In the isothermal case $(\gamma =0, D\varTheta =0)$ , $\textit{Ro}_R=-1$ and (5.51) reduces to

(5.52) \begin{equation} -DW\big(k_r^2r^2 + m^2\big)\big[4\textit{Ro} \varOmega ^2 m^2r^4 \textit{Re}_{\infty }^2 + \big(k_r^2r^2 + m^2\big)^3\big]=0, \end{equation}

which yields expression (5.28) determining $\textit{Re}_{\infty }(m)$ for the isothermal flows.

As $m \to \infty$ , the positive root of (5.51) asymptotically behaves as

(5.53) \begin{equation} \textit{Re}_{\infty }(m) \sim \frac { m^2}{2\varOmega r^2\sqrt {(1 + \textit{Ro}_R-Ro)(1-\gamma \varTheta )}}. \end{equation}

In contrast, as $m \to 0$ , it follows that $\textit{Re}_{\infty }(m) \sim m^{-1}$ . This behaviour indicates the existence of a minimum value of $\textit{Re}_{\infty }$ at some finite value of $m$ .

Solving simultaneously the equation $q(m,\textit{Re}_{\infty })=0$ given by (5.51) and $\partial _m q(m,\textit{Re}_{\infty })=0$ , we can find the pairs of $m_{\textit{min}}$ and $\textit{Re}_{\infty }(m_{\textit{min}})$ determining the minimiser and the minimal value of $\textit{Re}_{\infty }$ , which is the horizontal asymptote to the envelope in the non-isothermal case. For instance, choosing the parameters of figure 5 corresponding to the BCF with $\eta =0.8$ , $\textit{Pr}=5.5$ , $\gamma =0.0004$ and $k_r=4\sqrt {2}$ , evaluating (5.51) at the mean geometric radius (A4) and applying the procedure just described, we find the pairs shown in figure 6.

Figure 6. The BCF with $ W(r)$ defined by (A12) and (A13), and parameters $\eta = 0.8$ , $\textit{Pr} = 5.5$ , $\gamma = 0.0004$ and $k_r = 4\sqrt {2}$ , (red thick line) the asymptotic values of $\textit{Re}$ as $|Gr| \to \infty$ are shown for the neutral stability curves parameterised by $m$ , based on (5.51), for the following cases: (a) Rayleigh-unstable flow, $\mu = 0$ ; (b) modified Rayleigh line, $\mu = \mu _R \approx 0.63935$ (from (5.44)); and (c) Rayleigh-stable flow, $\mu = 0.8$ . The minimum of $\textit{Re}_{\infty }(m)$ for $m \lt 0$ corresponds to the horizontal asymptote of the upper branch of the envelope, while the minimum for $m \gt 0$ corresponds to the horizontal asymptote of the lower branch of the envelope, as illustrated in figure 5. The thin green parabola represents the approximation from (5.53).

5.2.2. Visco-thermodiffusive oscillatory McIntyre instability

The critical Reynolds number and the modified Rayleigh line for purely azimuthal flow: in the absence of axial flow ( $DW = 0$ ), we consider the equation $ a_1a_2 = a_0$ , derived from the stability condition (5.1c ), as a polynomial in $ k_z$ . Computing its discriminant provides the envelope of the neutral stability curves for oscillatory instability in the $ (\textit{Pr}, \textit{Re})$ plane, expressed as $ \textit{Re} = \textit{Re}_0^{\prime}$ , where

(5.54) \begin{equation} { \textit{Re}_0^{\prime} =\frac {3\sqrt {3}k_r^2}{2\varOmega } \frac {1+\textit{Pr} }{ \textit{Pr}} \frac {1}{\sqrt {4(Ro_{R}^{\prime}-Ro)(1-\gamma \varTheta )}} } \end{equation}

and

(5.55) \begin{equation} { Ro_{R}^{\prime}=-1 + \frac {1+\textit{Pr}}{8\textit{Pr}} \frac {\,r\,\gamma \, D\varTheta }{ 1-\gamma \varTheta } } \end{equation}

defines the modified Rayleigh line. Notably, $ Ro_{R}^{\prime} = Ro_R$ and $ \textit{Re}_{0}^{\prime} = 2\textit{Re}_0$ at $ \textit{Pr} = 1$ .

Rewriting (5.54) in terms of the Taylor number $ Ta = ( {2\textit{Re} \varOmega }/{3\sqrt {3}k_r^2})$ exactly reproduces the result of Kirillov & Mutabazi (2017), obtained for CCF with a radial temperature gradient.

At the geometric mean radius (A4), the expression (5.55) takes a form well suited for Couette–Taylor applications:

(5.56) \begin{equation} \mu _{R}^{\prime} = \eta ^{2}\frac {4\textit{Pr}(\gamma - 2)\ln \eta + \left(1-\frac {1}{\eta }\right)\gamma (\textit{Pr} + 1)}{4\textit{Pr} (\gamma - 2)\ln \eta + (1 - \eta )\gamma (\textit{Pr} + 1)}. \end{equation}

Similar to (5.44), the presence of a temperature gradient can cause $ \mu _{R}^{\prime}$ to deviate significantly from the isothermal Rayleigh line $ \mu ={\eta }^2$ .

Figure 7. Green-shaded regions represent stationary (LELS-GSF) instability domains touching the thick blue envelope, while purple-shaded regions above them correspond to oscillatory (McIntyre) instability domains touching the thick red envelope. Dashed black lines indicate neutral stability curves for $k_z \gt 0$ , and solid black lines indicate neutral stability curves for $k_z \lt 0$ . (a) Rayleigh-unstable BCF ( $\eta = 0.8$ , $\mu = 0$ ) with radial wavenumber $k_r = 4\sqrt {2}$ , Prandtl number $\textit{Pr} = 5.5$ , outward heating ( $\gamma = 0.0004$ ) and $k_z = \pm 0.6, \pm 2, \pm 4$ ; (b) Rayleigh-stable BCF ( $\eta = 0.8$ , $\mu = 0.8$ ) with $k_r = 3\sqrt {2}$ , $\textit{Pr} = 5.5$ , $\gamma = 0.01$ and $k_z = \pm 0.1, \pm 0.5, \pm 0.6, \pm 1.5, \pm 3$ ; (c) Rayleigh-unstable BCF ( $\eta = 0.8$ , $\mu = 0.62$ ) with $k_r = 4\sqrt {2}$ , $\textit{Pr} = 4$ , $\gamma = 0.01$ and $k_z = \pm 0.6, \pm 2, \pm 4$ ; (d) Rayleigh-stable BCF with a quasi-Keplerian azimuthal velocity profile ( $\eta = 0.99$ , $\mu = \eta ^{3/2}$ ), $k_r = 2\sqrt {2}$ , $\textit{Pr} = 0.5$ , inward heating ( $\gamma = -0.01$ ) and $k_z = \pm 0.2, \pm 0.5, \pm 1, \pm 2, \pm 3$ .

The  envelope  of  oscillatory  instability  domains  for  visco-thermodiffusive  swirling flows: for $ DW \neq 0$ , we plot the neutral stability curves $ a_1a_2 - a_0 = 0$ , which bound the oscillatory instability domains in the $ (Gr, \textit{Re})$ plane for different fixed values of $ k_z$ . These we show alongside the neutral stability curves $ a_0 = 0$ , which bound the stationary instability domains in the same plane.

Figure 7 presents the results of this computation for a range of Rayleigh-unstable and Rayleigh-stable BCF. Both the stationary and oscillatory instability domains form families parameterised by $ k_z$ , with their envelopes depicted as thick red curves (oscillatory instability) and thick blue curves (stationary instability). Each envelope consists of two branches, which may intersect at $ \textit{Gr} = 0$ and $ \textit{Re} = \textit{Re}_0$ for the stationary instability and at $ \textit{Gr} = 0$ and $ \textit{Re} = \textit{Re}_{0}^{\prime}$ for the oscillatory instability. Here, $ \textit{Re}_0$ is given by (5.42) and $ \textit{Re}_{0}^{\prime}$ is defined by (5.54). Recall that (5.42) indicates that the intersection point $ (\textit{Gr} = 0, \textit{Re} = \textit{Re}_0)$ exists if $ \textit{Ro} \lt Ro_R$ (or $ \mu \lt \mu _R$ ), where $ Ro_R$ (or $ \mu _R$ ) is determined from (5.43) (or (5.44)).

Similarly, the intersection point $ (\textit{Gr} = 0, \textit{Re} = \textit{Re}_{0}^{\prime})$ of the branches forming the envelope of the oscillatory instability domains exists, as follows from (5.54), if $ \textit{Ro} \lt Ro_{R}^{\prime}$ (or $ \mu \lt \mu _{R}^{\prime}$ ), where $ Ro_{R}^{\prime}$ (or $ \mu _{R}^{\prime}$ ) is determined from (5.55) (or (5.56)). At this intersection point the linear approximation to the branches is given by

(5.57) \begin{align} \frac {\textit{Re}}{\textit{Re}_{0}^{\prime}} &= 1 \pm \frac {\sqrt {2}}{243} \nonumber\\&\quad\times\frac {\gamma ^2 \left ( \textit{Pr} \!-\! 1 \right ) \textit{Re}_{0}^{{\prime}{2}} \varOmega ^2 r\varTheta D\varTheta DW \!+\! 9k_r^{2} \left ( 1+\textit{Pr} \right ) \left ( RoD\varTheta \!-\! 3 \gamma \varTheta k_r^{2} DW \right )}{ \textit{Ro} k_r^{6} \left (1+ \textit{Pr}\right )^2 } PrGr. \end{align}

Note that (5.57) presents an analytical expression for the onset of a visco-diffusive McIntyre oscillatory instability in a visco-thermodiffusive swirling flow with a radial temperature gradient and natural gravity, valid for an arbitrary Prandtl number, which, to the best of our knowledge, has not previously been reported in the literature.

By comparing the values of $ \textit{Re}_0$ and $ \textit{Re}_{0}^{\prime}$ , we can predict which type of instability will define the critical $ \textit{Re}$ near $ \textit{Gr} = 0$ . For the Rayleigh-unstable flow with an outward heating ( $\gamma \gt 0$ ) in figure 7(a), we find that $ \textit{Re}_0 \approx 187.14$ and $ \textit{Re}_{0}^{\prime} \approx 221.10$ , indicating that the flow is more likely to experience stationary instability. This is confirmed by the domain of oscillatory instability, shown in purple in figure 7(a), being contained entirely within the green-shaded domain of stationary instability. This result agrees with the numerical and experimental results by Yoshikawa et al. (2013) and Guillerm et al. (2015), respectively.

In contrast, the Rayleigh-unstable BCF with an outward heating in figure 7(c) exhibits $ \textit{Re}_{0}^{\prime} \approx 1049.11 \lt \textit{Re}_0 \approx 1326.27$ , indicating that oscillatory instability dominates over stationary instability within a finite range of $ Gr$ near $ \textit{Gr} = 0$ . Notably, the red envelope of the oscillatory instability domains intersects with the blue envelope of the stationary instability domains, forming two codimension-2 points for $\textit{Gr} \lt 0$ and $\textit{Gr} \gt 0$ , as shown in figure 7(c). This suggests that a sufficiently large temperature gradient may favour the stationary modes.

For the Rayleigh-stable BCF with an outward heating in figure 7(b), the stationary and oscillatory instability domains are distinctly separated. Specifically, for the parameters used in figure 7(b), $ \mu = 0.8 \gt \mu _{R}^{\prime} \approx 0.6381$ ( $ \mu = 0.8 \gt \mu _R \approx 0.6224$ ), indicating that $ \textit{Ro} \gt Ro_{R}^{\prime}$ in (5.54) and $ \textit{Ro} \gt Ro_R$ in (5.42). Consequently, the self-intersection points for both envelopes do not exist at $ \textit{Gr} = 0$ . The domain of oscillatory instability lies entirely within the domain of stationary instability, with the latter being dominant for the BCF shown in figure 7(b).

Finally, the Rayleigh-stable quasi-Keplerian BCF with inward heating ( $ \gamma \lt 0$ ) shown in figure 7(d) is characterised by the stationary instability domain being split into two parts, as $ \mu \approx 0.9850 \gt \mu _R \approx 0.9825$ . In contrast, the envelope of the oscillatory instability domains exhibits a self-intersection point at $ \textit{Gr} = 0$ , as $ \mu \approx 0.9850 \lt \mu _{R}^{\prime} \approx 0.9874$ . As a result, oscillatory instability dominates over stationary instability within a finite range near $ \textit{Gr} = 0$ , and the envelopes of the oscillatory and stationary instability domains intersect, forming two codimension-2 points, as shown in figure 7(d). This behaviour qualitatively agrees with findings by Cotrell & McFadden (2005) for SPF with a radial temperature gradient; see also (Kirillov 2025).

Figure 8. (a) Growth rates and (b) frequencies of stationary (green curves) and oscillatory McIntyre (purple curves) instabilities for the Rayleigh-stable BCF with a quasi-Keplerian azimuthal velocity profile ( $\eta = 0.99$ , $\mu = \eta ^{3/2}$ ), $k_r = 2\sqrt {2}$ , $\textit{Pr} = 0.5$ , inward heating ( $\gamma = -0.01$ ), $Gr=500$ and $k_z =-1$ .

Figure 8 illustrates the behaviour of the growth rate $\textrm{Re}(\lambda )$ and frequencies ${\textrm{Im}}(\lambda )$ for the Rayleigh-stable BCF with a quasi-Keplerian azimuthal velocity profile ( $\eta = 0.99$ , $\mu = \eta ^{3/2}$ ), $k_r = 2\sqrt {2}$ , $\textit{Pr} = 0.5$ , inward heating ( $\gamma = -0.01$ ), $\textit{Gr} = 500$ and $k_z = -1$ , corresponding to the stability map shown in figure 7(d).

Figure 8(a) shows that, as the azimuthal Reynolds number increases, the green curve corresponding to ${\textrm{Im}}(\lambda ) = 0$ rises from negative values, crosses the line $\textrm{Re}(\lambda ) = 0$ at some critical $\textit{Re}$ , becomes positive (indicating instability), reaches a maximum and then decreases again to negative values (restoring stability).

This interval of instability corresponds to the crossing of the green centrifugal instability domain in figure 7(d), which has both lower and upper critical Reynolds numbers. According to (4.6), $\textrm{Re} (\lambda ) = \sigma$ , i.e. the growth rate of the instability, and ${\textrm{Im}}(\lambda) = 0$ corresponds to the frequency of the unstable mode given by $\omega = -m\varOmega - ( {k_zW}/{S})$ . Thus, the marginal mode is stationary in a reference frame attached to the helical streamlines around the annular axis (Kirillov & Mutabazi 2024).

With a further increase in $\textit{Re}$ , the green curve in figure 8(a) attains increasingly negative values until it meets another branch, forming a double negative real root ( $\textrm{Re}(\lambda ) \lt 0$ , ${\textrm{Im}}(\lambda ) = 0$ ). This root then splits into two complex conjugate $\lambda$ roots with $\textrm{Re}(\lambda ) \lt 0$ , which, upon a continued increase in $\textit{Re}$ , eventually intersect the line $\textrm{Re}(\lambda ) = 0$ , as shown in figure 8.

Hence, the unstable mode acquires a non-trivial frequency, in contrast to the critical mode of the centrifugal instability, whose frequency is merely the Doppler shift. The growth rate and frequency of the complex conjugate pair are shown as purple curves in figure 8, and it is evident that the oscillatory instability is weaker than the stationary one, as its growth rate is an order of magnitude smaller, which is typical for dissipation-induced instabilities (McIntyre 1970; Labarbe & Kirillov 2021; Kirillov 2021). The growth rate of the oscillatory instability remains positive as long as the system remains within the purple domain of oscillatory instability; see figure 7(d).

As established in the earlier work by McIntyre (1970), the visco-diffusive oscillatory instability is significantly enhanced when the Prandtl number deviates from unity ( $\textit{Pr} \ne 1$ ). In figure 9 a reduction of the Prandtl number to $\textit{Pr} = 0.03$ reveals a richer and more intricate geometry of the envelope delineating the neutral stability boundaries of oscillatory instability domains, shown in purple. The two thick red branches of this envelope, which define the critical values of the control parameters $\textit{Re}$ and $Gr$ , intersect not only at $\textit{Gr} = 0$ and $\textit{Re} = \textit{Re}^{\prime}_{0}$ as in figure 7(a,c), but also at two additional non-zero values of $Gr$ , symmetrically located with respect to the vertical axis.

Figure 9. Green-shaded regions represent stationary (LELS-GSF) instability domains touching the thick blue envelope, while purple-shaded regions above them correspond to oscillatory (McIntyre) instability domains touching the thick red envelope for the Rayleigh-stable BCF with $\eta = 0.8$ , $k_r = 3\sqrt {2}$ , $\textit{Pr} = 0.03$ , inward heating ( $\gamma = -0.01$ ) and (a) $\mu = 0.692$ , $k_z = \pm 0.04, \pm 1.5, \pm 3, \pm 10, \pm 20$ ; and (b) $\mu = 0.8$ , $k_z = \pm 0.02, \pm 1.5, \pm 3, \pm 5, \pm 15$ . Dashed black lines indicate neutral stability curves for $k_z \gt 0$ and solid black lines indicate neutral stability curves for $k_z \lt 0$ .

This behaviour is illustrated in figure 9(a), corresponding to a Rayleigh-stable BCF configuration with $\eta = 0.8$ , $k_r = 3\sqrt {2}$ , $\textit{Pr} = 0.03$ , inward heating ( $\gamma = -0.01$ ) and $\mu = 0.692$ . Moreover, the thick red branches corresponding to oscillatory instability domains intersect the thick blue branches associated with stationary centrifugal instability domains depicted in greenish hues. As a result, the white stability domain enclosed by these envelopes takes on a distinctive diamond-like shape with five codimension-2 points (figure 9(a)) – a phenomenon that, to the best of our knowledge, has not been previously reported in the literature.

Increasing $\mu$ from 0.692 to 0.8 in figure 9(b), while keeping all other parameters identical to those in figure 9(a), reveals that the two self-intersection points of the thick red branches – defining the envelope of the oscillatory instability domains – persist even though the central intersection at $\textit{Gr} = 0$ and $\textit{Re} = \textit{Re}^{\prime}_{0}$ disappears. In figure 9(b) the two thick red branches exhibit distinct asymptotic directions as $\textit{Re} \rightarrow \infty$ . This indicates that one branch governs the stability boundary for lower values of $\textit{Re}$ before the codimension-2 point, while the other dominates after the intersection, in the limit $\textit{Re} \rightarrow \infty$ . This asymptotic divergence complicates the formulation of an analytical stability criterion analogous to the LELS-GSF criterion for centrifugal instability (5.46), making the development of a corresponding theoretical framework significantly more challenging.

Indeed, the equation $a_1 a_2 = a_0$ , derived from the stability condition (5.1c ) and governing the neutral stability curves that bound the domains of oscillatory instability, can be expressed as a polynomial in $k_z$ :

(5.58) \begin{align} c_1k_z^3\textit{Re}^3 &+ k_z^2\big(c_2k_z^4 + c_3k_z^2 + c_4\big)\textit{Re}^2 + k_z\big(c_5k_z^6 + c_6k_z^4 + c_7k_z^2 + c_8\big)\textit{Re}\nonumber \\ &+ c_9 k_z^{10} + c_{10}k_z^8 + c_{11} k_z^6 + c_{12}k_z^4 + c_{13}k_z^2 + c_{14}=0. \end{align}

By computing the discriminant of the polynomial (5.58), we find that in the Rayleigh-stable regime the asymptotic directions of the envelope bounding the domains of oscillatory instability in the $(Gr, \textit{Re})$ plane, as $\textit{Re} \rightarrow \infty$ , are determined by the roots of the equation in terms of the swirl parameter $S = {\textit{Re}}/{Gr}$ , i.e.

(5.59) \begin{equation} 6912c_9^3c_1^4c_2^7\big(27c_1^2c_{14}^2 - 18c_1c_4c_8c_{14} + 4c_1c_8^3 + 4c_4^3c_{14} - c_4^2c_8^2\big)=0, \end{equation}

where

(5.60) \begin{align} c_1 &= 2\overline {DW}\varOmega ^3\textit{Pr}^2S\varTheta \gamma ^2k_r r \big[4\overline {DW}^2\varOmega ^2\big(1-\gamma \varTheta\big)\nonumber\\&\quad- \big(1+\overline {DW}^2\big)\big(N_{\varOmega }^2(1-\gamma \varTheta ) + N^2\big)\big],\nonumber \\ c_2 &= -\,\varOmega ^2 \textit{Pr} S \gamma r \big(1+\overline {DW}^2\big)^2\nonumber \\ &\quad\times \big[\big(1+\overline {DW}^2\big)\big(2 \textit{Pr}(1-\gamma \varTheta )N_{\varOmega }^2 + (\textit{Pr} + 1)N^2\big) - 8\overline {DW}^2\varOmega ^2 \textit{Pr}(1-\gamma \varTheta )],\nonumber \\ c_4 &= -\varOmega \textit{Pr} k_r^2 \{\varOmega S \gamma k_r^2\big(2\textit{Pr}\big(1+\overline {DW}^2\big)\big(1-\gamma \varTheta\big)N_{\varOmega }^2\nonumber \\ &\quad+ \big[4\overline {DW}^2\varOmega ^2\varTheta ^2\gamma ^2 + (1+\overline {DW}^2)N^2\big]\big(\textit{Pr} + 1\big) - 8\overline {DW}^2\varOmega ^2 \textit{Pr} (1-\gamma \varTheta )\big)r \nonumber\\&\quad+ 2N^2\overline {DW}\textit{Pr}\},\nonumber \\ c_8 &= -\,k_r^5\big(\textit{Pr} + 1\big)\big(2S\varOmega ^3\varTheta \gamma ^2\overline {DW}k_r^2\big(3\textit{Pr} + 1\big)r + N^2\textit{Pr}\big),\nonumber \\ c_9 &= -\,2S\big(\textit{Pr} + 1\big)^2\big(1+\overline {DW}^2\big)^5 r\gamma \varOmega ^2,\nonumber \\ c_{14} &= -\,2S\big(\textit{Pr} + 1\big)^2 k_r^{10}r\gamma \varOmega ^2. \end{align}

The prefactors in (5.59) yield the conditions $c_1 = 0$ and $c_2 = 0$ , which correspond to possible asymptotic directions. These lead to the expressions

(5.61) \begin{equation} { \frac {N^2}{\varOmega ^2} +(1-\gamma \varTheta )\left (\frac {N_{\varOmega }^2}{\varOmega ^2}- \frac {4\overline {DW}^2}{1+\overline {DW}^2}\right ) =0 } \end{equation}

and

(5.62) \begin{equation} { \frac {\textit{Pr} + 1}{2\textit{Pr}}\frac {N^2}{\varOmega ^2}+(1-\gamma \varTheta )\left (\frac {N_{\varOmega }^2}{\varOmega ^2} - \frac {4\overline {DW}^2}{1+\overline {DW}^2}\right ) =0, } \end{equation}

respectively. Additional possible asymptotic directions arise from the vanishing of the bracketed term in (5.59). Which of these asymptotic directions is realised in a given regime depends on the specific parameter values entering the coefficients defined in (5.60).

This complex structure simplifies considerably in the asymptotic limit $\textit{Pr} \rightarrow 0$ , relevant for astrophysical and geophysical contexts. In particular, such values of $\textit{Pr}$ are typical in solar and stellar convection zones, where $\textit{Pr} \sim 10^{-6}$ , as well as in certain liquid metals, where $\textit{Pr} \sim 10^{-3}$ (Käpylä 2021). In this limit, the envelope of the oscillatory instability domains reduces to a single curve in the Rayleigh-unstable case, and to two symmetric, non-intersecting curves in the Rayleigh-stable case. Under these conditions, the bracketed term in (5.59) determines the unique asymptotic line as $\textit{Re} \rightarrow \infty$ , yielding the following criterion for the onset of the oscillatory McIntyre instability:

(5.63) \begin{equation} \frac {N_{\varOmega }^2}{\varOmega ^2}-\frac {4\overline {DW}^2}{1+\overline {DW}^2}\lt \frac {1+\overline {DW}^2}{16\overline {DW}^2(1-\gamma \varTheta )\gamma ^2\varTheta ^2}\left [\frac {N^2}{\varOmega ^2}-\frac {4\overline {DW}^2}{1+\overline {DW}^2}\gamma ^2\varTheta ^2\right ]^2. \end{equation}

6. Conclusion

Adapting the geometrical optics method to visco-thermodiffusive flows, we have developed a unified analytical framework for the local linear short-wavelength instabilities of both isothermal and non-isothermal swirling flows.

Within this framework, we showed that isothermal viscous swirling flows support only stationary centrifugal instabilities. We derived a closed-form expression for the growth rate and an explicit analytical stability criterion that generalises the classical LELS criterion by incorporating viscosity. This result recovers known limits, including the Rayleigh and Eckhardt–Yao criteria, and provides a theoretical basis for the observed transition between Rayleigh-stable and Rayleigh-unstable regimes.

A key innovation of our work lies in identifying and analytically constructing the envelope of the neutral stability curves using the connection between envelopes and polynomial discriminants. This geometric insight enables precise determination of terminal wavenumbers and critical states relevant to both experiments and numerical simulations.

Extending the analysis to non-isothermal swirling flows, we uncovered the existence of both stationary instability – arising from a combination of the LELS and GSF instabilities – and an oscillatory instability. The latter, a visco-thermodiffusive generalisation of the McIntyre instability, is governed by a new analytical criterion unifying the effects of radial stratification, rotation and diffusion. We further developed an algorithm to predict the dominant instability mode as a function of flow parameters, particularly the Prandtl number.

Our results align well with previous experimental and numerical studies of spiral Couette, spiral Poiseuille and baroclinic Couette flows, and provide a theoretical tool for guiding future investigations and experimental design.