1. Introduction
Eddies, along with fronts, i.e. regions of relatively sharp horizontal density gradients, are common occurrences in the atmosphere and the ocean. In the atmosphere, baroclinic and barotropic instabilities result in the formation of cyclones and anticyclones (Emanuel, Fantini & Thorpe Reference Emanuel, Fantini and Thorpe1987), and the resulting advection of temperature leads to the formation of fronts (Hoskins & Bretherton Reference Hoskins and Bretherton1972). Similar processes occur in the ocean too. For example, tropical instability near the equator can result in the formation of tropical instability vortices which are associated with strong equatorial sea surface temperature fronts (Holmes et al. Reference Holmes, Thomas, Thompson and Darr2014). Alternatively, processes such as wind and topographical variations can lead to the development of submesoscale eddies such as the structures which are associated with frontal regions (Mahadevan Reference Mahadevan2006). In addition, stirring by mesoscale (around 
$100$ km or more, corresponding to a Rossby number much smaller than unity) and submesoscale (1–100 km, corresponding to a Rossby number around unity) eddies can form frontal regions too (Shcherbina et al. Reference Shcherbina2015; Sarkar et al. Reference Sarkar, Pham, Ramachandran, Nash, Tandon, Buckley, Lotliker and Omand2016). The dynamics of these eddies, apart from significantly influencing the transport of various quantities such as heat, salinity and nutrients (Zhang, Wang & Qiu Reference Zhang, Wang and Qiu2014), can also have significant implications for large-scale phenomena like the El Niño southern oscillation and the equatorial heat budget (Holmes et al. Reference Holmes, Thomas, Thompson and Darr2014). Furthermore, the dynamics of submesoscale coherent vortices is an important consideration to understand transport and mixing in many parts of the ocean (McWilliams Reference McWilliams1985).
Large-scale geophysical eddies, apart from their vortical structure, often have a vertical (out of plane) gradient of momentum which is in thermal wind balance with the horizontal density variation (Pedlosky Reference Pedlosky1987). Factors such as background rotation and vertical stratification also affect the dynamics of these vortices. In addition, curvature effects have also been shown to have an effect on frontogenesis and the adjustment of fronts near vortices (Shakespeare Reference Shakespeare2016). Instabilities in these vortices represent an important pathway towards small-scale turbulence and mixing (Thomas, Tandon & Mahadevan Reference Thomas, Tandon and Mahadevan2008), thus playing a key role in parameterization of small-scale processes in large-scale climate models (Fox-Kemper, Ferrari & Hallberg Reference Fox-Kemper, Ferrari and Hallberg2008; Bachman et al. Reference Bachman, Fox-Kemper, Taylor and Thomas2017). In addition, their instability characteristics can also help in interpreting field observations (Ruddick Reference Ruddick1992; D'Asaro et al. Reference D'Asaro, Lee, Rainville, Harcourt and Thomas2011; Thompson et al. Reference Thompson, Lazar, Buckingham, Garabato, Damerell and Heywood2016) and large-scale climate model outputs (Fox-Kemper et al. Reference Fox-Kemper, Danabasoglu, Ferrari, Griffies, Hallberg, Holland, Maltrud, Peacock and Samuels2011), which in turn motivate fine-scale process studies (Mahadevan & Tandon Reference Mahadevan and Tandon2006).
Instabilities in geophysical vortices could range from purely inertial, such as the centrifugal instability in unstratified planar vortices (Kloosterziel & Van Heijst Reference Kloosterziel and Van Heijst1991), to purely convective, resulting from statically unstable vertical stratification (Chandrasekhar Reference Chandrasekhar1961). The intermediate regimes, where both inertial and buoyancy effects are present, can be conducive to what is called the symmetric instability (Hoskins Reference Hoskins1974; Xui & Clark Reference Xui and Clark1985; Thomas et al. Reference Thomas, Taylor, Ferrari and Joyce2013). Purely centrifugal instabilities are generally suppressed by Coriolis forces (Godeferd, Cambon & Leblanc Reference Godeferd, Cambon and Leblanc2001) but could play an important role in regions with anticyclonic vorticity and near the equator, where Rossby numbers are near unity (Haine & Marshall Reference Haine and Marshall1998). Symmetric instability can result in slantwise convection along tilted constant angular momentum surfaces when the potential vorticity is negative in an inviscid, adiabatic environment (Hoskins Reference Hoskins1974). Recently, curvature effects on the symmetric instability have been studied, and their importance in oceanic eddies has been highlighted (Buckingham, Gula & Carton Reference Buckingham, Gula and Carton2021a,Reference Buckingham, Gula and Cartonb). In addition, hyperbolic instabilities could also manifest in regions of high strain near the edges of vortices with hyperbolic stagnation points (Leblanc Reference Leblanc1997). Apart from the aforementioned local instabilities, global instabilities such as the Kelvin–Helmholtz (Peltier & Caulfield Reference Peltier and Caulfield2003) and baroclinic instabilities (Pierrehumbert & Swanson Reference Pierrehumbert and Swanson1995) are also a consideration in stratified eddies.
Vortices in thermal wind balance can also become unstable due to unequal mass and momentum diffusivities. McIntyre (Reference McIntyre1970), using normal mode theory in a locally Cartesian region of an eddy in thermal wind balance, showed that a Prandtl number away from unity modifies the symmetric instability criterion, and also introduces a local oscillatory instability. Laboratory experiments have demonstrated the occurrence of such an oscillatory instability, and resulting density layers (Baker Reference Baker1971; Meunier et al. Reference Meunier, Miquel, Le Dizès, Chowdhury and Alam2014). Studies regarding instabilities in long-lived mesoscale oceanic eddies have considered this diffusive instability mechanism to explain intrusions found in density profiles (Ruddick Reference Ruddick1992; Kuzmina & Zhurbas Reference Kuzmina and Zhurbas2000). Experimental studies have also observed the formation of density layers at a Prandtl number less than unity, which was further suggested to explain the layers on the scale of 
$1\ \text {m}$ seen in certain parts of the ocean. (Calman Reference Calman1977; Weber Reference Weber1980). An energetics-based approach to derive the results of McIntyre (Reference McIntyre1970) in a parallel shear flow was presented by Kloosterziel & Carnevale (Reference Kloosterziel and Carnevale2007). Here, we use the local stability approach (Lifschitz & Hameiri Reference Lifschitz and Hameiri1991) with small diffusivities in momentum and density to investigate short-wavelength instabilities in stratified eddies in thermal wind balance. Previous studies using the local stability approach have captured the effects of differential diffusion (in mass and momentum) on local instabilities in radially stratified circular Couette flow (Kirillov & Mutabazi Reference Kirillov and Mutabazi2017) and vertically stratified planar vortices (Singh & Mathur Reference Singh and Mathur2019).
In this study, we explore local instabilities in horizontally and vertically stratified vortices in thermal wind balance, including the effects of differential diffusion in mass and momentum. The local stability approach is first demonstrated to be efficient in capturing symmetric and double-diffusive instabilities. Apart from being based on a different approach compared with McIntyre (Reference McIntyre1970), our study is also wider in scope compared with that of McIntyre (Reference McIntyre1970) in three significant aspects: (i) exploration in the centrifugally unstable regime, (ii) viscous instability characteristics away from the neutral stability boundaries in both the centrifugally stable and unstable regimes and (iii) curvature effects. The local stability framework and the base flow configuration are described in § 2 followed by our results on the instability characteristics in § 3. The importance of our results, particularly the curvature effects, are discussed in three different contexts in § 4, following which the conclusions are presented in § 5.
2. Theory
The governing equations of motion for an incompressible flow with background rotation 
$\varOmega {\boldsymbol {e}}_z$ (
${\boldsymbol {e}}_z$ is the unit vector along the positive 
$z$-axis), within the Boussinesq approximation, are
$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}$$
$$\begin{gather}\frac{\textrm{D}\boldsymbol{u}}{\textrm{D}t} + 2\varOmega{\boldsymbol{e}}_z\times\boldsymbol{u} ={-}\frac{\boldsymbol{\nabla} p}{\rho_0} - g\frac{\rho}{\rho_0}\boldsymbol{e}_z +\nu\nabla^2\boldsymbol{u}, \end{gather}$$
$$\begin{gather}\frac{\textrm{D}\rho}{\textrm{D}t} = \kappa\nabla^2\rho, \end{gather}$$
where 
$t$ is time, 
$\rho$, 
$\boldsymbol {u}$ and 
$p$ the density, velocity and pressure, respectively, and gravity 
$\boldsymbol {g} = -g\boldsymbol {e}_z$ with 
$g>0$; 
$\rho _0$, 
$\nu$ and 
$\kappa$ are the reference density, kinematic viscosity and diffusion coefficient associated with density, respectively.
2.1. Base flow
We consider a base flow field that describes an axisymmetric vortex in stratified surroundings. Specifically, the velocity field is given by
\begin{equation} \boldsymbol{u}_{{{B}}} = V(r,z)\boldsymbol{e}_\omega, \end{equation}
where 
$(r,\omega ,z)$ represent the cylindrical polar coordinates, and the subscript 
$B$ denotes base flow quantities. The corresponding density field is given by 
$\rho _{B}(r,z)$, with its vertical and radial gradients characterized by the Brunt Väisälä frequency 
$N_z = \sqrt {-(\partial \rho _{B}/\partial z)(g/\rho _0)}$ (
$N_z^2>0$ throughout the current study), and 
$N_r = \sqrt {-(\partial \rho _{B}/\partial r)(g/\rho _0)}$. Finally, the pressure associated with the base flow is denoted as 
$p_{B}(r,z)$. The base flow field satisfies the mass conservation (2.1), the azimuthal component of the inviscid momentum equation ((2.2) with 
$\nu =0$) and the diffusion-free incompressibility condition ((2.3) with 
$\kappa =0$). Furthermore, requiring the base flow to satisfy the radial component of (2.2) gives
\begin{equation} 2\left(\varOmega + \frac{{V}}{r}\right)\frac{\partial{V}}{\partial z} ={-}\frac{g}{\rho_0}\frac{\partial \rho_{B}}{\partial r} (= N_r^2), \end{equation}
where it has been assumed that the vertical component of (2.2) is described by hydrostatic balance. In other words, (2.5) represents a modified thermal wind relation (Shakespeare Reference Shakespeare2016), where the centripetal, Coriolis and pressure gradient forces are in balance. In the limit of 
$V/r \to 0$, (2.5) recovers the well-known geostrophic balance (Vallis Reference Vallis2017). In this study, however, the curvature effects are retained in the relation between the vertical gradient of momentum and the horizontal stratification. A schematic of the base flow in thermal wind balance is shown in figure 1.
Figure 1. (a) A schematic of the axisymmetric base flow (
$\boldsymbol {u}_{\boldsymbol {B}} = V(r,z)\boldsymbol {{e}_\omega }$) considered in this study. The colour indicates density, and the black curves denote streamlines. The radial spacing between neighbouring streamlines is inversely proportional to the magnitude of the local azimuthal velocity. The flow is stably stratified along 
$z$. The radial density gradient and the vertical (along 
$z$) gradient of the azimuthal velocity are in thermal wind balance (2.5). (b) Schematic showing the perturbation wave vector 
$\boldsymbol {k}$ and the angle 
$\theta$ it makes with the 
$z$-axis.
2.2. Local stability equations
To study the stability properties of the base flow in § 2.1, we superimpose small perturbations to write the net flow field as
\begin{equation} \boldsymbol{u} = \boldsymbol{u}_{{{B}}{}} + \boldsymbol{u'},\quad {\rho} = {\rho}_{{{B}}{}} + {\rho'},\quad {p} = {p}_{{{B}}{}} + {p'}, \end{equation}where the prime denotes perturbation quantities. The linearized, viscous governing equations for the perturbation flow field are
$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u'} = 0, \end{gather}$$
$$\begin{gather}\frac{\textrm{d}\boldsymbol{u}'}{\textrm{d}t} + \boldsymbol{\nabla} \boldsymbol{u}_{{{B}}{}}\boldsymbol{\cdot}\boldsymbol{u}' + 2\varOmega{\boldsymbol{e}}_z\times\boldsymbol{u}' + \boldsymbol{\nabla} p' = \nu\nabla^2\boldsymbol{u'} - b'e_z, \end{gather}$$
$$\begin{gather}\frac{\textrm{d} b'}{\textrm{d}t} + {\boldsymbol{u}}'\boldsymbol{\cdot}\boldsymbol{\nabla} {\rho}_{{{B}}{}} \frac{g}{\rho_o} = \kappa\nabla^2 b', \end{gather}$$
where 
$b' = g\rho '/\rho _0$, and the base flow has been assumed to satisfy the inviscid governing equations. In other words, (2.7)–(2.9) represent the governing equations for small, viscous perturbations superimposed on the inviscid base flow described in § 2.1. Here, 
$\textrm {d}/\textrm {d}t = \partial /\partial t + \boldsymbol {u}_B\boldsymbol {\cdot }\boldsymbol {\nabla }$ is the material time derivative in the base flow 
${\boldsymbol {u}}_B$.
In accordance with the local stability approach (Lifschitz & Hameiri Reference Lifschitz and Hameiri1991), we assume short-wavelength perturbations of the form
$$\begin{gather} \boldsymbol{u}' = \exp(\textrm{i}\phi({\boldsymbol{x}},t)/\delta)[\boldsymbol{a}({\boldsymbol{x}},t) + \delta \boldsymbol{a}_{\delta}({\boldsymbol{x}},t) + \ldots], \end{gather}$$
$$\begin{gather}p' = \exp(\textrm{i}\phi({\boldsymbol{x}},t)/\delta)[{\rm \pi}({\boldsymbol{x}},t) + \delta {\rm \pi}_{\delta}({\boldsymbol{x}},t) + \ldots], \end{gather}$$
$$\begin{gather}b' = \exp(\textrm{i}\phi({\boldsymbol{x}},t)/\delta)[b({\boldsymbol{x}},t) + \delta b_{\delta}({\boldsymbol{x}},t) + \ldots], \end{gather}$$
where 
$\delta$ is a small ‘short-wavelength’ parameter, and 
$\phi ({\boldsymbol {x}},t)$ is a real scalar field that contains the wave-like behaviour in the perturbations. The perturbation wave vector 
$\boldsymbol {k}'$ is given by 
$\boldsymbol {k}' = \boldsymbol {\nabla } \phi /\delta = \boldsymbol {k}/\delta$, thus rendering the perturbation to be of short wavelength since 
$\delta$ is a small parameter; (
$\boldsymbol {a},{\rm \pi} ,b$) are the complex leading-order amplitudes of the velocity, pressure and buoyancy perturbations, respectively. Similarly, 
$({\boldsymbol {a}}_\delta ,{\rm \pi} _\delta ,b_\delta )$ are the complex perturbation amplitudes at 
${O}(\delta )$.
Assuming the two diffusivities 
$\nu$ and 
$\kappa$ to be of 
${O}(\delta ^2)$, and substituting the perturbation forms in (2.10)–(2.12) in the governing equations (2.7)–(2.9), results at 
${O}(\delta ^0)$ in the local stability equations
\begin{gather} \frac{\textrm{d}\boldsymbol{k}}{\textrm{d}t} ={-}\boldsymbol{\nabla}\boldsymbol{u}_{{{B}}{}}^\textrm{T}\boldsymbol{\cdot}\boldsymbol{k},\end{gather}
\begin{gather} \frac{\textrm{d}\boldsymbol{a}}{\textrm{d}t} ={-}\boldsymbol{\nabla}\boldsymbol{u}_{{{B}}{}}\boldsymbol{\cdot}\boldsymbol{a} - 2\varOmega{\boldsymbol{e}}_z\times\boldsymbol{a} - b\boldsymbol{e}_z + \frac{\boldsymbol{k}}{|\boldsymbol{k}|^2}[2(\boldsymbol{\nabla}\boldsymbol{u}_{{{B}}{}}\boldsymbol{\cdot} \boldsymbol{a})\boldsymbol{\cdot}\boldsymbol{k} + 2(\varOmega{\boldsymbol{e}}_z\times\boldsymbol{a})\boldsymbol{\cdot}\boldsymbol{k} + bk_z] \nonumber\\ \hspace{-19.5pc}- c_\nu|\boldsymbol{k}|^2\boldsymbol{a}, \end{gather}
\begin{gather} \frac{\textrm{d}b}{\textrm{d}t} ={-}(\boldsymbol{a}\boldsymbol{\cdot}\boldsymbol{\nabla}\rho_{{{B}}{}})\frac{g}{\rho_0} - c_\kappa|\boldsymbol{k}|^2b.\end{gather}
In (2.13)–(2.15), 
$\textrm {d}/\textrm {d}t = \partial /\partial t + \boldsymbol {u}_B\boldsymbol {\cdot }\boldsymbol {\nabla }$ is the material time derivative in the base flow 
${\boldsymbol {u}}_B$, 
$k_z$ the 
$z$-component of 
$\boldsymbol {k}$, and 
$(c_\nu ,c_\kappa ) = (\nu ,\kappa )/\delta ^2$. The assumptions of 
$\nu$ and 
$\kappa$ being of 
${O}(\delta ^2)$ are motivated by the second-order spatial gradients in the viscous terms, which then contribute at 
${O}(\delta ^0)$ in (2.8)–(2.9). Previous studies on local instabilities have incorporated viscous effects in this manner. For example, (2.13)–(2.14) without background rotation and buoyancy were considered by Landman & Saffman (Reference Landman and Saffman1987), who showed that viscous terms suppress the growth rates associated with elliptic instability in an unstratified, planar elliptical vortex. In addition to elucidating viscous damping of inviscid instabilities, the 
${O}(\delta ^2)$ scaling of momentum and mass diffusivities have helped the local approach unearth double-diffusive instabilities in stratified vortices as well (Kirillov & Mutabazi Reference Kirillov and Mutabazi2017; Singh & Mathur Reference Singh and Mathur2019). Finally, the 
${O}(\delta ^2)$ scaling for temperature and salinity diffusivities also leads to the reduction of local stability equations to the classical salt fingering instability equation for a quiescent base flow (Appendix A).
In the absence of diffusive effects (
$c_\nu = c_\kappa = 0$), (2.13)–(2.15) coincide with the equations considered by Miyazaki (Reference Miyazaki1993) for a stably stratified elliptical vortex with no in-plane stratification; further, if background rotation is also absent, (2.13)–(2.15) reduce to those considered by Miyazaki & Fukumoto (Reference Miyazaki and Fukumoto1992). Finally, if background rotation and buoyancy are both absent, (2.13)–(2.14) coincide with those derived by Lifschitz & Hameiri (Reference Lifschitz and Hameiri1991). The equations governing the evolution of 
$\boldsymbol {a}$ and 
$b$ are solved on closed fluid particle trajectories in the base flow for periodic wave vectors, to subsequently obtain corresponding growth rates. We proceed to obtain the conditions under which perturbation wave vectors are periodic upon integrating the wave vector evolution equation (2.13) over a closed trajectory in the base flow given by (2.4). For a wave vector 
$\boldsymbol {k}= k_r\boldsymbol {e}_r + k_\omega \boldsymbol {e}_\omega + k_z\boldsymbol {e}_z$, (2.13) can be written in component form as
$$\begin{gather} \frac{\textrm{d}k_r}{\textrm{d}t} = k_\omega\frac{{V}}{r} - \frac{\partial {V}}{\partial r}k_\omega, \end{gather}$$
$$\begin{gather}\frac{\textrm{d}k_\omega}{\textrm{d}t} = 0, \end{gather}$$
$$\begin{gather}\frac{\textrm{d}k_z}{\textrm{d}t} ={-}\frac{\partial {V}}{\partial z}k_\omega. \end{gather}$$
Equation (2.17) indicates that 
$k_\omega$ is a constant along the trajectory, and a non-zero value for the constant would lead to a monotonically changing 
$k_r$ and 
$k_z$. Therefore, for each of 
$k_r$, 
$k_\omega$ and 
$k_z$ to be periodic upon integration over one closed trajectory, 
$k_\omega$ has to necessarily be zero. With 
$k_\omega = 0$, 
$k_r$ and 
$k_z$ become constant along the trajectory. It is noteworthy that 
$k_\omega = 0$ represents zero azimuthal wavenumber, i.e. axisymmetric perturbations in the normal mode approach, although the perturbation amplitudes are still allowed to evolve along the circular base flow trajectory in the local stability approach. In summary, any wave vector given by 
$\boldsymbol {k}= k_r\boldsymbol {e}_r + k_z\boldsymbol {e}_z$ represents a periodic wave vector; the amplitude evolution equations (2.14)–(2.15) are solved for such periodic wave vectors. The periodicity of the wave vector enables the use of Floquet theory to obtain growth rates, and ensures uniqueness of the perturbation wave vector at a given spatial location. The latter aspect is potentially relevant for establishing a connection between local and global stability approaches.
For the base flow in (2.4), considering periodic wave vectors, (2.14)–(2.15) can be written in cylindrical polar coordinates as
\begin{equation} \frac{\textrm{d}(\boldsymbol{a}|b)}{\textrm{d}t} = \boldsymbol{A}(\boldsymbol{a}|b), \end{equation}
where 
$(\boldsymbol {a}|b) = [a_r, a_\omega , a_z, b]^\textrm {T}$ and 
$\boldsymbol {A}$ is a 
$4\times 4$ matrix. Substituting 
$\boldsymbol {a} = a_r\boldsymbol {e}_r + a_\omega \boldsymbol {e}_\omega + a_z\boldsymbol {e}_z$ and a periodic wave vector 
$\boldsymbol {k} = k_r\boldsymbol {e}_r + k_z\boldsymbol {e}_z$ in (2.14)–(2.15), we obtain
\begin{align} {\boldsymbol{A}} =\left[\begin{array}{cccc} -c_\nu|\boldsymbol{k}|^2 & 2\left(\dfrac{{V}}{r}+\varOmega\right)\left(1-\dfrac{k_r^2}{|\boldsymbol{k}|^2}\right) & 0 & \dfrac{k_rk_z}{|\boldsymbol{k}|^2} \\ -\left(\dfrac{\partial {V}}{\partial r}+\dfrac{{V}}{r}+2\varOmega\right) & -c_\nu|\boldsymbol{k}|^2 & {- \dfrac{\partial {V}}{\partial z}} & 0 \\ 0 & - \dfrac{2k_zk_r}{|\boldsymbol{k}|^2}\left(\dfrac{{V}}{r}+\varOmega\right) & -c_\nu|\boldsymbol{k}|^2 & - \left(1-\dfrac{k_z^2}{|\boldsymbol{k}|^2}\right)\\ {N_r^2} & 0 & N_z^2 & -c_\kappa|\boldsymbol{k}|^2 \end{array}\right], \end{align} where 
$\textrm {d}\boldsymbol {e}_r/\textrm {d}t = \boldsymbol {e}_\omega V/r$ and 
$\textrm {d} \boldsymbol {e}_\omega /\textrm {d}t = -\boldsymbol {e}_r V/r$ have been used to account for circular fluid particle trajectories with an angular speed 
$V/r$. Owing to the constancy of 
$r$ and 
$V$ on a streamline, and 
$k_r$, 
$k_z$, 
$|\boldsymbol {k}|$ also being invariant along a streamline for periodic wave vectors, the coefficient matrix 
$\boldsymbol {A}$ in (2.20) is time independent.
3. Results
The eigenvalues (
$\lambda$) of 
$\boldsymbol {A}$ (2.20) represent the growth rates, and are governed by
\begin{align} &{(c_\nu|\boldsymbol{k}|^2+\lambda)}\left\{ (c_\nu|\boldsymbol{k}|^2+\lambda)^2(c_\kappa|\boldsymbol{k}|^2+\lambda) + (c_\nu|\boldsymbol{k}|^2+\lambda) \left[N_z^2\left(1-\frac{k_z^2}{|\boldsymbol{k}|^2}\right){-N_r^2\frac{k_rk_z}{|\boldsymbol{k}|^2}}\right]\right. \nonumber\\ &\quad + (c_\kappa|\boldsymbol{k}|^2+\lambda)\left[{-}2\frac{\partial {V}}{\partial z}\left(\frac{{V}}{r}+\varOmega\right)\frac{k_rk_z}{|\boldsymbol{k}|^2}\right.\nonumber\\ &\quad \left.\left. +2\left(\frac{\partial{V}}{\partial r}+\frac{{V}}{r}+2\varOmega\right)\left(\frac{{V}}{r}+\varOmega\right) \left(1-\frac{k_r^2}{|\boldsymbol{k}|^2}\right)\right] \right\} = 0. \end{align}
The growth rate 
$\sigma$ is given by the maximum real part of 
$\lambda$, which we analyse in detail in the rest of this paper. Instabilities with zero or non-zero imaginary part in the corresponding eigenvalues are henceforth referred to as monotonic or oscillatory, respectively. In the presence of diffusion terms 
$(c_\nu >0,c_\kappa >0)$, one of the eigenvalues of 
$\boldsymbol {A}$ is 
$-c_\nu |\boldsymbol {k}|^2$, while the other three eigenvalues of 
$\boldsymbol {A}$ can be given by a cubic equation. Dividing by 
$|\varPhi /2|^{3/2}$, (3.1) is non-dimensionalized to give the governing equation for the other three eigenvalues as
\begin{align}
&(\tilde{\lambda} + \tilde{c})^3 +
\tilde{c}\left(\frac{1-Sc}{Sc}\right)(\tilde{\lambda}+\tilde{c})^2
+ 2\hat{\varPhi}(\cos^2\theta +
\tan\varTheta\tan\varGamma\sin^2\theta \notag\\
&\quad
-\tan\varGamma\sin2\theta)(\tilde{\lambda}+\tilde{c})
+\tilde{c}\left(\frac{1-Sc}{Sc}\right)2\hat{\varPhi}\left(\cos^2\theta
- \frac{\tan\varGamma}{2}\sin2\theta\right) = 0,
\end{align}where
\begin{equation} \varPhi = 2\left(\frac{V}{r}+\varOmega\right)\left(\frac{\partial V}{\partial r} + \frac{V}{r} + 2\varOmega\right), \end{equation}
where 
$\tilde {\lambda } = \lambda /\sqrt {|\varPhi /2|}$, 
$\hat {\varPhi } = \varPhi /|\varPhi |$ and 
$Sc = \nu /\kappa = c_\nu /c_\kappa$ is the Schmidt number. As shown in § 3.2, 
$\hat {\varPhi } = -1$ corresponds to centrifugal instability in the inviscid limit. With this choice of non-dimensionalization, the important limit of weak (but finite) viscous effects is easily recovered without the divergence of non-dimensionalized growth rate.
The non-dimensional parameters that appear in (3.2) are
$$\begin{gather} {\tan{\varGamma} = \frac{N_r^2}{\varPhi}}, \quad \tan{\varTheta} = \frac{N_z^2}{N_r^2}, \end{gather}$$
$$\begin{gather}\theta = \cos^{{-}1}\frac{k_z}{|\boldsymbol{k}|}, \quad {\tilde{c} = \frac{c_\nu|\boldsymbol{k}|^2}{\sqrt{|\varPhi/2|}}} , \end{gather}$$
where 
$(\varGamma ,\varTheta )$ in (3.4a,b) describe the base flow, and 
$(\theta ,\tilde {c})$ in (3.5a,b) describe the perturbation wave vector. 
$\varPhi$, as defined in (3.3), is the Rayleigh discriminant, i.e. twice the product between the total angular velocity and the absolute vertical vorticity in the base flow. Specifically, 
$\varGamma$ is a measure of the vertical velocity gradient (recall that 
$N_r^2 = 2(\varOmega + {{V}}/{r}){\partial {V}}/{\partial z}$ due to thermal wind balance) with respect to the strength of centrifugal instability in the inviscid, unstratified limit (see § 3.2), and 
$\varTheta$ is the inclination of the isopycnals with respect to the 
$z$-axis in the base flow. In (3.5a,b), 
$\theta$ is the inclination of the perturbation wave vector with respect to the 
$z$-axis, and 
$\tilde {c}$ is the inverse of a perturbation Reynolds number, which quantifies the viscous effects in the perturbation with respect to the inertia in the base flow. It is worth noting that the non-dimensional parameters in (3.4a,b)–(3.5a,b) are defined such that the inviscid limit can be recovered without any singularities. In addition to the four parameters in (3.4a,b)–(3.5a,b), 
$Sc$ is also a parameter that appears in (3.2), thus influencing the instability characteristics. As shown in § 3.4, for given values of 
$(\tan \varGamma , \tan \varTheta , \theta , \tilde {c}, Sc)$, monotonic and oscillatory instabilities cannot occur simultaneously.
3.1. Comparison with McIntyre (Reference McIntyre1970)
McIntyre (Reference McIntyre1970) considered a base flow similar to the current study, described by a purely azimuthal velocity field 
$V_M(r,z)\boldsymbol {e}_\omega$ and a potential temperature field 
$T_M(r,z)$, under the influence of background rotation 
$\varOmega \boldsymbol {e}_z$, gravity 
$-g\boldsymbol {e}_z$ and in thermal wind balance. Defining 
$C_M(r,z) = 2\varOmega r + {{V}_M}(r,z)$, McIntyre (Reference McIntyre1970) identified the following non-dimensional base flow parameters:
\begin{equation} \tan{\varGamma}_M = \frac{{{C}}_{M,z}}{{{C}}_{M,r}}, \quad \tan{\varTheta}_M = \frac{{T}_{M,z}}{{T}_{M,r}}, \end{equation}
with the Schmidt number being given by 
$Sc = \nu /\kappa$. In (3.6a,b), the subscripts 
$z$ and 
$r$ denote differentiation with respect to 
$z$ and 
$r$, respectively.
Operating in the regime 
$N_z^2 = g\alpha T_{M,z}>0$ (statically stable, with 
$\alpha$ being the thermal expansion coefficient), 
$2\varOmega C_{M,r}>0$ (centrifugally stable in locally Cartesian coordinates), McIntyre (Reference McIntyre1970) considered a locally Cartesian (curvature effects ignored) unbounded normal mode theory (Wentzel–Kramers–Brillouin–Jeffreys, i.e. WKBJ-type approximation where spatial gradients of the base flow are assumed constant). Considering perturbations of the form 
$\exp [\textrm {i}k_M(r\sin \theta +z\cos \theta ) + \omega _M t]$, with 
$\theta$ being the angle between the wave vector and the 
$z$-axis, the author derived the governing equation for the non-dimensional growth rate as
\begin{equation} (\omega_M + k_M^2)^2(\omega_M + Sc^{{-}1}k_M^2) + (G + J)\omega_M + (G + Sc^{{-}1}J)k_M^2 = 0, \end{equation}
where 
$\textrm {Re}(\omega _M)$ is the growth rate, 
$k_M$ is the magnitude of the wave vector, 
$G = s(\tan {\varTheta _M}\sin ^2\theta - \sin \theta \cos \theta )$, 
$J = s(\cot {\varGamma _M}\cos ^2\theta - \sin \theta \cos \theta )$ and 
$s = \text {sgn}(\tan {\varGamma _M}) = \text {sgn}(\tan {\varTheta _M})$. The criteria for monotonic and oscillatory instabilities were then shown to be 
${\tan {\varTheta _M}}/{\tan {\varGamma _M}} < {(1+Sc)^2}/{4Sc}$ and 
${\tan {\varTheta _M}}/{\tan {\varGamma _M}} < {(3Sc+1)^2}/{8Sc(Sc+1)}$, respectively. These criteria can be recovered as the zero curvature limit of our results in § 3.3.3.
To explore the relation between our local stability framework, and McIntyre's theory, we re-write (3.2) in the limit of zero curvature and non-dimensionalization using the length and time scales of McIntyre (Reference McIntyre1970). Specifically, defining 
$L = (\nu ^2/|N_r^2| )^{1/4}$ and 
$T = (|N_r^2|)^{-1/2}$ as the length and time scales, (3.2) can be written as
\begin{align} &\left(\hat{\lambda}+\frac{{\hat{k}}^2}{\delta^2}\right)^2\left(\hat{\lambda}+Sc^{{-}1}\frac{{\hat{k}}^2}{\delta^2}\right) + s[\tan{\varTheta_M}\sin^2\theta + \cot{\varGamma_M}\cos^2\theta - 2\sin\theta\cos\theta]\hat{\lambda}\nonumber\\ &\quad + s[\tan{\varTheta_M}\sin^2\theta - \sin\theta\cos\theta+ Sc^{{-}1}(\cot{\varGamma_M}\cos^2\theta - \sin\theta\cos\theta)]\frac{{\hat{k}}^2}{\delta^2} = 0, \end{align}
where 
$\hat {\lambda } = \lambda T$, 
${\hat {k}} = \boldsymbol {k}L$ and 
$\nu = c_\nu \delta ^2$ is the kinematic viscosity as defined in § 2.2; 
$\hat {\varPhi } = 1$ has been assumed. We note that 
$\hat {\varPhi } = 1$ is equivalent to 
$2\varOmega C_{M,r}>0$, i.e. the study of McIntyre (Reference McIntyre1970) is restricted to the first and third quadrants of the 
$(\tan \varGamma ,\tan \varTheta )$ plane. Equation (3.8) is the same as (3.7), with 
$k_M = {\hat {k}}/\delta$. Thus, we have shown that, in the limit of zero curvature, the local stability equation for the growth rate concurs with the locally Cartesian unbounded normal mode theory. This connection between local stability and normal mode approaches builds on previous such connections demonstrated for other flows: inviscid, unstratified axisymmetric vortices with an axial flow (Le Duc & Leblanc Reference Le Duc and Leblanc1999; Mathur et al. Reference Mathur, Ortiz, Dubos and Chomaz2014), viscous unstratified elliptical vortex (Landman & Saffman Reference Landman and Saffman1987) and viscous, in-plane stratified axisymmetric flow with an arbitrary Prandtl number (Kirillov & Mutabazi Reference Kirillov and Mutabazi2017).
Apart from our study being based on a different approach compared with McIntyre (Reference McIntyre1970), our results are also wider in scope compared with the results of McIntyre (Reference McIntyre1970) in three significant aspects: (i) exploration in the centrifugally unstable regime (second and fourth quadrants of 
$(\tan \varGamma ,\tan \varTheta )$ plane), (ii) viscous instability characteristics away from the neutral stability boundaries in both the centrifugally stable and unstable regimes and (iii) curvature effects. With respect to point (ii), we note that the length scale 
$L$ used by McIntyre (Reference McIntyre1970) for non-dimensionalization makes it difficult to explore the regime of small but finite viscous effects. As discussed in § 3.3, significant growth rates of both monotonic and oscillatory instabilities can occur in the weakly viscous regime, where even the most unstable perturbations can lie. The curvature effects are discussed more in detail in § 4.1.
3.2. Inviscid limit
Substituting the inviscid limit of 
$c_\nu = c_\kappa = 0$ in (3.1), the non-trivial roots are obtained as
\begin{equation} \lambda^2 ={-}\frac{1}{2}(N_z^2 + \varPhi) + \frac{1}{2}(N_z^2 - \varPhi)\cos2\theta + \left(\frac{N_r^2}{2} + \frac{\partial V}{\partial z}\left(\frac{V}{r} + \varOmega\right)\right)\sin2\theta, \end{equation}
where 
$\theta \in [0,{\rm \pi} ]$ is the angle made by the wave vector 
$\boldsymbol {k}$ with the 
$z$-axis, such that 
$k_z/|\boldsymbol {k}| = \cos \theta$ and 
$k_r/|\boldsymbol {k}| = \sin \theta$. Maximizing 
$\lambda ^2$ with respect to 
$\theta$, and requiring it to be positive results in the following inviscid instability criterion:
\begin{equation} {N_z^2\varPhi} < {2{N_r^2\left(\frac{V}{r}+\varOmega\right)\frac{\partial V}{\partial z}} ( = N_r^4)}. \end{equation}
As a result, 
$\varPhi <0$ (in other words, 
$\hat {\varPhi }=-1$) guarantees instability (recall that 
$N_z^2>0$ is assumed). The maximum growth rate occurs at
\begin{equation}
\theta^* = \begin{cases} \dfrac{\rm \pi}{4} -
\frac{1}{2}\tan^{{-}1}{(N_z^2-{\varPhi})}/{2N_r^2}, &
\text{if } N_r^2\ge0 \\
\dfrac{3{\rm \pi}}{4} -
\frac{1}{2}\tan^{{-}1}{(N_z^2-{\varPhi})}/{2N_r^2}, &
\text{if } N_r^2<0,
\end{cases}
\end{equation}
where the range of 
$\tan ^{-1}$ is taken to be 
$[-{\rm \pi} /2,{\rm \pi} /2]$. We recall from (2.5) that thermal wind balance requires 
$N_r^2 = 2({V}/{r}+\varOmega )({\partial V}/{\partial z})$. The ranges of unstable 
$\theta$ are given by
$$\begin{gather} [\theta_1,\theta_2], \text{ if } N_r^2,\varPhi>0, \quad [0,\theta_1]\cup[\theta_2+{\rm \pi},{\rm \pi}], \text{ if } N_r^2>0,\varPhi<0, \end{gather}$$
$$\begin{gather}[\theta_2+{\rm \pi},\theta_1+{\rm \pi}], \text{ if } N_r^2<0,\varPhi>0, \quad [0,\theta_2]\cup[\theta_1+{\rm \pi},{\rm \pi}], \text{ if } N_r^2,\varPhi<0, \end{gather}$$where
\begin{equation} \theta_{1,2} = \tan^{{-}1}\left\{\frac{N_r^2}{N_z^2}\left[1 \mp \sqrt{1-{\frac{N_z^2\varPhi}{N_r^4}}}\right]\right\}, \end{equation}
where we have again taken the range of 
$\tan ^{-1}$ as 
$[-{\rm \pi} /2,{\rm \pi} /2]$.
The results in (3.10) and (3.11), in the limit of no curvature effects, i.e. 
$V/r=0$, are in agreement with the previously known two-dimensional inviscid symmetric instability criteria derived based on the normal mode approach (Ooyama Reference Ooyama1966; Hoskins Reference Hoskins1974). To the best of our knowledge, recovery of the symmetric instability criteria in the local stability framework has not been pointed out before. Furthermore, (3.10)–(3.14) provide corrections due to finite curvature effects, which are difficult to capture in the normal mode approach (Buckingham et al. Reference Buckingham, Gula and Carton2021a) although they are negligible only in the limit of 
$|V/r|\ll |\varOmega |$. Equation (3.10), which is in agreement with the symmetric instability criterion based on Lagrangian parcel arguments (Solberg Reference Solberg1936; Emanuel Reference Emanuel1979) and the recent studies of Buckingham et al. (Reference Buckingham, Gula and Carton2021a,Reference Buckingham, Gula and Cartonb), also establishes the local stability framework as an alternative approach to understand curvature effects. Finally, an insightful interpretation of (3.10) is that the square of the absolute circulation (
$2{\rm \pi} (\varOmega r^2 + Vr)$) decreases with radius when moving along iso-density surfaces (Emanuel Reference Emanuel1979); this form of the symmetric instability criterion has been previously used to interpret instabilities in pancake vortices (Negretti & Billant Reference Negretti and Billant2013; Yim & Billant Reference Yim and Billant2016; Yim, Billant & Ménesguen Reference Yim, Billant and Ménesguen2016).
In the limit of no radial stratification, i.e. 
$N_r=0$, the instability criterion in (3.10) reduces to 
$\varPhi < 0$, with the maximum growth rate being given by 
$\lambda ^* = \sqrt {-\varPhi }$ occurring at 
$\theta ^*=0$. This no-radial-stratification inviscid limit was highlighted using the local stability approach by Singh & Mathur (Reference Singh and Mathur2019), and is consistent with the large axial wavenumber limit of Billant & Gallaire (Reference Billant and Gallaire2005). In this unstable regime of 
$\varPhi <0$, the range of unstable 
$\theta$ reduces from 
$[0,{\rm \pi} /2]$ for 
$N_z=0$ to 
$[0,{\tan ^{-1}\sqrt {-\varPhi /N_z^2}}]$ for a finite 
$N_z$. Finally, in the limit of a homogeneous fluid, i.e. 
$N_r=N_z=0$, (3.9) gives the classical inviscid centrifugal instability criterion (Kloosterziel & Van Heijst Reference Kloosterziel and Van Heijst1991). This limit of centrifugal instability in unstratified vortices with background rotation has been considered in previous local stability studies too (Sipp & Jacquin Reference Sipp and Jacquin2000; Nagarathinam, Sameen & Mathur Reference Nagarathinam, Sameen and Mathur2015).
In terms of the non-dimensional base flow parameters in (3.4a,b), the inviscid instability criterion in (3.10) can be written as
\begin{equation} \frac{\tan\varTheta}{\tan\varGamma}<1, \end{equation}
with the most unstable 
$\theta$ given by
\begin{equation}
\theta^* = \begin{cases} \dfrac{\rm \pi}{4} -
\frac{1}{2}\tan^{{-}1}{\{\tan\varTheta/2 -
1/(2\tan\varGamma)\}}, & \text{if } N_r^2\ge0 \\
\dfrac{3{\rm \pi}}{4} -
\frac{1}{2}\tan^{{-}1}{\{\tan\varTheta/2 -
1/(2\tan\varGamma)\}}, & \text{if } N_r^2<0. \end{cases}
\end{equation}It is worth highlighting that
\begin{equation} \frac{\tan\varTheta}{\tan\varGamma} = {\frac{\varPhi N_z^2}{N_r^4}} = \frac{N_z^2}{(\partial V/\partial z)^2}\frac{(\partial V/\partial r +V/r + 2\varOmega)}{2(\varOmega + V/r)} \end{equation}
represents a modified gradient Richardson number. The inviscid instability regimes on the 
$(\tan \varGamma ,\tan \varTheta )$ plane are shown in figure 2. In the first and third quadrants, 
$\tan \varTheta = \tan \varGamma$ separates the stable and unstable regions, whereas the entire second and fourth quadrants (
$\varPhi <0 =$ 
$\hat {\varPhi }=-1$ in these quadrants) are unstable. It is worth recalling that all inviscid instabilities are monotonic in nature (zero imaginary part for the growth rate).
Figure 2. Inviscidly stable (white background) and unstable (grey background) regions on the 
$(\tan \varGamma ,\tan \varTheta )$ plane. The origin 
$O$ is at 
$(\tan \varGamma ,\tan \varTheta ) = (0,0)$.
3.3. Viscous regime
In the limit of 
$Sc=1$, the solutions of (3.2) are 
$\tilde {\lambda } = -\tilde {c}$, 
$\pm [-\hat {\varPhi }\{(1+\tan \varTheta \tan \varGamma ) + (1-\tan \varTheta \tan \varGamma )\cos 2\theta -2\tan \varGamma \sin 2\theta \}]^{1/2}-\tilde {c}$. The corresponding dimensional values of 
$\lambda$ are 
$\lambda = -c_\nu |\boldsymbol {k}|^2$, 
$\lambda = \lambda _i - c_\nu |\boldsymbol {k}|^2$, where 
$\lambda _i$ are the inviscid eigenvalues given by (3.9). In other words, for 
$Sc = 1$, all the eigenvalues from the inviscid limit get reduced by 
$c_\nu |\boldsymbol {k}|^2$, thus resulting in the inviscid growth rates being reduced by viscous effects, and the occurrence of two stable eigenvalues 
$(=-c_\nu |\boldsymbol {k}|^2)$. As a consequence, the eigenvalues for 
$Sc=1$ and 
$|\boldsymbol {k}|=0$ coincide with the inviscid eigenvalues, thus establishing an equivalence between inviscid instabilities and instabilities at 
$Sc=1$.
For 
$Sc\ne 1$, three non-trivial eigenvalues governed by (3.2) occur. While two of them could possibly be a continuous extension of the two non-trivial eigenvalues from the inviscid limit, the counterpart of the third eigenvalue was trivial and stable in the inviscid limit. For 
$Sc\ne 1$, we proceed to plot the non-dimensionalized growth rate (
$\tilde {\sigma }=\sigma /\sqrt {|\varPhi /2|}$) as a function of the perturbation parameters (3.5a,b) for given base flow parameters (3.4a,b). As highlighted earlier in this section, our choice of non-dimensionalization allows us to easily explore the important regime of weak but finite viscous effects.
Figure 3 shows the growth rate (in log scale) plotted as a function of 
$\theta$ and 
$\tilde {c}$ for 
$\tan \varGamma = 1$, and three different values of 
$\tan \varTheta$ 
$(1.2, 1$ and 
$0.4)$, which progressively takes the system from stable to neutrally stable to unstable in the inviscid limit (see first quadrant of figure 2). We recall from (3.15) that 
$\tan \varTheta /\tan \varGamma <1$ is the inviscid instability criterion, which also holds for 
$Sc = 1$. It is also worth noting that for 
$Sc = 1$, the instability is always monotonic in nature. In figure 3, panels (a,b,c) and (d,e,f) correspond to 
$Sc=0.1$ and 5, respectively.
Figure 3. Growth rate (
$\tilde {\sigma }$, on a log scale) as a function of 
$\theta$ and 
$\tilde {c}$ for 
$Sc = 0.1$ (a–c) and 
$Sc = 5$ (d–f), with 
$\tan \varTheta = 1.2, 1$ and 
$0.4$ in (a,d), (b,e) and (c, f), respectively. All the plots correspond to 
$\tan \varGamma = 1$, thus belonging to the first quadrant of figure 2. The green and brown regions correspond to zero (monotonic instability) and non-zero (oscillatory instability) imaginary part for the growth rate. In each plot, the unfilled circle corresponds to the location of maximum growth rate on the entire (
$\theta$, 
$\tilde {c}$) plane, with the red colour indicating a location away from 
$\tilde {c} = 0$ axis. No instabilities occur for 
$\theta >{\rm \pi} /2$ in all the plots in this figure. The values of 
$\tan \varTheta /\tan \varGamma$ for (a,b,d,e) correspond to the inviscidly stable regime, and hence the growth rate is zero on the 
$\tilde {c} = 0$ axis.
For 
$\tan \varTheta = 1.2$ (figure 3a), i.e. the inviscidly stable regime, we observe unstable growth rates at 
$Sc = 0.1$ though no instabilities occur for 
$Sc = 1$. The monotonic and oscillatory (non-zero imaginary part for the growth rate) instability regions are shown in green and brown, respectively. The monotonic instability dominates the oscillatory instability on the 
$(\theta ,\tilde {c})$ plane, with the overall maximum growth rate occurring at 
$(\theta ^*,\tilde {c}^*) = (1.05\ \text {rad},0.14)$ within the monotonic instability region. The monotonic instability is completely suppressed at a threshold value of 
$\tilde {c}^t = 0.54$. For the oscillatory instability, the maximum growth rate occurs at 
$(\theta ^*,\tilde {c}^*) = (0.49\ \text {rad},0.035)$, and it is completely suppressed beyond 
$\tilde {c}= 0.09$. In terms of the range of 
$\theta$, monotonic and oscillatory instabilities occur within [0.8,1.45] rad and [0.19,0.67] rad, respectively. The width of these unstable 
$\theta$ ranges are used later (§ 3.3.2) to characterize the instability region on the 
$(\theta ,\tilde {c})$ plane.
For 
$\tan \varTheta = 1$ (figure 3b), i.e. on the inviscid neutral stability boundary, monotonic and oscillatory instabilities are again observed at 
$Sc=0.1$, with an overall increase in the growth rates compared with 
$\tan \varTheta = 1.2$. The monotonic and oscillatory instability regions are also slightly larger for 
$\tan \varTheta = 1$ than for 
$\tan \varTheta = 1.2$. For 
$\tan \varTheta = 0.4$ and 
$Sc=0.1$ (figure 3c), i.e. the inviscidly unstable regime, the regions of instability and the associated growth rates are noticeably larger. Furthermore, the maximum growth rate on the entire 
$(\theta ,\tilde {c})$ plane now occurs at 
$(\theta ^*,\tilde {c}^*) = (0.93\ \text {rad}, 0)$, still within the monotonic instability region. For reference, the monotonic instability for 
$(\tan \varGamma ,\tan \varTheta ) = (1,0.4)$ at 
$Sc = 1$ is also strongest at the same 
$(\theta ^*,\tilde {c}^*)$ location (see (3.16)).
Figure 3(d–f) corresponds to the same parameters as panels (a,b,c), except for 
$Sc$, which is now set at 
$5$. In the inviscidly stable regime of 
$\tan \varTheta = 1.2$ (figure 3d), monotonic instability is observed while oscillatory instability is absent. The overall growth rates are smaller in magnitude than for 
$Sc = 0.1$, and the maximum growth rate occurs at 
$(\theta ^*,\tilde {c}^*) = (0.5\ \text {rad},0.35)$. With an increase in 
$\tan \varTheta$ to unity (figure 3e), the growth rate magnitudes increase, and oscillatory instability appears. In contrast to 
$Sc = 0.1$, oscillatory instability occurs at larger 
$\theta$ than that for the monotonic instability at 
$Sc = 5$. Finally, in the inviscidly unstable regime of 
$\tan \varTheta = 0.4$ (figure 3f), oscillatory instability again appears alongside the monotonic instability, the region of instability on the (
$\theta ,\tilde {c}$) plane is substantially larger than for smaller 
$\tan \varTheta$, and the maximum growth rate occurs at the same 
$(\theta ^*,\tilde {c}^*) = (0.93\ \text {rad},0)$ as that for the monotonic instability at 
$Sc = 1$.
In summary of figure 3, for 
$Sc\ne 1$, monotonic instability appears in the inviscidly stable regime too, while a new oscillatory instability appears in inviscidly stable as well as unstable regimes. The maximum growth rate on the 
$(\theta ,\tilde {c})$ plane, however, always lies in the monotonic instability region, irrespective of 
$Sc$, 
$\tan \varGamma$ and 
$\tan \varTheta$. We show later in §§ 3.3.3 and 3.4 that the monotonic instability observed for 
$Sc\ne 1$ is a continuous extension of the monotonic symmetric instability that occurs at both 
$Sc = 1$ and the inviscid limit. We proceed to investigate (in figure 4) the instabilities for 
$Sc\ne 1$ in the fourth quadrant of figure 2, the entirety of which is inviscidly unstable due to 
$\varPhi <0$.
Figure 4. Growth rate (
$\tilde {\sigma }$, on a log scale) as a function of 
$\theta$ and 
$\tilde {c}$ for (a) 
$Sc = 0.01$, (b) 
$Sc = 1$, (c) 
$Sc = 5$, with 
$(\tan \varGamma ,\tan \varTheta ) = (1,-1)$, and thus belonging to the fourth quadrant in figure 2. The green and brown regions correspond to zero (monotonic instability) and non-zero (oscillatory instability) imaginary part for the growth rate. In each plot, the unfilled circle corresponds to the location of maximum growth rate on the entire (
$\theta$, 
$\tilde {c}$) plane. In (c), the growth rate corresponding to oscillatory instability vanishes on the 
$\tilde {c} = 0$ axis.
Figure 4(a,b) shows that only monotonic instability occurs at 
$Sc = 0.1$ and 
$Sc = 1$ for 
$(\tan \varGamma ,\tan \varTheta ) = (1,-1)$. The range of unstable 
$\theta$ is noticeably larger for 
$Sc = 0.1$ compared with 
$Sc = 1$. At 
$Sc = 5$ (figure 4c), oscillatory instability appears along with the monotonic instability. Finally, we note that the most unstable 
$(\theta ,\tilde {c})$ is 
$(2.75 \ \text {rad},0)$ , which lies within the monotonic instability region, for all three 
$Sc$ in figure 4. Furthermore, we find that the largest growth rate is always associated with the monotonic instability for every unstable (
$\tan \varGamma ,\tan \varTheta ,Sc$). Towards obtaining neutral stability boundaries in terms of the base flow parameters and 
$Sc$, we proceed to plot the most unstable growth rate in the monotonic and oscillatory instability regions of plots such as in figures 3 and 4, on the (
$\tan \varGamma ,\tan \varTheta$) plane.
3.3.1. Dominant instability characteristics
The maximum growth rate on the (
$\theta ,\tilde {c}$) plane, and the corresponding location at which it occurs, are plotted as a function of 
$\tan \varGamma$ and 
$\tan \theta$ for 
$Sc=0.1$ and 
$5$ in the panels (a,b,c) and (g,h,i) of figure 5, respectively. We recall that the maximum growth rate always corresponds to the monotonic instability. For reference, the equivalent plots for 
$Sc = 1$ are shown in panels (d,e,f), with the black lines in (a,b,d,g,h,i) representing the boundary separating stable and unstable regions for 
$Sc = 1$. The panels (a,d,g), (b,e,h) and (c,f,i) show the distributions of maximum growth rate 
$\tilde {\sigma }^{max}$ (on a log scale), and the most unstable (
$\theta ,\tilde {c}) = (\theta ^*,\tilde {c}^*$), respectively.
Figure 5. Maximum growth rate (
$\log \tilde {\sigma }^{max}$) on the (
$\theta$, 
$\tilde {c}$) plane, plotted as a function of 
$\tan \varGamma$ and 
$\tan \varTheta$ for (a) 
$Sc = 0.1$, (d) 
$Sc = 1$ and (g) 
$Sc = 5$. Corresponding locations on the (
$\theta$, 
$\tilde {c}$) plane where the maximum growth rate occurs are plotted as 
$\theta ^*$ (b,e,h) and 
$\tilde {c}^*$ (c, f,i). The black line in all the plots represents the neutral stability (
$\tilde {\sigma }^{max} = 0$) boundary for 
$Sc = 1$.
For 
$Sc = 0.1$, monotonic instability is observed in the inviscidly stable region of the first quadrant of the (
$\tan \varGamma ,\tan \varTheta$) plane, with the neutral stability boundary being steeper than that for 
$Sc = 1$ (figure 5a). In contrast to 
$Sc=1$, some part of the monotonic instability region for 
$Sc = 0.1$ has the maximum growth rate occurring at finite viscous effects, i.e. 
$\tilde {c}^*>0$ (figure 5c). Interestingly, 
$\tilde {c}^*>0$ occurs for the whole of the unstable region that was inviscidly stable, and some part of the inviscidly unstable region as well. The 
$\tilde {c}^* = 0$ boundary is coincident with a distinct 
$\theta ^*$ contour (figure 5b), on which the value of 
$\theta ^*$ is equal to the inviscid value. The value of 
$\theta ^*$ seems to abruptly increase as we enter the 
$\tilde {c}^*>0$ region. In the third quadrant of the (
$\tan \varGamma ,\tan \varTheta$) plane, we observe similar characteristics as in the first quadrant, but with the 
$\theta ^*$ variations occurring in the neighbourhood of 
$\theta ^* = {\rm \pi}$ instead of 
$\theta ^* = 0$. Finally, the characteristics of dominant instability in the second and fourth quadrants of the (
$\tan \varGamma ,\tan \varTheta$) plane for 
$Sc = 0.1$ are similar to those for 
$Sc = 1$.
For 
$Sc = 5$, some part of the inviscidly stable regions in the first and third quadrants of the (
$\tan \varGamma ,\tan \varTheta$) plane becomes unstable (figure 5g), as observed for 
$Sc = 0.1$. In contrast to 
$Sc = 0.1$, however, 
$\tilde {c}^*>0$ seems to occur predominantly in the inviscidly stable region for 
$Sc = 5$. Also, in the first quadrant of the (
$\tan \theta ,\tan \varGamma$) plane for 
$Sc = 5$, 
$\theta ^*$ decreases as we move into the 
$\tilde {c}^*>0$ region. As observed for 
$Sc = 0.1$, the dominant instability characteristics in the second and fourth quadrants seem similar to 
$Sc = 1$ for 
$Sc = 5$. In summary of figure 5, some parts of inviscidly stable regions become monotonically unstable for 
$Sc\ne 1$, with the extent of instability depending on the value of 
$Sc$. 
$Sc\ne 1$ also introduces the occurrence of the maximum growth rate at finite viscous effects, i.e. 
$\tilde {c}^*>0$, with such regions showing 
$\theta ^*$ distributions that are different from the inviscid limit, i.e. 
$Sc = 1$.
Figure 6 shows the dominant instability characteristics from within the oscillatory instability regions (brown regions in figures 3 and 4), plotted on the (
$\tan \varGamma ,\tan \varTheta$) plane for 
$Sc = 0.1$ (a–c) and 
$5$ (d–f). We recall that oscillatory instability (non-zero imaginary part for the growth rate) does not occur for 
$Sc = 1$, and monotonic instability dominates oscillatory instability for each (
$\tan \varGamma ,\tan \varTheta ,Sc$). Figure 6(a) shows that the entire inviscidly unstable region in the first quadrant is susceptible to oscillatory instability for 
$Sc = 0.1$, with the corresponding 
$\theta ^*_o$ being around 
${\rm \pi} /2$ close to the origin, and approaching zero for large (
$\tan \varGamma ,\tan \varTheta$) (figure 6b). Figure 6(c) shows that oscillatory instability occurs at finite viscous effects (
$\tilde {c}^*_o>0$), with rapid variations observed within the inviscidly unstable region. Some part of the inviscidly stable region is also seen to be susceptible to oscillatory instability, with a growth rate that is much smaller than that for monotonic instability. The extent of oscillatory instability is also smaller than that of monotonic instability, as indicated by the neutral stability boundary (dashed line) in figure 6(a). The oscillatory instability features in the third quadrant are similar to those in the first quadrant, with the 
$\theta _0^*$ variations occurring around 
${\rm \pi}$ rather than 
$0$. Finally, the inviscidly unstable second and fourth quadrants are not affected by oscillatory instability for 
$Sc = 0.1$.
Figure 6. Maximum growth rate (
$\log \tilde {\sigma }^{max}_o$) in the oscillatory instability region on the (
$\theta$, 
$\tilde {c}$) plane, plotted as a function of 
$\tan \varGamma$ and 
$\tan \varTheta$ for (a) 
$Sc = 0.1$ and (d) 
$Sc = 5$. Corresponding locations on the (
$\theta$, 
$\tilde {c}$) plane where the maximum oscillatory instability growth rate occurs are plotted as 
$\theta ^*_o$ (b,e) and 
$\tilde {c}^*_o$ (c, f). The solid black line in all the plots represents the 
$\tilde {\sigma }^{max} = 0$ boundary for 
$Sc = 1$. The black dashed line in all the plots indicates the neutral stability boundary (
$\tilde {\sigma }^{max} = 0$, see figure 6) for the respective Sc.
For 
$Sc = 5$, the inviscidly unstable region in the first quadrant is again entirely susceptible to the oscillatory instability (figure 6d), with 
$\tilde {\sigma }^{max}_o$, 
$\theta _o^*$ and 
$\tilde {c}^*_o$ all being relatively large close to the 
$\tan \varTheta = 0$ axis. In addition, there is also a small portion of the inviscidly stable region where oscillatory instability is observed, with the corresponding 
$\theta _o^*$ and 
$\tilde {c}^*_o$ being close to zero. Finally, in contrast to 
$Sc = 0.1$, the inviscidly unstable second and fourth quadrants are entirely susceptible to the oscillatory instability for 
$Sc = 5$. The abrupt changes in 
$\theta _o^*$ close to the 
$\tan \varTheta = 0$ axis in the second and fourth quadrants of figure 6(e) are attributed to the switch of the location of maximum oscillatory instability growth rate from one unstable region to another on the (
$\theta$, 
$\tilde {c}$) plane (see brown regions in figure 4(c) for an example of two disconnected oscillatory instability regions on the (
$\theta$, 
$\tilde {c}$) plane). For each 
$(\tan \varGamma, \tan \varTheta )$, we proceed to identify the range of unstable 
$\theta$ on the (
$\theta ,\tilde {c}$) plane for each of the monotonic and oscillatory instabilities. The width of such an unstable range, indicative of the range of perturbation wave vectors that are unstable, is then plotted as a function of 
$\tan \varGamma$ and 
$\tan \varTheta$ for various 
$Sc$, as discussed in § 3.3.2.
3.3.2. Range of unstable perturbations
Figure 7 shows the width of unstable 
$\theta$ range on the 
$(\tan \varGamma ,\tan \varTheta )$ plane for the monotonic (a–c) and oscillatory (d–f) instabilities. The panels (a,d), (b,e) and (c,f) correspond to 
$Sc = 0.1$, 
$1$ and 
$5$, respectively. For 
$Sc = 1$ (figure 7b), monotonic instability occurs in a small range around 
$\theta ^*$ (previously shown in figure 5e) in the region close to the neutral stability boundary (black line) in the first and third quadrants. The unstable 
$\theta$ width increases as we move away from the neutral stability boundary towards the 
$\tan \varTheta = 0$ axis. In the second and fourth quadrants, a significant range of 
$\theta$ is monotonically unstable near the 
$\tan \varTheta = 0$ and 
$\tan \varGamma = 0$ axes. In regard to oscillatory instability, as expected, the entire plane of 
$(\tan \varGamma ,\tan \varTheta )$ is stable for 
$Sc = 1$ (figure 7e).
Figure 7. Width of the unstable 
$\theta$ range in the monotonic instability (a–c) and oscillatory instability (d–f) regions on the (
$\theta$, 
$\tilde {c}$) plane, plotted as a function of 
$\tan \varGamma$ and 
$\tan \varTheta$ for (a,d) 
$Sc = 0.1$, (b,e) 
$Sc=1$ and (c, f) 
$Sc = 5$. The solid black line in all the plots represents the neutral stability boundary (
$\tilde {\sigma }^{max} = 0$) for 
$Sc = 1$.
In comparison with 
$Sc = 1$, the width of unstable 
$\theta$ range for monotonic instability is noticeably larger throughout the (
$\tan \varGamma ,\tan \varTheta$) plane for 
$Sc = 0.1$ (figure 7a). A similar conclusion holds for 
$Sc = 5$ too (figure 7c), with the extent of increase in the width of the unstable 
$\theta$ range being relatively smaller compared with 
$Sc = 0.1$. For oscillatory instability, the width of unstable 
$\theta$ is, in general, smaller than that for monotonic instability for both 
$Sc = 0.1$ (figure 7d) and 
$5$ (figure 7f). For 
$Sc = 0.1$, the width of the unstable 
$\theta$ range decreases as we move away from the origin (figure 7d). For 
$Sc = 5$, oscillatory instability appears in all the four quadrants of the (
$\tan \varGamma ,\tan \varTheta$) plane (figure 7f), with relatively large width for the unstable 
$\theta$ range occurring in a region slightly shifted away from the 
$\tan \varTheta = 0$ axis. The regions close to 
$\tan \varGamma = 0$ in the second and fourth quadrants correspond to relatively large width for the unstable 
$\theta$ range of the oscillatory instability, which is similar to what we observed for the monotonic instability in figure 7(c).
The neutral stability boundaries in figure 5 (
$\tilde {\sigma }^{max} = 0$) and the 
$\tilde {\sigma }^{max}_o = 0$ boundaries for the oscillatory instability in figure 6 are all observed to be straight lines passing through the origin on the (
$\tan \varGamma ,\tan \varTheta$) plane, indicating their dependence only on the ratio 
$\tan \varTheta /\tan \varGamma$. As highlighted in § 3.3, 
$\tan \varTheta /\tan \varGamma$ represents a modified gradient Richardson number. We proceed to investigate the variation (with 
$Sc$) of the critical value of 
$\tan \varTheta /\tan \varGamma$ that separates stable and unstable regions.
3.3.3. Neutral stability boundaries
Figure 8(a) shows the regions of monotonic and oscillatory instabilities on the plane of 
$Sc$ on the 
$x$-axis and 
$\tan \varTheta /\tan \varGamma$ on the 
$y$-axis, with the corresponding neutral stability boundaries shown in blue and red, respectively. The neutral stability boundary for the monotonic instability (blue dashed line) also represents the overall neutral stability boundary, thus rendering the region above the blue dashed line stable. The threshold 
$\tan \varTheta /\tan \varGamma$ above which the flow is stable increases as 
$Sc$ goes away from unity. At 
$Sc = 0.1$ or 
$10$, the threshold 
$\tan \varTheta /\tan \varGamma$ is around 
$3.1$ times larger than the value at 
$Sc = 1$. The neutral stability boundary for the oscillatory instability (red dashed line) lies below the overall neutral stability boundary at all 
$Sc\ne 1$. In contrast to 
$Sc<1$. the red dashed line for 
$Sc>1$ always remains close to (but above) 
$\tan \varTheta /\tan \varGamma = 1$, reaching a value of 
$9/8$ (see Appendix D) as 
$Sc\to \infty$. The red solid line captures the lower threshold of 
$\tan \varTheta /\tan \varGamma$ below which oscillatory instability does not occur. For all 
$Sc<1$, the red solid boundary occurs below the red dashed boundary, and there exists no red solid boundary for 
$Sc>1$, indicating that oscillatory instability occurs for all 
$Sc\ne 1$. For 
$Sc<1$, oscillatory instability occurs between 
$\tan \varTheta /\tan \varGamma =0$ and the red dashed line, whereas the entire region below the red dashed line is susceptible to oscillatory instability for 
$Sc>1$. In other words, a finite (
$Sc<1$) or an infinite (
$Sc>1$) range of 
$\tan \varTheta /\tan \varGamma$ becomes susceptible to oscillatory instability as soon as 
$Sc$ becomes different from unity. In the centrifugally unstable regime in the limit of no radial stratification (
$\varPhi <0$, 
$N_r^2 = 0$), oscillatory instability was previously reported only for 
$Sc>1$ (Singh & Mathur Reference Singh and Mathur2019), which is consistent with the 
$\tan \varTheta /\tan \varGamma \to -\infty$ limit in figure 8(a).
Figure 8. (a) Neutral stability boundaries on the (
$Sc$, (
$\tan \varTheta /\tan \varGamma$)) plane for the overall instability (blue dashed line) and the oscillatory instability (red dashed and solid lines). The blue background indicates region with monotonic instability only, and magenta background indicates regions with monotonic and oscillatory instabilities. (b) The wave vector orientation 
$\theta ^*$ at which the maximum growth rate occurs on the neutral stability boundaries shown in (a). The cyan dashed and the cyan solid lines indicate the 
$\theta ^*$ corresponding to monotonic instability along the neutral stability boundaries (red dashed and red solid lines, respectively) for the oscillatory instability.
In figure 8(b), the blue line shows the variation of 
$\theta ^*$, the location of maximum growth rate, along the overall neutral stability boundary, which is the same as that for the monotonic instability. For 
$Sc\ll 1$, it is around 
$1.1$ rad, and it decreases towards the inviscid value (3.11) of 
${\rm \pi} /4$ at 
$Sc=1$, before decreasing further towards 
$0.03$ rad for 
$Sc\gg 1$. In contrast, along the neutral stability boundary (red dashed line in figure 8a) for the oscillatory instability, 
$\theta _o^*$ increases from around 
$0.1$ rad for 
$Sc\ll 1$ towards 
$0.9$ rad for 
$Sc\gg 1$. On the lower neutral stability boundary (red solid line in figure 8a) for the oscillatory instability, 
$\theta _o^*$ is around 
$0.68$ rad for 
$Sc<1$. The variation of 
$\theta ^*$ corresponding to monotonic instability (cyan dashed and cyan solid lines in figure 8b) along the neutral stability boundaries for the oscillatory instability (red dashed and red solid lines in figure 8a) indicate that the monotonic and oscillatory instabilities are separated in 
$\theta$ when they occur simultaneously. While not shown in figure 8, the separation in 
$\theta$ between monotonic and oscillatory instabilities has been verified in all the regions where they occur simultaneously. Finally, with respect to 
$\tilde {c}^*$ and 
$\tilde {c}^*_o$, they are found to be zero along the neutral stability boundaries for the monotonic (blue dashed line in figure 8a) and oscillatory (red dashed line in figure 8a) instabilities, respectively. Along the lower neutral stability boundary for the oscillatory instability (red solid line in figure 8a), however, 
$\tilde {c}^*_o$ is finite and varies significantly.
As shown in Appendix D, it is possible to derive closed form expressions for the neutral stability boundaries on the (
$Sc$, 
$(\tan \varTheta /\tan \varGamma )$) plane. The resulting unstable regions are
\begin{equation} \frac{\tan\varTheta}{\tan\varGamma} < \frac{(1+Sc)^2}{4Sc} \end{equation}for the monotonic instability, and
\begin{equation}
\left.\begin{array}{rl} 0 <
\dfrac{\tan\varTheta}{\tan\varGamma} <
\dfrac{(3Sc+1)^2}{8Sc(Sc+1)}, & \text{if } Sc<1 \\
\dfrac{\tan\varTheta}{\tan\varGamma} <
\dfrac{(3Sc+1)^2}{8Sc(Sc+1)}, & \text{if } Sc>1
\end{array}\right\}
\end{equation}
for the oscillatory instability. In the centrifugally stable regime 
$(\tan \varTheta /\tan \varGamma >0)$ with zero curvature effects, (3.18)–(3.19) coincide with the neutral stability boundaries of McIntyre (Reference McIntyre1970).
3.4. Reduction of parameters
We have thus far presented our results in terms of five different non-dimensional parameters: two of them describing the base flow (3.4a,b), two more describing the perturbations (3.5a,b) and 
$Sc$. It is, however, possible to choose an alternate non-dimensionalization of (3.1) to further reduce the number of governing non- dimensional parameters, albeit at the cost of not separating the base flow and perturbation characteristics. Specifically, (3.2) can be re-written as
\begin{equation} {\lambda}^{*3} + \left(\frac{1-Sc}{Sc}\right){\lambda}^{*2} + [\alpha_1+\alpha_2-\beta_1-\beta_2]{\lambda}^* + \left(\frac{1-Sc}{Sc}\right)[\alpha_1-\beta_1]=0, \end{equation}
where 
${\lambda }^* = (\lambda +c_\nu |\boldsymbol {k}|^2)/(c_\nu |\boldsymbol {k}|^2$), and
\begin{equation} \left.\begin{gathered} \alpha_1 = \frac{\varPhi\cos^2\theta}{(c_\nu|\boldsymbol{k}|^2)^2}, \quad \alpha_2 = \frac{N_z^2\sin^2\theta}{(c_\nu|\boldsymbol{k}|^2)^2}, \\ \beta_1 = \frac{\partial{V}}{\partial z}\frac{(\varOmega + {V}/{r})\sin2\theta}{(c_\nu|\boldsymbol{k}|^2)^2}, \quad \beta_2 = \frac{N_r^2\sin2\theta}{2(c_\nu|\boldsymbol{k}|^2)^2}. \end{gathered}\right\} \end{equation}
In terms of 
$\lambda ^*$, instability occurs if 
$\textrm {Re}(\lambda ^*)>1$. We also define a non-dimensional growth rate 
$\bar {\sigma } = \sigma /(c_\nu |\boldsymbol {k}|^2)$ 
$= \max (\textrm {Re}(\lambda ^*) - 1)$, where 
$\sigma$ is the maximum real part of the eigenvalues 
$\lambda$. Unlike the non-dimensionalization used in (3.2), the inviscid limit of 
$c_\nu = 0$ represents a singularity in (3.20), and hence the weakly viscous limit of 
$\nu \to 0$ is not easily recovered. Additionally, the thermal wind relation in (2.5) requires 
$\beta _1 = \beta _2 = \beta$.
Equation (3.20) shows that instability characteristics are dependent only on three different parameters, namely 
$\alpha _1-\beta$, 
$\alpha _2-\beta$ and 
$Sc$. As shown in Appendix B, the neutral stability boundaries for (3.20) are
$$\begin{gather} \frac{\alpha_1-\beta}{Sc} + \alpha_2-\beta+\frac{1}{Sc} = 0,\quad \text{if } s(\alpha_1-\beta) < s\frac{1+Sc}{1-Sc}, \end{gather}$$
$$\begin{gather}2(\alpha_1-\beta) + (\alpha_2-\beta)\left(1+\frac{1}{Sc}\right)+2 \left(1+\frac{2}{Sc}+\frac{1}{Sc^2}\right) = 0, \quad \text{if}\ s(\alpha_1-\beta) > s\frac{1+Sc}{1-Sc}, \end{gather}$$
where 
$s = {\textrm {sgn} (1-Sc)}$. In the limit of 
$Sc=1$, the neutral stability boundaries are governed by (3.22) only.
In the unstable neighbourhood 
$(({\alpha _1-\beta })/{Sc} + \alpha _2-\beta +{1}/{Sc} \lesssim 0)$ of the neutral stability boundary in (3.22), the unstable eigenvalue is real, and represents a continuous extension of the inviscid symmetric instability that occurs for 
$\alpha _1 + \alpha _2 - 2\beta + 1 < 0$. The other two eigenvalues (stable) are either (i) real, if 
$s(\alpha _1-\beta ) > s/(4Sc(Sc-1))$ or (ii) complex conjugates, if 
$s(\alpha _1-\beta ) < s/(4Sc(Sc-1))$. In the unstable neighbourhood 
$(2(\alpha _1-\beta ) + (\alpha _2-\beta )(1+{1}/{Sc})+2(1+{2}/{Sc}+{1}/{Sc^2}) \lesssim 0)$ of the neutral stability boundary in (3.23), the unstable eigenvalues are complex conjugates, and do not represent a continuous extension of the inviscid or 
$Sc=1$ instability; the corresponding instability is an oscillatory instability, based on the non-zero imaginary parts for the unstable eigenvalues. The third eigenvalue is real and stable. Finally, as shown in Appendix C, monotonic and oscillatory instabilities cannot simultaneously occur for given values of 
$\alpha _1-\beta$, 
$\alpha _2-\beta$ and 
$Sc$. As a consequence, even for given values of 
$(\tan \varGamma ,\tan \varTheta ,\theta ,\tilde {c},Sc)$, since 
$\alpha _1$, 
$\alpha _2$ and 
$\beta$ are uniquely specified (
$N_z^2>0$ assumed), monotonic and oscillatory instabilities cannot occur simultaneously.
Figure 9 shows the neutral stability boundaries on the (
$(\alpha _1-\beta )$, 
$(\alpha _2-\beta )$) plane for 
$Sc = 0.1$ (figure 9a) and 
$Sc = 5$ (figure 9b). For reference, the neutral stability boundary for 
$Sc = 1$ (or inviscid instability) is shown in black in figure 9(a,b). For a given 
$Sc$, the neutral stability boundary switches from the monotonic instability (magenta boundary) to the oscillatory instability (red boundary) at 
$\alpha _1-\beta = (1+Sc)/(1-Sc)$. The normalized growth rate associated with the monotonic instability (
$\bar {\sigma }_m$, calculated numerically from (3.20)) increases somewhat rapidly as we move way from its corresponding neutral stability boundary. In contrast, the normalized growth rate associated with the oscillatory instability (
$\bar {\sigma }_o$, calculated numerically from (3.20)) remains relatively small even as we move away from the corresponding neutral stability boundary. Inviscidly stable regions (to the right of the black boundary) (1s) can succumb to the monotonic instability if 
$s(\alpha _1-\beta )<-s$, or (2s) will remain stable if 
$-s<(\alpha _1-\beta )< s(Sc^2+3Sc+1)/(Sc(1-Sc))$, or (3s) can succumb to the oscillatory instability if 
$s(\alpha _1-\beta )>s(Sc^2+3Sc+1)/(Sc(1-Sc))$. In the inviscidly unstable regions, we observe qualitatively three different behaviours: (1u) monotonic instability only, (2u) both monotonic and oscillatory instabilities being absent and (3u) oscillatory instability only. The stable triangular region (region 2u) is bounded by 
$[\alpha _1-\beta , \alpha _2-\beta ] = [-1,0],[(Sc^2+3Sc+1)/(Sc(1-Sc)),(4Sc+1)/(Sc(Sc-1))]$, 
$[(1+Sc)/(1-Sc),2/(Sc(Sc-1))]$. The equation describing the boundary between regions (1u) and (3u) is derived in Appendix C (see (C5)).
Figure 9. Neutral stability boundaries (3.22)–(3.23) on the (
$(\alpha _1-\beta )$, 
$(\alpha _2-\beta )$) plane for (a) 
$Sc = 0.1$, (b) 
$Sc = 5$. The magenta solid lines separate the stable region from the monotonic instability region (growth rate 
$\bar {\sigma }_m$ shown using the green colour map). The red dashed lines separate the stable region from the oscillatory instability region (growth rate 
$\bar {\sigma }_o$ shown using the brown colour map). In both the plots, the magenta circle denotes the location where the neutral stability boundary switches from monotonic to oscillatory instability, and the black line denotes the neutral stability boundary for 
$Sc = 1$. In (a), labels (1s, 2s, 3s) correspond to inviscidly stable regions and (1u, 2u, 3u) correspond to inviscidly unstable regions.
In summary of this section, we identified new parameters 
$(\alpha _1-\beta , \alpha _2-\beta , Sc)$ based on an alternate non-dimensionalization, which allowed us to derive analytical expressions for various instability characteristics. Specifically, expressions for the neutral stability boundaries and the interface between monotonic and oscillatory instabilities on the (
$(\alpha _1-\beta )$, 
$(\alpha _2-\beta )$) plane were derived; in addition, it was also shown that monotonic and oscillatory instabilities cannot simultaneously occur for the same values of the parameters. It is, however, noteworthy that the parameters 
$\alpha _1$, 
$\alpha _2$ and 
$\beta$ are a combination of base flow and perturbation characteristics. The instability characteristics in terms of physically more relevant parameters, which also separate the effects of base flow and perturbations, were earlier presented in §§ 3.2–3.3.
4. Discussion
In this section, we address two different goals, namely (i) to highlight the effects of curvature, an aspect that is relatively less understood, and (ii) to explore typical oceanic regimes and comment on the possibility of the instabilities discussed in § 3.
4.1. Effects of curvature
As shown in § 3.3.3, the neutral stability boundaries for the monotonic and oscillatory instabilities are given by (3.18) and (3.19), respectively. In the absence of curvature effects, the neutral stability boundaries are obtained by replacing 
$\varGamma$ and 
$\varTheta$ by 
$\varGamma _M$ and 
$\varTheta _M$, where 
$\varGamma _M$ and 
$\varTheta _M$ are as defined in (3.6a,b). We recall that the subscript 
$M$ refers to the variables defined in the study of McIntyre (Reference McIntyre1970), which did not account for curvature effects (see § 3.1); (
$\varGamma _M,\varTheta _M$) are related to 
$(\varGamma ,\varTheta )$ by
\begin{equation} \tan{\varGamma}_M = \tan\varGamma (1+R_c),\quad \tan{\varTheta}_M = \tan\varTheta, \end{equation}
where 
$R_c = ({V}/r)/(\partial {V}/\partial r + 2\varOmega )$ is a non-dimensional parameter that quantifies the curvature effects. The unstable region (on the (
$Sc,(\tan \varTheta _M/\tan \varGamma _M)$) plane) corresponding to monotonic instability (3.18) can now be written in terms of 
$\tan \varGamma _M,\tan \varTheta _M$ and 
$R_c$ as
\begin{equation} \frac{\tan\varTheta_M(1+R_c)}{\tan\varGamma_M} < \frac{(1+Sc)^2}{4Sc}. \end{equation}Correspondingly, the unstable region for oscillatory instability (3.19) can be written as
\begin{equation}
\left.\begin{array}{rl} 0 <
\dfrac{\tan\varTheta_M(1+R_c)}{\tan\varGamma_M} <
\dfrac{(3Sc+1)^2}{8Sc(Sc+1)}, & \text{if }
Sc<1 \\
\dfrac{\tan\varTheta_M(1+R_c)}{\tan\varGamma_M} <
\dfrac{(3Sc+1)^2}{8Sc(Sc+1)}, & \text{if }
Sc>1.
\end{array}\right\}
\end{equation} Setting 
$R_c = 0$ in (4.2)–(4.3) recovers the limit of zero curvature exactly. For 
$R_c\to \infty$, the unstable regions asymptotically approach 
$\tan \varTheta _M/\tan \varGamma _M<0$ for monotonic instability at all 
$Sc$ and oscillatory instability for 
$Sc>1$; oscillatory instability does not occur for 
$Sc<1$ at 
$R_c = \infty$. For 
$R_c\to -\infty$, the unstable regions asymptotically approach 
$\tan \varTheta _M/\tan \varGamma _M>0$ for monotonic instability at all 
$Sc$ and oscillatory instability for 
$Sc>1$; oscillatory instability again does not occur for 
$Sc<1$ at 
$R_c = -\infty$. In other words, for strongly negative 
$R_c$, ignoring the curvature effects can result in qualitatively and quantitatively misleading conclusions in terms of which regions on the (
$Sc,(\tan \varGamma _M/\tan \varTheta _M)$) plane are unstable. Such misleading conclusions are a result of a change in the sign of 
$(1+R_c)$ at 
$R_c = -1$. To quantitatively visualize the effects of curvature, we proceed to plot the neutral stability boundaries on the (
$Sc,(\tan \varTheta _M/\tan \varGamma _M)$) plane for various values of 
$R_c$.
Figure 10 shows the neutral stability boundaries on the plane of 
$Sc$ and 
$\tan \varTheta _M/\tan \varGamma _M$ for different values of 
$R_c$. While the region below the boundaries is unstable for 
$R_c>-1$, the region above the corresponding neutral stability boundary is unstable for 
$R_c<-1$. In addition, for all 
$R_c$, there is a neutral stability boundary for oscillatory instability at 
$\tan \varTheta _M/\tan \varGamma _M = 0$ for 
$Sc<1$ in figure 10(b). In the absence of curvature effects, i.e. 
$R_c = 0$, the black neutral stability boundaries in figure 10(a,b) are exact. For 
$R_c\ne 0$, the neutral stability boundaries move away from the black boundaries; however, a theory that ignores curvature effects (McIntyre (Reference McIntyre1970), for example) always predicts the black boundaries as the neutral stability boundaries. For 
$R_c>0$ (red boundaries), the regions of both monotonic and oscillatory instabilities are smaller than what a theory with no curvature effects would predict. On the other hand, for 
$-1< R_c<0$ (blue boundaries), the actual regions of instability are larger than at 
$R_c = 0$. For a given 
$R_c>-1$, the distance between the exact and the curvature-effect-free neutral stability boundaries increases as 
$Sc$ moves away from unity. For 
$R_c<-1$, the regions of instability are above the corresponding neutral stability boundaries (shown in green for 
$R_c = -2$), which now lie in the negative half of the plane. Therefore, as mentioned earlier, a curvature-free theory would be completely misleading if 
$R_c<-1$.
Figure 10. Neutral stability boundaries on the (
$Sc,(\tan \varTheta _M/\tan \varGamma _M)$) plane (see (3.6a,b) for the definitions of 
$\tan \varTheta _M$ and 
$\tan \varGamma _M$) for (a) monotonic instability (4.2), and (b) oscillatory instability (4.3). The green, blue, black and red colours correspond to 
$R_c = -2$, 
$-0.5$, 0, 1, respectively, where 
$R_c$ is as defined in (4.1a,b).
4.2. Relevant ocean regimes
In this subsection, we discuss representative ocean regimes to explore the likelihood of occurrence of the instabilities described in this study. As shown in figure 8(a), the values of 
$Sc$ and 
$\tan \varTheta /\tan \varGamma$ determine if an eddy is susceptible to any short-wavelength instability, be it inviscid or due to differential diffusion between momentum and mass. With respect to 
$Sc$, molecular diffusion gives 
$Sc\approx 7$ for temperature in water, and 
$Sc\approx 700$ for salinity in water. In highly turbulent flows, the effective diffusion rates for momentum and scalar fields are often assumed to be nearly equal, resulting in a turbulent Schmidt number value near unity (Kays & Crawford Reference Kays and Crawford1993; Kays Reference Kays1994). The turbulent Schmidt number being close to unity, however, is not necessarily valid in stably stratified turbulent flows, as highlighted by numerical simulations (Venayagamoorthy & Stretch Reference Venayagamoorthy and Stretch2010) and field data (Salehipour et al. Reference Salehipour, Peltier, Whalen and MacKinnon2016). Specifically, 
$Sc$ can vary from around 
$0.7$ at small Richardson numbers, to around 
$20$ at large Richardson numbers (Elliott & Venayagamoorthy Reference Elliott and Venayagamoorthy2011). In addition, the buoyancy Reynolds number also plays a role in determining the value of turbulent 
$Sc$ (Salehipour & Peltier Reference Salehipour and Peltier2015).
Based on the definition in (3.3), 
$\varTheta$ can be physically interpreted as the orientation of isopycnals with respect to the vertical direction. Therefore, its value can range from being very small at front locations (steep isopycnals) to being very large far from fronts (horizontal isopycnals). As a consequence, the magnitude of 
$\tan \varTheta /\tan \varGamma$ is close to zero at sharp fronts, and is infinitely large when horizontal stratification is not present, i.e. far from the fronts. Further, in a stable stratification 
$(N_z^2>0)$, the sign of 
$\tan \varTheta /\tan \varGamma$ would be the same as the sign of 
$\varPhi$. Given the possibility of 
$\tan \varTheta /\tan \varGamma$ taking any value in 
$(-\infty ,\infty )$, it is important to accurately estimate the effective 
$Sc$ in various settings. As noted in figure 8(a), the departure of 
$Sc$ from unity can both increase the range of unstable 
$\tan \varTheta /\tan \varGamma$ and introduce the possibility of oscillatory instability.
It is also worthwhile to discuss the importance of curvature effects in various oceanic settings. In mesoscale eddies, for example, the current speed and radius are, respectively, of the order of 
$1\ \textrm {m}\ \textrm {s}^{-1}$ and 
$100$ km. These values give estimates for 
$V/r$ and 
$\partial V/\partial r$ as 
$10^{-5}\ \textrm {s}^{-1}$. For cyclonic eddies, as the latitude varies from 
$0$ (equator) to 
$90^{\circ }$ (poles), the curvature parameter in (4.1a,b) varies from 
$1$ to 
$0.064$. Correspondingly, for anticyclonic eddies, it varies from 
$1$ at the equator to 
$\infty$ at a latitude of around 
$4^\circ$N, and then becomes strongly negative slightly northward; 
$R_c$ then reaches a value of 
$-1$ at around 
$8^{\circ }$N, before asymptotically increasing towards 
$-0.074$ at 
$90^{\circ }$N. The latitude at which 
$R_c$ reverses sign will be larger for stronger current speeds or smaller eddies. As shown in figure 10, 
$R_c$ can significantly affect the stability criterion, and a curvature-effect-free theory would therefore be error prone.
Another feature of importance in sub-mesoscale oceanic eddies is their finite extent in both the horizontal and vertical directions. Specifically, intra-thermocline eddies tend to be around 10–20 km and 200–600 m in horizontal and vertical extents, respectively (Nguyen et al. Reference Nguyen, Hua, Schopp and Carton2012), and are often modelled as pancake vortices. Depending on the vorticity distribution, Reynolds number, Froude number, Burger number and Schmidt number, these pancake vortices can be in different instability regimes, including barotropic and baroclinic instabilities (Yim et al. Reference Yim, Billant and Ménesguen2016; Storer, Poulin & Ménesguen Reference Storer, Poulin and Ménesguen2018), short-wavelength instabilities such as symmetric and double-diffusive instabilities (Negretti & Billant Reference Negretti and Billant2013; Yim & Billant Reference Yim and Billant2016). Some of these instabilities are recognized as driving mechanisms for density layering observed in oceanic eddies (Hua et al. Reference Hua, Ménesguen, Le Gentil, Schopp, Marsset and Aiki2013). Of specific relevance to our study are the results of Yim & Billant (Reference Yim and Billant2016), who performed global stability analysis to investigate short-wavelength and other instabilities in Gaussian pancake vortices. Yim & Billant (Reference Yim and Billant2016) also presented results for 
$Sc \ne 1$, to show that oscillatory instability (McIntyre Reference McIntyre1970) occurs at 
$Sc = 700$ but not at 
$Sc = 7$. To recall from our results (figure 8), for 
$Sc>1$, the entire centrifugally unstable regime is conducive to the oscillatory instability as well; the Gaussian pancake vortex, however, contains centrifugally stable and unstable regions. It would therefore be insightful to evaluate the local stability criteria (and identify the most unstable perturbations corresponding to monotonic and oscillatory instabilities) on all the streamlines in the Gaussian pancake vortex considered by Yim & Billant (Reference Yim and Billant2016), to subsequently establish a relation between the local and global stability results.
5. Conclusions
In this study, we performed a local stability analysis of a stratified, axisymmetric vortex in thermal wind balance, which represents a model for large-scale vortices in the ocean and atmosphere. The results are also relevant for fronts that appear in the vicinity of eddies. Using the local stability approach, which is alternate to the conventional normal mode approach of McIntyre (Reference McIntyre1970), our results extend those of McIntyre (Reference McIntyre1970) to include both curvature effects and centrifugally unstable regimes. In addition, the limit of small but finite viscous effects are easily recovered with the non-dimensionalization used in our study; this allowed us to explore the entire plane of perturbation characteristics at base flow conditions that were far from the neutral stability boundaries. The local stability framework was first shown to recover the well-known symmetric instability in the inviscid limit (§ 3.2). The study then focused on the effects of Schmidt number, which represents the ratio between momentum and mass diffusivities.
In the viscous, double-diffusive regime, the instability characteristics were shown to depend on five non-dimensional parameters: two base flow parameters, two perturbation parameters and 
$Sc$ (§ 3.3). Two different instabilities, monotonic and oscillatory, were shown to occur for 
$Sc\ne 1$. Monotonic instability, which is the same as symmetric instability at 
$Sc = 1$, can occur for inviscidly stable (absence of symmetric instability) parameter values if 
$Sc\ne 1$; correspondingly, inviscidly unstable parameter values can become stable for 
$Sc\ne 1$. Oscillatory instability, characterized by a non-zero imaginary part in the growth rate, occurs only for 
$Sc\ne 1$. The modification of monotonic and oscillatory instability regimes due to 
$Sc\ne 1$ were then discussed in detail in the entire parameter space, including the centrifugally unstable regime (§§ 3.3.1–3.3.2). For a given base flow and 
$Sc$, the monotonic instability dominates the oscillatory instability in terms of the growth rate, although both instabilities cannot occur simultaneously for the same perturbations. Neutral stability boundaries on the plane of 
$Sc$ and a modified gradient Richardson number 
$(\tan \varTheta /\tan \varGamma )$ were then identified for both monotonic and oscillatory instabilities (§ 3.3.3). The monotonic instability was shown to significantly extend into the 
$\tan \varTheta /\tan \varGamma >1$ region when 
$Sc$ moves away from unity (for both 
$Sc<1$ and 
$Sc>1$). The oscillatory instability, which does not occur at 
$Sc = 1$, was shown to occupy a significant portion of the monotonically unstable region in 
$\tan \varTheta /\tan \varGamma >0$ for 
$Sc<1$. For 
$Sc>1$, oscillatory instability extends into the entirety of inviscidly unstable 
$(-\infty <\tan \varTheta /\tan \varGamma <1)$ region, with a small excursion to the 
$\tan \varTheta /\tan \varGamma >1$ region.
In § 3.4, it was shown that it is possible to further reduce the number of governing non-dimensional parameters to three by choosing an alternate non-dimensionalization of (3.1), albeit at the cost of not separating the base flow and perturbation characteristics (§ 3.4). Each of the two newly identified parameters other than 
$Sc$ represent a combination of base flow and perturbation characteristics. Neutral stability boundaries for both monotonic and oscillatory instabilities were then analytically derived in terms of the aforementioned three parameters. Additionally, it was also shown that monotonic and oscillatory instabilities cannot simultaneously occur for the same values of the parameters. Finally, in §§ 4.1–4.2, the implications of our overall results were discussed from two different view points. Specifically, the curvature effects were quantified in terms of a non-dimensional parameter, and it was shown that ignoring curvature effects can be misleading in typical oceanic settings. In addition to curvature effects, non-traditional effects (inclusion of the complete Coriolis force, Gerkema et al. Reference Gerkema, Zimmerman, Maas and Van Haren2008) could also be important, especially near the equator. Indeed, recent studies have shown that the horizontal component of background rotation can significantly affect the symmetric instability (Itano & Maruyama Reference Itano and Maruyama2009; Kloosterziel, Carnevale & Orlandi Reference Kloosterziel, Carnevale and Orlandi2017; Zeitlin Reference Zeitlin2018).
In the future, it would be worthwhile investigating the manifestation of the instabilities discussed in this paper in direct numerical simulations, including the nonlinear flow states that would result from these instabilities. To explore short-wavelength instability regimes in realistic oceanic settings, the instability criteria derived in this study could be evaluated in field data (in situ and satellite measurements). The results in this study are also potentially useful in interpreting outputs from large-scale climate models, which often use some constant 
$Sc$ value, unity or otherwise (Danabasoglu et al. Reference Danabasoglu, Bates, Briegleb, Jayne, Jochum, Large, Peacock and Yeager2012). An interesting follow-up to this study would be to investigate the effects of non-isotropic diffusivity, which is often the case in stratified oceanic settings, and triple-diffusive effects resulting from the different diffusivities of momentum, temperature and salinity. Finally, while the current study focused on short-wavelength instabilities, other instabilities such as Kelvin–Helmholtz and baroclinic instabilities have to be considered too to construct an overall instability regimes diagram.
Acknowledgements
The authors would like to thank A. Tandon for his input, and B. Fox-Kemper for pointing us to some of the recent literature.
Funding
The authors also acknowledge the Monsoon Mission Grant MM/2014/IND-002, the FIST grant SR/FST/ET-II/2017/109 and the IIT Madras research initiative titled Geophysical Flows Lab.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Classical double-diffusive instability (salt fingering)
In this appendix, we show that the local stability framework recovers the classical salt fingering instability using a 
${O}(\delta ^2)$ scaling for the temperature, salinity and momentum diffusivities. For a quiescent base flow (
$U_B = 0$) with a salinity profile 
$S_B(z)$ and a temperature profile 
$T_B(z)$, the local stability equations can be written as
$$\begin{gather} \frac{\textrm{d}\boldsymbol{a}}{\textrm{d}t} ={-}b\hat{e}_z + b\boldsymbol{k}\frac{k_z}{|\boldsymbol{k}|^2} - c_\nu|\boldsymbol{k}|^2\boldsymbol{a}, \end{gather}$$
$$\begin{gather}\frac{\textrm{d}s}{\textrm{d}t} ={-}\boldsymbol{a}\boldsymbol{\cdot}\boldsymbol{\nabla} S_B - c_s|\boldsymbol{k}|^2s, \end{gather}$$
$$\begin{gather}\frac{\textrm{d}\tau}{\textrm{d}t} ={-}\boldsymbol{a}\boldsymbol{\cdot}\boldsymbol{\nabla} T_B - c_\tau|\boldsymbol{k}|^2\tau, \end{gather}$$
where 
$s,\tau$ denote the leading-order perturbation amplitudes in salinity and temperature, respectively (salinity and temperature perturbations are written as a series similar to (2.10)–(2.12)). The diffusivities of momentum, salinity and temperature are given by 
$\nu =c_\nu \delta ^2,\kappa _s=c_s\delta ^2$ and 
$\kappa _\tau =c_\tau \delta ^2$, where 
$c_\nu , c_s$ and 
$c_\tau$ are constants. Here, 
$b$ refers to the leading-order buoyancy perturbation amplitude
\begin{equation} b = g[\beta s - \alpha\tau], \end{equation}
where 
$\beta ,\alpha$ are expansion coefficients in salinity and temperature, respectively. Writing (A1) in component form, and substituting (A4), (A1)–(A3) become
$$\begin{gather} \frac{\textrm{d}a_x}{\textrm{d}t} = (\beta s - \alpha\tau)g\frac{k_xk_z}{|\boldsymbol{k}|^2} - c_\nu|\boldsymbol{k}|^2a_x, \end{gather}$$
$$\begin{gather}\frac{\textrm{d}a_y}{\textrm{d}t} = (\beta s - \alpha\tau)g\frac{k_yk_z}{|\boldsymbol{k}|^2} - c_\nu|\boldsymbol{k}|^2a_y, \end{gather}$$
$$\begin{gather}\frac{\textrm{d}a_z}{\textrm{d}t} ={-}(\beta s - \alpha\tau)g + (\beta s - \alpha\tau)g\frac{k_z^2}{|\boldsymbol{k}|^2} - c_\nu|\boldsymbol{k}|^2a_z, \end{gather}$$
$$\begin{gather}\frac{\textrm{d}s}{\textrm{d}t} ={-}a_z\frac{\textrm{d}S_B}{\textrm{d}z} - c_s|\boldsymbol{k}|^2s, \end{gather}$$
$$\begin{gather}\frac{\textrm{d}\tau}{\textrm{d}t} ={-}a_z\frac{\textrm{d}T_B}{\textrm{d}z} - c_\tau|\boldsymbol{k}|^2\tau, \end{gather}$$
where 
$\boldsymbol {a} = a_x\hat {e}_x + a_y\hat {e}_y + a_z\hat {e}_z$, and 
$\boldsymbol {k} = k_x\hat {e}_x + k_y\hat {e}_y + k_z\hat {e}_z$.
For a quiescent base flow, (2.13) suggests that any wave vector remains constant. Thus, (A5)–(A9) are a set of constant coefficient, linear, first-order ordinary differential equations. The non-trivial growth rates 
$\lambda$ can then be shown to be governed by
\begin{equation} \hat{\lambda}^3 + (1+Pr + Pr_s)\hat{\lambda}^2 + [\hat{T}- \hat{S} + Pr_s + Pr + Pr_s Pr]\hat{\lambda} - \hat{S}Pr + \hat{T} + Pr_s Pr = 0, \end{equation}
where 
$\hat {\lambda } = \lambda /c_s|\boldsymbol {k}|^2$, 
$Pr_s = c_\nu /c_s$, 
$Pr = c_\tau /c_s$, 
$\hat {S} = \beta g(1-k_z^2\mathbin {{/|\boldsymbol {k}|^2}})(\textrm {d}S_B/\textrm {d}z)/(c_s^2|\boldsymbol {k}|^4)$, 
$\hat {T} = \alpha g (1-k_z^2\mathbin {{/|\boldsymbol {k}|^2}})(\textrm {d}T_B/\textrm {d}z)/(c_s^2|\boldsymbol {k}|^4)$. Equation (A10) is identical to the well-known growth rate equation for double-diffusive salt fingering instability (Turner Reference Turner1973). In summary, we have shown that the 
${O}(\delta ^2)$ scaling for diffusivities in the local stability equations allows us to recover the growth rates associated with the double-diffusive salt fingering instability.
Appendix B. Neutral stability boundaries in terms of (
$\alpha _1$, $\alpha _2$, $\beta$, $Sc$)
In this appendix, we derive the neutral stability boundaries on the (
$\alpha _1-\beta$), (
$\alpha _2-\beta$) plane for various 
$Sc$. As shown in (3.20), 
${\lambda }^* = (\lambda +c_\nu |\boldsymbol {k}|^2)/(c_\nu |\boldsymbol {k}|^2)$ is governed by
\begin{equation} {\lambda}^{*3} + \left(\frac{1-Sc}{Sc}\right){\lambda}^{*2} + [\alpha_1+\alpha_2-2\beta]{\lambda}^* + \left(\frac{1-Sc}{Sc}\right)[\alpha_1-\beta]=0, \end{equation}
with the neutral stability boundaries occurring at 
$\textrm {Re}(\lambda ^*) = 1$. To identify these boundaries, we assume 
$\lambda ^* = 1+{{\varDelta }} + \textrm {i}c$ (
${{\varDelta }}$, 
$c$ are real), where 
${{\varDelta }}$ is small in the neighbourhood of neutral stability boundaries. Substituting this 
$\lambda ^*$ in (B1), and ignoring terms that are 
${O}({{\varDelta }}^2$) or higher, we get the imaginary part of the equation as
\begin{equation} -\textrm{i}c^3+3\textrm{i}c+ 6\textrm{i}c{{{\varDelta}}} + \left(\frac{1}{Sc}-1\right)2\textrm{i}c(1+{{{\varDelta}}}{}) + \textrm{i}c(\alpha_1+\alpha_2-2\beta) = 0. \end{equation}
Solving (B2) for 
$c$, we get
\begin{equation} c=0 \text{ or } c^2 = 1+\frac{2}{Sc}+\alpha_1 + \alpha_2-2\beta + 2{{{\varDelta}}}\left(2+\frac{1}{Sc}\right). \end{equation} For 
$c=0$, the real part of (B1) reduces to
\begin{equation} \frac{\alpha_1-\beta}{Sc} + \alpha_2-\beta+\frac{1}{Sc} + {{{\varDelta}}}\left[\alpha_1+\alpha_2-2\beta+1+\frac{2}{Sc}\right]=0, \end{equation}
which shows that 
$\lambda ^* = 1$ is a root of (B1) on the line
\begin{equation} \frac{\alpha_1-\beta}{Sc} + \alpha_2-\beta+\frac{1}{Sc} = 0, \end{equation}
if 
$\alpha _1 + \alpha _2 - 2\beta + 1 + 2/Sc \ne 0$. Furthermore, the line (B5) is a neutral stability boundary only if the other two roots of (B1) are stable on it. The roots of (B1) other than 
$\lambda ^*_1 = 1+{{\varDelta }}$ are given by
\begin{equation} {\lambda}^*_{2,3} ={-}\frac{1}{2}\left(\frac{1}{Sc}+{{{\varDelta}}}\right) \pm\frac{1}{2}\sqrt{\left(\frac{1}{Sc}+{{{\varDelta}}}\right)^2 - 4\left(\frac{1+{{{\varDelta}}}}{Sc} + {{{\varDelta}}} + \alpha_1+\alpha_2-2\beta\right)}. \end{equation}
The roots 
${\lambda }^*_{2,3}$ are real (with 
${{\varDelta }} = 0$) if
\begin{equation} s(\alpha_1-\beta) > s\frac{1}{4Sc(Sc-1)}, \end{equation}and stable (less than unity) if
\begin{equation} s(\alpha_1-\beta) < s\frac{1+Sc}{1-Sc}, \end{equation}
where 
$s = \textrm {sgn}(1-Sc)$. If conditions (B7) and (B8) are satisfied, then the eigenvalue that goes from stable to unstable across the neutral stability boundary in (B5) is real, while the remaining two stable eigenvalues are also real in the neighbourhood of the neutral stability boundary. For 
$\lambda ^*_1= 1$, the other two eigenvalues (B6) are complex on the boundary in (B5) if
\begin{equation} s(\alpha_1-\beta) < s\frac{1}{4Sc(Sc-1)}, \end{equation}
and are always stable since their real part is negative for 
${{\varDelta }}=0$. In other words, if (B9) is satisfied, the straight line in (B5) is a neutral stability boundary, with the eigenvalue that goes from stable to unstable across the neutral stability boundary being real and the remaining two stable eigenvalues being complex. In summary, (B5) represents the monotonic (since 
$c = 0$) neutral stability boundary for either conditions (B7) and (B8) being satisfied or condition (B9) being satisfied.
For the case of 
$c^2 = 1+{2}/{Sc}+\alpha _1+ \alpha _2-2\beta + 2{{\varDelta }}(2+{1}/{Sc})$ in (B3), the real part of (B1) becomes
\begin{align}
& 2(\alpha_1-\beta) +
(\alpha_2-\beta)\left(1+\frac{1}{Sc}\right)+2\left(1+\frac{2}{Sc}+\frac{1}{Sc^2}\right) \notag\\
&\quad ={-}2{{{\varDelta}}}\left(\alpha_1+ \alpha_2-2\beta + 5 +
\frac{1}{Sc^2} + \frac{6}{Sc} \right),
\end{align}
provided 
$1+{2}/{Sc}+\alpha _1 + \alpha _2-2\beta + 2{{\varDelta }}(2+{1}/{Sc})>0$. Thus, the line
\begin{equation} 2(\alpha_1-\beta) + (\alpha_2-\beta)\left(1+\frac{1}{Sc}\right)+2\left(1+\frac{2}{Sc}+\frac{1}{Sc^2}\right)=0 \end{equation}
is a neutral stability boundary with 
$c^2>0$ and 
${{\varDelta }}=0$ (as the third root 
$\lambda ^*_3 = -1-1/Sc$ is stable) if
\begin{equation} s(\alpha_1-\beta) < s\frac{1+Sc}{1-Sc}. \end{equation}
The condition in (B12) is obtained based on the requirement of 
$c^2>0$ on the line in (B11). In summary, (B11) is the oscillatory neutral stability boundary for those (
$\alpha _1-\beta$) satisfying (B12).
Appendix C. Monotonic–oscillatory instability interface
In this appendix, we derive closed form expressions for the interface between monotonic and oscillatory instabilities within the unstable region on the (
$\alpha _1-\beta$), (
$\alpha _2-\beta$) plane (see figure 9). Approaching the interface from the oscillatory instability side, we seek solutions of the form 
$\lambda ^* = (a+\textrm {i}b)$, 
$(a-\textrm {i}b)$, 
$c$ for (3.20), where 
$a$, 
$b$ and 
$c$ are all real. Comparing the coefficients of (3.20) and the cubic given by 
$(\lambda ^*-(a+ \textrm {i}b))(\lambda ^*-(a-\textrm {i}b))(\lambda ^*-c)$, we get 
$2a+c= - (1-Sc)/Sc$. If monotonic and oscillatory instabilities were to occur simultaneously, we require 
$a>1$ and 
$c>1$. (recall that 
$\textrm {Re}{(\lambda ^*)} > 1$ for instability). As a result, 
$Sc$ has to satisfy 
$(1-Sc)/Sc>-3$, which is not possible for 
$Sc>0$. Hence, for all real values of 
$\alpha _1$, 
$\alpha _2$, 
$\beta$ and 
$Sc$, the monotonic and oscillatory instabilities cannot occur simultaneously.
Now, we seek the interface between monotonic and oscillatory instabilities as the boundary where the unstable roots 
$a+\textrm {i}b$ have 
$b\to 0$ while we approach it from the oscillatory instability side. Upon substituting 
$\lambda ^* = a + \textrm {i}b$, the real and imaginary parts of (3.20) become
$$\begin{gather} a^3 - 3ab^2 + \frac{1-Sc}{Sc}(a^2-b^2) + [\alpha_1+\alpha_2 - 2\beta]a + \frac{1-Sc}{Sc}(\alpha_1-\beta) = 0, \end{gather}$$
$$\begin{gather}-\textrm{i}b^3 + 3a^2\textrm{i}b + \frac{1-Sc}{Sc}2\textrm{i}ab + [\alpha_1+\alpha_2-2\beta]\textrm{i}b = 0. \end{gather}$$
Solving for 
$b$ in (C2), and requiring the 
$b\ne 0$ solution to tend towards zero gives
\begin{equation} a_\pm \to -\frac{1-Sc}{Sc} \pm \sqrt{\frac{(1-Sc)^2}{9Sc^2} - \frac{\alpha_1+\alpha_2-2\beta}{3}}. \end{equation}
For instability, we require 
$a>1$, which can be satisfied only by 
$a_+$ in (C3). Requiring 
$a_+>1$ gives
\begin{equation} \alpha_1+\alpha_2-2\beta <{-}1-\frac{2}{Sc}. \end{equation}
Noting that 
$\alpha _1+\alpha _2-2\beta = -1-{2}/{Sc}$ intersects the monotonic and oscillatory neutral stability boundaries (3.22)–(3.23) at a common location (
$\alpha _1-\beta = (Sc+1)/(Sc-1)$, 
$\alpha _2 - \beta = -2/Sc(Sc-1))$, and comparing their slopes, it can be shown that the monotonic–oscillatory instability interface satisfies (C4). The equation describing the monotonic–oscillatory instability interface can then be derived by substituting 
$a_+$ in (C1) with 
$b=0$ to get
\begin{align} &\frac{2(1-Sc)^3}{9Sc^3} + \frac{1-Sc}{3Sc}(2(\alpha_1-\beta) +(\alpha_2-\beta)) \nonumber\\ &\quad -\left(\frac{2(1-Sc)^2}{9Sc^2} - \frac{2(\alpha_1+\alpha_2-2\beta)}{3}\right) \sqrt{\frac{(1-Sc)^2}{9Sc^2} - \frac{\alpha_1+\alpha_2-2\beta}{3}} = 0. \end{align}Appendix D. Neutral stability boundaries in terms of physical parameters
To derive the neutral stability boundaries on the plane of base flow parameters, we first write 
$\alpha _1-\beta$ and 
$\alpha _2-\beta$ (§ 3.4) in terms of the physical parameters in (3.4a,b)–(3.5a,b)
\begin{equation} \alpha_1-\beta = \hat{\varPhi}\frac{2\cos^2\theta - \tan\varGamma\sin 2\theta}{\tilde{c}^2}, \quad \alpha_2-\beta = \hat{\varPhi}\frac{\tan\varGamma(2\tan\varTheta\sin^2\theta - \sin 2\theta)}{\tilde{c}^2}. \end{equation}
In this appendix, we assume 
$\hat {\varPhi } = 1$, i.e. the centrifugally stable regime.
In the neighbourhood of the monotonic neutral stability boundary on the (
$\alpha _1-\beta$), (
$\alpha _2-\beta$) plane, the growth rate 
${{\varDelta }}$ is governed by (B4). For given base flow parameters, we first maximize 
${{\varDelta }}$ with respect to the perturbation parameters, and then identify the 
${{\varDelta }} = 0$ boundary on the plane of base flow parameters. To maximize 
${{\varDelta }}$ with respect to 
$\theta$, we differentiate (B4) with respect to 
$\theta$ to get
\begin{align} &\frac{1}{Sc}\left(\frac{\partial \alpha_1}{\partial\theta}-\frac{\partial\beta}{\partial\theta}\right) + \frac{\partial \alpha_2}{\partial\theta}-\frac{\partial\beta}{\partial\theta} + \frac{\partial{{\varDelta}}{}}{\partial \theta}\left[\alpha_1+\alpha_2-2\beta+1+\frac{2}{Sc}\right]\notag\\ &\quad + {{{\varDelta}}}{}\left[\frac{\partial\alpha_1}{\partial \theta}+\frac{\partial \alpha_2}{\partial \theta}-2\frac{\partial\beta}{\partial \theta}\right]=0, \end{align}
and substitute 
$\partial {{\varDelta }}/\partial \theta = 0$. Furthermore, 
${{\varDelta }} = 0$ on the neutral stability boundary we are seeking, thus reducing (D2) to
\begin{equation} \frac{1}{Sc}\left(\frac{\partial \alpha_1}{\partial\theta}-\frac{\partial\beta}{\partial\theta}\right) + \frac{\partial \alpha_2}{\partial\theta}-\frac{\partial\beta}{\partial\theta}=0. \end{equation}
Substituting the expressions in (D1) in (D3) gives the expression for the most unstable 
$\theta$ for monotonic instability as
\begin{equation} \tan2\theta^* = \frac{(Sc+1)\tan\varGamma}{Sc\tan\varTheta\tan\varGamma-1}. \end{equation}
Recalling from § 3.3.1 that 
$\tilde {c} = 0$ on the neutral stability boundaries in the plane of (
$Sc,(\tan \varTheta /\tan \varGamma )$) and substituting (D4) in (B5), we obtain the monotonic neutral stability boundary as
\begin{equation} \frac{\tan\varTheta}{\tan\varGamma} = \frac{(1+Sc)^2}{4Sc}. \end{equation} In the neighbourhood of the oscillatory neutral stability boundary on the (
$\alpha _1-\beta$), (
$\alpha _2-\beta$) plane, the growth rate 
${{\varDelta }}{}$ is governed by (B10). For given base flow parameters, we first maximize 
${{\varDelta }}{}$ with respect to the perturbation parameters, and then identify the 
${{\varDelta }}{} = 0$ boundary on the plane of base flow parameters. To maximize 
${{\varDelta }}{}$ with respect to 
$\theta$, we differentiate (B10) with respect to 
$\theta$ to get
\begin{align} &2\left(\frac{\partial \alpha_1}{\partial\theta}-\frac{\partial\beta}{\partial\theta}\right) + \left(\frac{\partial \alpha_2}{\partial\theta}-\frac{\partial\beta}{\partial\theta}\right) \left(1+\frac{1}{Sc}\right) + 2\frac{\partial{{{\varDelta}}}{}}{\partial \theta}\left[\alpha_1+\alpha_2-2\beta+5+\frac{1}{Sc^2}+\frac{6}{Sc}\right] \nonumber\\ &\quad +2{{{\varDelta}}}{}\left[\frac{\partial\alpha_1}{\partial \theta}+\frac{\partial \alpha_2}{\partial \theta}-2\frac{\partial\beta}{\partial \theta}\right]=0, \end{align}
and substitute 
$\partial {{\varDelta }}{}/\partial \theta = 0$. Furthermore, 
${{\varDelta }}{} = 0$ on the neutral stability boundary we are seeking, thus reducing (D6) to
\begin{equation} 2\left(\frac{\partial \alpha_1}{\partial\theta}-\frac{\partial\beta}{\partial\theta}\right) + \left(\frac{\partial \alpha_2}{\partial\theta}-\frac{\partial\beta}{\partial\theta}\right) \left(1+\frac{1}{Sc}\right)=0. \end{equation}
Substituting the expressions in (D1) in (D7) gives the expression for the most unstable 
$\theta$ for oscillatory instability as
\begin{equation} \tan2\theta^*_o = \frac{3Sc+1}{(1+Sc)\tan\varTheta-4Sc\cot\varGamma}. \end{equation}
Recalling from § 3.3.1 that 
$\tilde {c} = 0$ on the neutral stability boundaries in the plane of (
$Sc,(\tan \varTheta /\tan \varGamma )$) and substituting (D8) in (B11), we obtain the oscillatory neutral stability boundary as
\begin{equation} \frac{\tan\varTheta}{\tan\varGamma} = \frac{(3Sc+1)^2}{8Sc(Sc+1)}. \end{equation}
Open access
