© J. Park et al. 2021

1. Introduction

Stellar rotation is one of the key physical processes to build a modern picture of stellar evolution (e.g., Maeder 2009, and references therein). Indeed, it triggers the transport of angular momentum and of chemicals, which drives the rotational and chemical evolution of stars, respectively (e.g., Zahn 1992; Maeder & Zahn 1998; Mathis & Zahn 2004). This has a major impact on the late stages of their evolution (e.g., Hirschi et al. 2004), their magnetism (e.g., Brun & Browning 2017), their winds and mass losses (e.g., Ud-Doula et al. 2009; Matt et al. 2015), and the interactions with their planetary and galactic environment (e.g., Gallet et al. 2017; Strugarek et al. 2017).

A robust ab-initio evaluation of the strength of each (magneto-)hydrodynamical mechanism that transports momentum and chemicals is thus mandatory to understand astrophysical observations. In particular, our knowledge of the internal rotation of stars and its evolution has been revolutionized thanks to space-based asteroseismology with the Kepler space mission (NASA; Aerts et al. 2019, and references therein). It has been proven that stars mostly host weak differential rotation in the whole Hertzsprung-Russell diagram, such as our Sun. Stellar interiors are thus the seat of efficient mechanisms that transport angular momentum all along their evolution. These mechanisms have not been identified yet even if several candidates have been proposed, such as stable and unstable magnetic fields (e.g., Moss 1992; Charbonneau & MacGregor 1993; Spruit 1999, 2002; Fuller et al. 2019), stochastically-excited internal gravity waves (e.g., Talon & Charbonnel 2005; Rogers 2015; Pinçon et al. 2017), and mixed gravito-acoustic modes (Belkacem et al. 2015a,b). In this framework, improving our knowledge of the hydrodynamical turbulent transport induced by the instabilities of the stellar differential rotation is mandatory since it has been recently proposed as another potential efficient mechanism to transport angular momentum (Barker et al. 2020; Garaud 2020) while it constitutes one of the cornerstones of the theory of the rotational transport and mixing along the evolution of stars (Zahn 1992).

Since stellar radiation zones are stably stratified, rotating regions, it is expected that the turbulent transport triggered there by the instabilities of the vertical and horizontal shear of the differential rotation should be anisotropic (Zahn 1992). As pointed out in Mathis et al. (2018) and Park et al. (2020), it is the vertical shear instabilities that have received important attention in the literature, in particular by taking into account the impact of the high thermal diffusion in stellar radiation zones (Zahn et al. 1983; Lignières 1999) and the interactions with the horizontal turbulence induced by the horizontal shear (Talon & Zahn 1997). The obtained prescriptions, mainly derived using phenomenological modelings, have been broadly implemented in state-of-the-art evolution models of rotating stars (e.g., Ekström et al. 2012; Marques et al. 2013; Amard et al. 2019). They are now tested using direct numerical simulations devoted to the stellar regime (Prat & Lignières 2013, 2014; Garaud et al. 2017; Prat et al. 2016; Gagnier & Garaud 2018; Kulenthirarajah & Garaud 2018). The horizontal turbulence induced by horizontal gradients of the differential rotation has only been examined in the stellar regime in a few works, which again use phenomenological arguments (Zahn 1992; Maeder 2003), results from lab experiments studying differentially rotating flows (Richard & Zahn 1999; Mathis et al. 2004), and first devoted numerical simulations (Cope et al. 2020). The systematic study of the combined effect of stable stratification, rotation, and thermal diffusion in the stellar context has only been recently undertaken.

In Park et al. (2020), we thus examined the behavior of shear instabilities sustained by a horizontal shear as a function of stratification, rotation, and thermal diffusion. However, in this first work, we neglected the horizontal projection of the rotation vector and the corresponding terms of the Coriolis acceleration, following the traditional approximation of rotation (TAR) which is often used to describe geophysical and astrophysical flows in stably stratified, rotating regions (Eckart 1960). However, recent works have disputed this approximation as it can fail in some cases. For instance, the dynamics of near-inertial waves is strongly influenced by the nontraditional effects in such a way that properties of the wave reflection and associated mixing are largely modified (Gerkema & Shrira 2005; Gerkema et al. 2008). Moreover, couplings between gravito-inertial waves, inertial waves, and wave-induced turbulence in stars can only be properly treated when the full Coriolis acceleration is considered (Mathis et al. 2014). Finally, Zeitlin (2018) demonstrated analytically that the instability of a linear shear flow can be significantly modified if the full Coriolis acceleration is considered. Nontraditional effects can be particularly important for stellar structure configurations where the Coriolis acceleration can compete with the Archimedean force in the direction of both entropy and chemical stratification, for instance during the formation of the radiative core of pre-main-sequence, low-mass stars or in the radiative envelope of rapidly rotating upper main-sequence stars. These regimes should be treated properly to build robust one- or two-dimensional (1D or 2D) secular models of the evolution of rotating stars (e.g., Ekström et al. 2012; Amard et al. 2019; Gagnier et al. 2019).

In this paper, we, therefore, continue our previous work that examined the effect of thermal diffusion on horizontal shear instabilities in stably stratified rotating stellar radiative zones (Park et al. 2020), but we now consider the full Coriolis acceleration. In particular, we investigate how the full Coriolis acceleration modifies the dynamics of two types of shear instabilities: the inflectional and inertial instabilities. In Sect. 2, we formulate linear stability equations derived from the Navier-Stokes equations when the fluid is stably stratified and thermally diffusive in the rotating nontraditional f-plane where both vertical and horizontal components of the rotation vector are taken into account. We consider a base shear flow in a hyperbolic tangent form as a canonical shear profile used in previous studies (Schmid & Henningson 2001; Deloncle et al. 2007; Griffiths 2008; Arobone & Sarkar 2012). In Sect. 3, we provide general numerical results on the inflectional and inertial instabilities. In Sect. 4, we use the Wentzel-Kramers-Brillouing-Jeffreys (WKBJ) approximation in the asymptotic limit of large vertical wavenumbers to investigate how the inertial instability is modified by nontraditional effects in the asymptotic nondiffusive and high-diffusivity cases. In Sect. 5, we analyze in detail numerical results in wide ranges of parameters for both the inflectional and inertial instabilities. In Sects. 6 and 7, we propose expressions for the horizontal turbulent viscosity and the characteristic time of the turbulent transport, which can be possibly applied to the 1D- and 2D-secular evolution models of rotating stars. Finally, in Sect. 8, we provide our conclusions and discussions about the nontraditional effects on the horizontal shear instabilities, and we propose the perspectives of this work for astrophysical flows.

2. Problem formulation

2.1. Navier-Stokes equations and base steady state

We consider the Navier-Stokes and heat transport equations under the Boussinesq approximation in a rotating frame. We define the local Cartesian coordinates (x, y, z) with x the longitudinal coordinate, y the latitudinal coordinate, and z the vertical coordinate:

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathit{u}=0,\end{array}$$(1)

$$\begin{array}{c}\hfill \frac{\partial \mathit{u}}{\partial t}+(\mathit{u}·\mathrm{\nabla })\mathit{u}+\mathit{f}\times \mathit{u}=-\frac{1}{{\rho }_0}\mathrm{\nabla }p-{\alpha }_{\mathrm{T}}\mathrm{\Theta }\mathit{g}+{\nu }_0{\mathrm{\nabla }}^2\mathit{u},\end{array}$$(2)

$$\begin{array}{c}\hfill \frac{\partial \mathrm{\Theta }}{\partial t}+\mathit{u}·\mathrm{\nabla }\mathrm{\Theta }={\kappa }_0{\mathrm{\nabla }}^2\mathrm{\Theta },\end{array}$$(3)

where u = (u,v,w) is the velocity, p is the pressure, Θ is the temperature deviation from the reference temperature T₀, f = (0,f_h, 0,f_v, 0) is the Coriolis vector with the horizontal (latitudinal) Coriolis component f_h, 0 = 2Ω₀ sin θ and the vertical Coriolis component f_v, 0 = 2Ω₀ cos θ where ${\mathrm{\Omega }}_0={(f_{\mathrm{h},0}^2+f_{\mathrm{v},0}^2)}^{1/2}/2$ is the stellar rotation rate, θ is the colatitude, ρ₀ is the reference density, g = (0, 0, −g) is the gravity, ν₀ is the reference viscosity, κ₀ is the reference thermal diffusivity, α_T is the thermal expansion coefficient that assumes a linear relation between the density deviation and the temperature deviation Θ, and ∇² = ∂²/∂x² + ∂²/∂y² + ∂²/∂z² denotes the Laplacian operator. Figure 1 illustrates the local coordinate system of the horizontal shear flow in the radiative zone when the nontraditional f-plane is considered at a given colatitude θ. For the base state, we consider a canonical example of the steady horizontal shear flow U = (U(y),0,0) in a hyperbolic tangent form:

$$\begin{array}{c}\hfill U(y)=U_{0}tanh\left(\frac{y}{L_0}\right),\end{array}$$(4)

Fig. 1. (a) Schematic of the radiative and convective zones, colored as yellow and orange, respectively, for the configuration of low-mass stars rotating with angular speed Ω0. (b) Horizontal shear flow U(y) in a local nontraditional f-plane at a colatitude θ in the radiative zone; zc denotes the transition altitude between the radiative and convective zones.

where U₀ and L₀ are the reference velocity and length scales, respectively. This hyperbolic tangent profile has an inflection point at y = 0. Therefore, it has been considered in many previous stability studies (see e.g., Schmid & Henningson 2001; Deloncle et al. 2007; Arobone & Sarkar 2012) to investigate the inflectional instability. We also adopt this profile to compare with other studies and extend our understanding of the inflectional instability for a horizontal shear when the full Coriolis acceleration is taken into account. In the presence of the horizontal Coriolis parameter f_h, 0, such a base flow is steady if a thermal-wind balance between U(y) and the base temperature $\overline{\mathrm{\Theta }}(\mathit{y},z)$ is satisfied as follows:

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{T}}g\frac{\partial \overline{\mathrm{\Theta }}}{\partial y}=-f_{\mathrm{h},0}\frac{\partial U}{\partial y}.\end{array}$$(5)

Such a thermal-wind balance has been adopted for cases where the rotation vector and the base shear are not aligned in stratified fluids; for instance, ageostrophic instability in vertical shear flows in stratified fluids in the traditional f-plane (Wang et al. 2014). In addition to the thermal-wind balance, the base temperature $\overline{\mathrm{\Theta }}(\mathit{y},z)$ is considered to be stably stratified with a linearly increasing profile in z; therefore, it has the form

$$\begin{array}{c}\hfill \overline{\mathrm{\Theta }}(y,z)=\frac{\mathrm{\Delta }{\overline{\mathrm{\Theta }}}_0}{\mathrm{\Delta }z}z-\frac{f_{\mathrm{h},0}}{{\alpha }_{\mathrm{T}}g}U(y),\end{array}$$(6)

where $\mathrm{\Delta }{\overline{\mathrm{\Theta }}}_0$ is the difference in base temperature along the vertical distance Δz at a given y.

2.2. Linearized equations

We consider the velocity perturbation $\check{\mathit{u}}=\mathit{u}-\mathit{U}=(\check{u},\check{\mathit{v}},\check{\mathit{w}})$, the pressure perturbation $\check{p}=p-P$ where P is the base pressure, and the temperature perturbation $\check{T}=\mathrm{\Theta }-\overline{\mathrm{\Theta }}$. We nondimensionalize variables in Eqs. (1)–(3) by considering the length scale L₀, the velocity scale U₀, the time scale as t₀ = L₀/U₀, the pressure scale as ${\rho }_0{U}_0^2$, and the temperature scale as $L_0\mathrm{\Delta }{\overline{\mathrm{\Theta }}}_0/\mathrm{\Delta }z$. For infinitesimally small perturbations, we have the following nondimensional linearized equations:

$$\begin{array}{c}\hfill \frac{\partial \check{u}}{\partial x}+\frac{\partial \check{v}}{\partial y}+\frac{\partial \check{w}}{\partial z}=0,\end{array}$$(7)

$$\begin{array}{c}\hfill \frac{\partial \check{u}}{\partial t}+U\frac{\partial \check{u}}{\partial x}+(U\prime -f)\check{v}+\stackrel{\sim }{f}\check{w}=-\frac{\partial \check{p}}{\partial x}+\frac{1}{\mathit{Re}}{\mathrm{\nabla }}^2\check{u},\end{array}$$(8)

$$\begin{array}{c}\hfill \frac{\partial \check{v}}{\partial t}+U\frac{\partial \check{v}}{\partial x}+f\check{u}=-\frac{\partial \check{p}}{\partial y}+\frac{1}{\mathit{Re}}{\mathrm{\nabla }}^2\check{v},\end{array}$$(9)

$$\begin{array}{c}\hfill \frac{\partial \check{w}}{\partial t}+U\frac{\partial \check{w}}{\partial x}-\stackrel{\sim }{f}\check{u}=-\frac{\partial \check{p}}{\partial z}+N^2\check{T}+\frac{1}{\mathit{Re}}{\mathrm{\nabla }}^2\check{w},\end{array}$$(10)

$$\begin{array}{c}\hfill \frac{\partial \check{T}}{\partial t}+U\frac{\partial \check{T}}{\partial x}-\frac{\stackrel{\sim }{f}U\prime }{N^2}\check{v}+\check{w}=\frac{1}{\mathit{Pe}}{\mathrm{\nabla }}^2\check{T},\end{array}$$(11)

where the prime symbol (′) denotes the derivative with respect to y, f = 2Ω cos θ and $\stackrel{\sim }{f}=2\mathrm{\Omega }sin\theta$ are the dimensionless vertical and horizontal Coriolis parameters, respectively, Ω is the dimensionless stellar rotation rate, Re is the Reynolds number

$$\begin{array}{c}\hfill Re=\frac{U_0{L}_0}{{\nu }_0},\end{array}$$(12)

N is the dimensionless Brunt-Väisälä frequency¹

$$\begin{array}{c}\hfill N=\sqrt{\frac{{\alpha }_{\mathrm{T}}g{L}_0^2}{U_0^2}\frac{\mathrm{\Delta }{\overline{\mathrm{\Theta }}}_0}{\mathrm{\Delta }z}},\end{array}$$(13)

and Pe is the Péclet number

$$\begin{array}{c}\hfill Pe=\frac{U_0{L}_0}{{\kappa }_0}.\end{array}$$(14)

In Eq. (11), the term $-\stackrel{\sim }{f}U\prime \widehat{v}/N^2$ on the left-hand side appears due to the thermal-wind balance (5).

To derive linear stability equations, we apply a normal mode expansion to perturbations as

$$\begin{array}{c}\hfill [\check{\mathit{u}},\check{p},\check{T}]=[\widehat{\mathit{u}}(y),\widehat{p}(y),\widehat{T}(y)]exp[\mathrm{i}(k_{x}x+k_{z}z)+\sigma t]+c.c.,\end{array}$$(15)

where $\widehat{\mathit{u}}=(\widehat{u},\widehat{v},\widehat{w})$, $\widehat{p}$ and $\widehat{T}$ are the mode shapes of velocity, pressure, and temperature perturbations, respectively, i² = −1, k_x is the streamwise (longitudinal) wavenumber, k_z is the vertical wavenumber, σ = σ_r + iσ_i is the complex growth rate where the real part σ_r is the growth rate and the imaginary part σ_i is the temporal frequency, and c.c. denotes the complex conjugate. Using the normal mode expansion, we obtain the following linear stability equations

$$\begin{array}{c}\hfill \mathrm{i}{k}_x\widehat{u}+\frac{\mathrm{d}\widehat{v}}{\mathrm{d}y}+\mathrm{i}{k}_z\widehat{w}=0,\end{array}$$(16)

$$\begin{array}{c}\hfill (\sigma +\mathrm{i}{k}_{x}U)\widehat{u}+(U\prime -f)\widehat{v}+\stackrel{\sim }{f}\widehat{w}=-\mathrm{i}{k}_x\widehat{p}+\frac{1}{\mathit{Re}}{\widehat{\mathrm{\nabla }}}^2\widehat{u},\end{array}$$(17)

$$\begin{array}{c}\hfill (\sigma +\mathrm{i}{k}_{x}U)\widehat{v}+f\widehat{u}=-\frac{\mathrm{d}\widehat{p}}{\mathrm{d}y}+\frac{1}{\mathit{Re}}{\widehat{\mathrm{\nabla }}}^2\widehat{v},\end{array}$$(18)

$$\begin{array}{c}\hfill (\sigma +\mathrm{i}{k}_{x}U)\widehat{w}-\stackrel{\sim }{f}\widehat{u}=-\mathrm{i}{k}_z\widehat{p}+N^2\widehat{T}+\frac{1}{\mathit{Re}}{\widehat{\mathrm{\nabla }}}^2\widehat{w},\end{array}$$(19)

$$\begin{array}{c}\hfill (\sigma +\mathrm{i}{k}_{x}U)\widehat{T}-\frac{\stackrel{\sim }{f}U\prime }{N^2}\widehat{v}+\widehat{w}=\frac{1}{\mathit{Pe}}{\widehat{\mathrm{\nabla }}}^2\widehat{T},\end{array}$$(20)

where ${\widehat{\mathrm{\nabla }}}^2=\frac{{\mathrm{d}}^2}{\mathrm{d}{\mathit{y}}^2}-k^2$ with $k^2=k_x^2+k_z^2$. By eliminating $\widehat{p}$ using the continuity Eq. (16), Eqs. (16)–(20) can be gathered into a matrix form as an eigenvalue problem:

$$\begin{array}{c}\hfill \sigma \mathcal{B}\left(\begin{array}{c}\widehat{u}\\ \widehat{v}\\ \widehat{w}\\ \widehat{T}\end{array}\right)=\mathcal{A}\left(\begin{array}{c}\widehat{u}\\ \widehat{v}\\ \widehat{w}\\ \widehat{T}\end{array}\right),\end{array}$$(21)

where 𝒜 and ℬ are the operator matrices expressed as

$$\begin{array}{c}\hfill \mathcal{A}=\left[\begin{array}{cccc}{\mathcal{A}}_{11}& {\mathcal{A}}_{12}& {\mathcal{A}}_{13}& {\mathcal{A}}_{14}\\ {\mathcal{A}}_{21}& {\mathcal{A}}_{22}& {\mathcal{A}}_{23}& {\mathcal{A}}_{24}\\ {\mathcal{A}}_{31}& {\mathcal{A}}_{32}& {\mathcal{A}}_{33}& {\mathcal{A}}_{34}\\ 0& {\mathcal{A}}_{42}& {\mathcal{A}}_{43}& {\mathcal{A}}_{44}\end{array}\right],\end{array}$$(22)

$$\begin{array}{c}\hfill \mathcal{B}=\left[\begin{array}{cccc}-k^2& \mathrm{i}{k}_x\frac{\mathrm{d}}{\mathrm{d}y}& 0& 0\\ 0& {\widehat{\mathrm{\nabla }}}^2& 0& 0\\ 0& \mathrm{i}{k}_z\frac{\mathrm{d}}{\mathrm{d}y}& -k^2& 0\\ 0& 0& 0& 1\end{array}\right],\end{array}$$(23)

where

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathcal{A}}_{11}=k_x(\mathrm{i}{k}^{2}U+k_z\stackrel{\sim }{f})-\frac{k^2}{\mathit{Re}}{\widehat{\mathrm{\nabla }}}^2,\hfill \\ \hfill & {\mathcal{A}}_{12}=k_z^2(U\prime -f)+k_x^{2}U\frac{\mathrm{d}}{\mathrm{d}y}+\frac{\mathrm{i}{k}_x}{\mathit{Re}}{\widehat{\mathrm{\nabla }}}^2\frac{\mathrm{d}}{\mathrm{d}y},\hfill \\ \hfill & {\mathcal{A}}_{13}=k_z^2\stackrel{\sim }{f},{\mathcal{A}}_{14}=k_x{k}_z{N}^2,\hfill \end{array}\end{array}$$(24)

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathcal{A}}_{21}=k^{2}f-\mathrm{i}{k}_z\stackrel{\sim }{f}\frac{\mathrm{d}}{\mathrm{d}y},\hfill \\ \hfill & {\mathcal{A}}_{22}=\mathrm{i}{k}_x(k^{2}U+U\prime -f\frac{\mathrm{d}}{\mathrm{d}y}-U\frac{{\mathrm{d}}^2}{\mathrm{d}{y}^2})+\frac{1}{\mathit{Re}}{\widehat{\mathrm{\nabla }}}^4,\hfill \\ \hfill & {\mathcal{A}}_{23}=\mathrm{i}{k}_x\stackrel{\sim }{f}\frac{\mathrm{d}}{\mathrm{d}y},{\mathcal{A}}_{24}=-\mathrm{i}{k}_z{N}^2\frac{\mathrm{d}}{\mathrm{d}y},\hfill \end{array}\end{array}$$(25)

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathcal{A}}_{31}=-k_x^2\stackrel{\sim }{f},\hfill \\ \hfill & {\mathcal{A}}_{32}=k_x{k}_z(U\frac{\mathrm{d}}{\mathrm{d}y}+f-U\prime )+\frac{\mathrm{i}{k}_z}{\mathit{Re}}{\widehat{\mathrm{\nabla }}}^2,\hfill \\ \hfill & {\mathcal{A}}_{33}=k_x(\mathrm{i}{k}^{2}U-k_z\stackrel{\sim }{f})-\frac{k^2}{\mathit{Re}}{\widehat{\mathrm{\nabla }}}^2,{\mathcal{A}}_{34}=-k_x^2{N}^2,\hfill \end{array}\end{array}$$(26)

$$\begin{array}{c}\hfill {\mathcal{A}}_{42}=\frac{\stackrel{\sim }{f}U\prime }{N^2},{\mathcal{A}}_{43}=-1,{\mathcal{A}}_{44}=-\mathrm{i}{k}_{x}U+\frac{1}{\mathit{Pe}}{\widehat{\mathrm{\nabla }}}^2.\end{array}$$(27)

We discretize numerically the operators 𝒜 and ℬ in the y-direction using the rational Chebyshev function that maps the Chebyshev domain y_cheb ∈ ( − 1, 1) onto the physical space y ∈ ( − ∞, ∞) via the mapping $\mathit{y}/L_{\mathrm{map}}={\mathit{y}}_{\mathrm{cheb}}/\sqrt{1-{\mathit{y}}_{\mathrm{cheb}}^2}$ where L_map is the mapping factor (Deloncle et al. 2007; Park et al. 2020). To distinguish physical and numerically-converged modes from spurious modes, we use a convergence criterion based on the residual of coefficients on the Chebyshev functions proposed by Fabre & Jacquin (2004). This technique allows us to separate highly-oscillatory spurious modes from physical modes that decays smoothly at boundaries as y → ±∞. Furthermore, we chose the number of collocation points in the y-direction from 100 to 200, which was sufficient in our study to confirm the convergence of physical modes. We impose vanishing boundary conditions as y → ±∞ by suppressing terms in the first and last rows of the operator matrices (Antkowiak 2005; Park 2012). Numerical results of the horizontal shear instability are compared and validated with results of Deloncle et al. (2007), Arobone & Sarkar (2012), Park et al. (2020) in stratified and rotating fluids for the traditional case when $\stackrel{\sim }{f}=0$.

While the stability of horizontal shear flows in stratified-rotating fluids in the traditional f-plane has one symmetry condition on σ with ±k_x and ±k_z (Deloncle et al. 2007; Park et al. 2020), an analogous symmetry condition does not exist due to the horizontal Coriolis parameter $\stackrel{\sim }{f}>0$. Instead, there exist two separate symmetry conditions on σ in terms of k_x and k_z as follows:

$$\begin{array}{c}\hfill \sigma (k_x,k_z)={\sigma }^{\ast }(-k_x,-k_z),\end{array}$$(28)

and

$$\begin{array}{c}\hfill \sigma (k_x,-k_z)={\sigma }^{\ast }(-k_x,k_z),\end{array}$$(29)

where the asterisk (*) denotes the complex conjugate. Therefore, in this paper, we investigate both negative and positive k_z for k_x ≥ 0.

2.3. Simplified equations for $\widehat{v}$

Equations (16)–(20) can be simplified into ordinary differential equations (ODEs) for $\widehat{v}$ in particular cases. For instance, in the inviscid and nondiffusive limits (Re → ∞ and Pe → ∞), we can derive the single second-order ODE for $\widehat{v}$ as follows:

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\widehat{v}}{\mathrm{d}{y}^2}+V_1\frac{\mathrm{d}\widehat{v}}{\mathrm{d}y}+V_0\widehat{v}=0,\end{array}$$(30)

where

$$\begin{array}{c}\hfill V_1=\frac{2s{k}_{x}U\prime -2\mathrm{i}{k}_{z}f\stackrel{\sim }{f}}{s^2-N^2-{\stackrel{\sim }{f}}^2}-\frac{Q\prime }{Q},\end{array}$$(31)

s = −iσ + k_xU is the Doppler-shifted frequency, $Q=k^2{s}^2-k_x^2{N}^2$,

$$\begin{array}{c}\hfill \begin{array}{cc}{V}_0=\hfill & -k_z^2\frac{\mathrm{\Gamma }}{s^2-N^2-{\stackrel{\sim }{f}}^2}-k_x^2-\frac{k_x}{s}(U″-\frac{Q\prime }{Q}U\prime )-\frac{k_{x}f}{s}\frac{Q\prime }{Q}\hfill \\ \hfill & -\frac{2k_x^{2}U\prime (U\prime -f)}{s^2-N^2-{\stackrel{\sim }{f}}^2}+\frac{\stackrel{\sim }{f}}{s^2-N^2-{\stackrel{\sim }{f}}^2}[\frac{fQ\prime }{Q}(\mathrm{i}{k}_z-\frac{k_x\stackrel{\sim }{f}}{s})-\stackrel{\sim }{f}{k}_x^2],\hfill \end{array}\end{array}$$(32)

and

$$\begin{array}{c}\hfill \mathrm{\Gamma }=s^2+f(U\prime -f).\end{array}$$(33)

For finite Pe in the inviscid limit, we can find for k_x = 0 the single fourth-order ODE as follows:

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^4\widehat{v}}{\mathrm{d}{y}^4}+{\mathcal{V}}_3\frac{{\mathrm{d}}^3\widehat{v}}{\mathrm{d}{y}^3}+{\mathcal{V}}_2\frac{{\mathrm{d}}^2\widehat{v}}{\mathrm{d}{y}^2}+{\mathcal{V}}_1\frac{\mathrm{d}\widehat{v}}{\mathrm{d}y}+{\mathcal{V}}_0\widehat{v}=0,\end{array}$$(34)

where

$$\begin{array}{c}\hfill {\mathcal{V}}_3=-\frac{\mathrm{i}{k}_z\stackrel{\sim }{f}(U\prime -2f)}{{\sigma }^2+{\stackrel{\sim }{f}}^2},\end{array}$$(35)

$$\begin{array}{c}\hfill {\mathcal{V}}_2=k_z^2(\frac{\mathrm{\Gamma }}{{\sigma }^2+{\stackrel{\sim }{f}}^2}-1)-\frac{3\mathrm{i}{k}_z\stackrel{\sim }{f}U″}{{\sigma }^2+{\stackrel{\sim }{f}}^2}-Pe\frac{\sigma ({\sigma }^2+{\stackrel{\sim }{f}}^2+N^2)}{{\sigma }^2+{\stackrel{\sim }{f}}^2},\end{array}$$(36)

$$\begin{array}{c}\hfill {\mathcal{V}}_1=k_z^3\frac{\mathrm{i}\stackrel{\sim }{f}(U\prime -2f)}{{\sigma }^2+{\stackrel{\sim }{f}}^2}+\frac{2k_z^{2}fU″}{{\sigma }^2+{\stackrel{\sim }{f}}^2}-\frac{3k_z\mathrm{i}\stackrel{\sim }{f}{U}^{″\prime }}{{\sigma }^2+{\stackrel{\sim }{f}}^2}-Pe\frac{2\mathrm{i}{k}_z\stackrel{\sim }{f}f\sigma }{{\sigma }^2+{\stackrel{\sim }{f}}^2},\end{array}$$(37)

$$\begin{array}{c}\hfill {\mathcal{V}}_0=-\frac{k_z^4\mathrm{\Gamma }}{{\sigma }^2+{\stackrel{\sim }{f}}^2}+\frac{k_z^3\mathrm{i}\stackrel{\sim }{f}U″}{{\sigma }^2+{\stackrel{\sim }{f}}^2}+\frac{k_z^{2}fU″}{{\sigma }^2+{\stackrel{\sim }{f}}^2}-\frac{k_z\mathrm{i}\stackrel{\sim }{f}{U}^{″″}}{{\sigma }^2+{\stackrel{\sim }{f}}^2}-Pe\frac{\sigma k_z^2\mathrm{\Gamma }}{{\sigma }^2+{\stackrel{\sim }{f}}^2}.\end{array}$$(38)

In Sect. 4, the above single ODEs (30) and (34) will be used for asymptotic analyses with the WKBJ approximation for large k_z to derive asymptotic dispersion relations for the inertial instability.

3. General stability results

In this section, we investigate examples of numerical stability results to see how the full Coriolis acceleration modifies the instabilities of the horizontal shear flow. Figure 2 shows spectra of the eigenvalue σ for various sets of (k_x, k_z) at Re = ∞,  N = 1, and Pe = 1 for Coriolis components f = 0.5 and $\stackrel{\sim }{f}=2$ (i.e., 2Ω ≃ 2.062 and θ ≃ 76°). In the eigenvalue spectra, the real part of the complex growth rate σ_r determines the stability while the imaginary part σ_i denotes the frequency of the mode. In Fig. 2a for (k_x, k_z) = (0.25, 0), we see different types of eigenvalues: widely distributed stable eigenvalues, clusters of neutral eigenvalues, and one unstable eigenvalue. The neutral or stable modes are continuous modes, which are not exponentially decreasing with |y| → ∞, or gravito-inertial waves that can possess critical points (Astoul et al. 2020). However, this paper only investigates unstable modes that are the most likely to cause turbulence with horizontal shear.

Fig. 2. Eigenvalue spectra in the space (σi, σr) at Re = ∞,  N = 1, Pe = 1, f = 0.5, and $\stackrel{\sim }{f}=2$ (i.e., Ω = 1.031 and θ = 76°) for (a) (kx, kz) = (0.25, 0), (b) (kx, kz) = (0, 6), and (c) (kx, kz) = (0.25, 6). Triangles denote the maximum growth rates.

The eigenfunction of the most unstable mode in Fig. 2a is plotted in Fig. 3 in panels a and b. The mode shape $\widehat{v}$ has maxima around y ≃ ±2.96, and the mode displayed in the physical space (x, y) shows that the velocity v(x, y) at z = 0 has an inclined structure in the direction opposite to the shear. This feature is a characteristic of the inflectional instability mode as previously studied for stratified fluids in the traditional f-plane (Arobone & Sarkar 2012; Park et al. 2020).

Fig. 3. Eigenfunctions of the most unstable modes in Fig. 2 at Re = ∞,  N = 1, Pe = 1, f = 0.5, and $\stackrel{\sim }{f}=2$: (a) the mode shape $\widehat{v}(\mathit{y})$ and (b) velocity v in the physical space (x, y) at z = 0 for (kx, kz) = (0.25, 0), and σ = 0.0712, (c) the mode shape $\widehat{v}(\mathit{y})$ and (d) velocity v in the physical space (y, z) at x = 0 for (kx, kz) = (0, 6), and σ = 0.250. In (a) and (c), black, blue, and red lines denote the absolute value, real part, and imaginary part of the mode shape $\widehat{v}$, respectively.

Figure 2b shows an eigenvalue spectrum that has nearly-neutral clustered eigenvalues, four neutral modes with frequency |σ_i| > 2, a stable eigenvalue, and one unstable eigenvalue. We note that f = 0.5 belongs to the inertial instability regime for k_z >  0 in the traditional approximation, and the instability is not from the inflectional point of the base flow since k_x = 0. For the unstable eigenvalue, we plot the corresponding mode in panels c and d of Fig. 3. We see that the mode shape $\widehat{v}$ is oscillatory in the region |y| < 10, and this wavelike behavior centered around y = 0 reminds us of the inertial instability mode (Park et al. 2020). We also observe clearly in Fig. 3d a characteristic alternating pattern of v when plotted in the physical space (y, z).

In Fig. 2c, we plot the eigenvalue spectra for (k_x, k_z) = (0.25, 6). They have weakly unstable or stable clustered eigenvalues due to the critical point where the Doppler-shifted frequency becomes zero (i.e., s = 0), separate neutral and stable eigenvalues, and one unstable eigenvalue. The mode arisen by the critical point is beyond the scope of this study, and it will be covered by the follow-up work (Astoul et al. 2020). At the wavenumbers k_z = 6 and k_x = 0.25, the inflectional instability is not present and we only have the inertial instability with the mode shape (not shown) similar to that in Fig. 3c.

Figure 4 shows contours of the growth rate σ_r for the most unstable mode in the parameter space (k_x, k_z) for different sets of parameters $(f,\stackrel{\sim }{f},Pe)$ at N = 1 and Re = ∞. We found that the most unstable modes have zero temporal frequency σ_i; therefore, we plot only the real part of the maximum growth rate σ_r. In Fig. 4a, we consider the nondiffusive limit Pe = ∞ in an inertially-stable regime at f = 0 according to the criterion in the traditional approximation (Arobone & Sarkar 2012). If we consider $\stackrel{\sim }{f}=0$ (i.e., 2Ω = 0), only the inflectional instability exists and its maximum growth rate σ_max = 0.1897 is attained at k_x = 0.445 and k_z = 0 (see also, Deloncle et al. 2007). As $\stackrel{\sim }{f}$ increases (i.e., on the equator for rotating stars with 2Ω >  0), we see that the regime of the inflectional instability is widened in the parameter space (k_x, k_z), and the higher maximum growth rate σ_max ≃ 0.277 is attained at k_z = 0 and at a higher value of k_x ≃ 0.54.

Fig. 4. Contours of the maximum growth rate max(σr) in the parameter space of (kx, kz) (solid lines) for Re = ∞ and N = 1 for (a) $(f,\stackrel{\sim }{f},Pe)=(0,2,\infty )$, (b) $(f,\stackrel{\sim }{f},Pe)=(0,2,1)$, and (c) $(f,\stackrel{\sim }{f},Pe)=(0.5,2,1)$. For (a) and (b), this corresponds to (Ω, θ) = (1, 90° ), and (Ω, θ) = (1.031, 76°) for (c). Black and blue solid lines are contours of the maximum growth rate of the inflectional and inertial instabilities, respectively, and gray dashed lines in the background are contours of the maximum growth rate for the same parameters but in the traditional f-plane approximation where $\stackrel{\sim }{f}=0$.

In Fig. 4b, we now consider a finite Péclet number Pe = 1. At $f=\stackrel{\sim }{f}=0$ (i.e., 2Ω = 0), the instability is stabilized as the thermal diffusivity increases (see also, Park et al. 2020) and the unstable regime in the parameter space (k_x, k_z) is slightly shrunk compared to that at Pe = ∞ in Fig. 4a. The maximum growth rate of the inflectional instability at Pe = 1 and $\stackrel{\sim }{f}=0$ still remains as σ_max = 0.1897 at k_x = 0.445. As $\stackrel{\sim }{f}$ increases, we see that the inflectionally-unstable regime at $\stackrel{\sim }{f}=2$ is even more shrunk than the unstable regime at $\stackrel{\sim }{f}=0$. The maximum growth rate is decreased to σ_max ≃ 0.0717 and is found at (k_x, k_z) ≃ (0.26, 0.02). This implies that a slightly three-dimensional inflectional instability is now more unstable than the two-dimensional inflectional instability at finite Pe and $\stackrel{\sim }{f}>0$. What is also very interesting in Fig. 4b is that we observe other unstable regimes. The growth rate contours for large |k_z| are reminiscent of those for the inertial instability, which has a maximum growth rate as |k_z| → ∞ at k_x = 0. It is remarkable to observe this inertial instability at f = 0 (i.e., in the inertially-stable regime in the traditional approximation) when the horizontal Coriolis component $\stackrel{\sim }{f}$ becomes positive and the fluid is thermally diffusive.

We now consider in Fig. 4c the case f = 0.5, which belongs to the inertially-unstable regime in the traditional f-plane. We see a clear difference between growth-rate contours at $\stackrel{\sim }{f}=0$ (i.e., on the pole with 2Ω = 0.5) and those at $\stackrel{\sim }{f}=2$ (i.e., 2Ω ≃ 2.062 and θ ≃ 76°). While the growth rate smoothly increases as |k_z| increases in the traditional case at $\stackrel{\sim }{f}=0$, the regime of the inflectional instability is well separated from the regime of the inertial instability for $\stackrel{\sim }{f}=2$. The inflectional instability at $\stackrel{\sim }{f}=2$ has a maximum growth rate σ_max ≃ 0.0725 around (k_x, k_z) ≃ (0.26, 0.07) while the inertial instability has a maximum growth rate as |k_z| → ∞ at k_x = 0.

In the traditional f-plane approximation, the inflectional instability reaches its maximum growth rate σ_max ≃ 0.1897 in the inviscid limit Re = ∞ at k_x ≃ 0.445 and k_z = 0 (also shown by the gray dashed line in Fig. 5). In this case, the maximum growth rate σ_max is independent from f and N (Park et al. 2020). But without the traditional approximation, the growth rate of the inflectional instability now becomes dependent of N and $\stackrel{\sim }{f}$ as shown in Fig. 5 at Pe = ∞ and k_z = 0. We verified numerically that the most unstable mode is reached for a finite k_x at k_z = 0 in the parameter space of (k_x, k_z) in the inviscid and nondiffusive limits (i.e., Re = Pe = ∞). In these limits, we have the following second-order ODE for $\widehat{v}$

$$\begin{array}{c}\hfill \begin{array}{cc}& \frac{{\mathrm{d}}^2\widehat{v}}{\mathrm{d}{y}^2}+(\frac{2s{k}_{x}U\prime }{s^2-N^2-{\stackrel{\sim }{f}}^2}-\frac{2s{k}_{x}U\prime }{s^2-N^2})\frac{\mathrm{d}\widehat{v}}{\mathrm{d}y}\hfill \\ \hfill & -(k_x^2+\frac{k_{x}U″}{s}+\frac{2k_x^2{U}^{\prime 2}}{s^2-N^2-{\stackrel{\sim }{f}}^2}-\frac{2k_x^2{U}^{\prime 2}}{s^2-N^2})\widehat{v}=0.\hfill \end{array}\end{array}$$(39)

Fig. 5. Inviscid growth rate of the inflectional instability for various parameter sets of $\stackrel{\sim }{f}$ and N at f = 0 (i.e., θ = 90°), Pe = ∞, and kz = 0.

We clearly see that Eq. (39) is still independent from f but it depends on N if $\stackrel{\sim }{f}>0$. For a fixed N, we see in Fig. 5 that the maximum growth rate increases with $\stackrel{\sim }{f}$, while it decreases with N at a fixed $\stackrel{\sim }{f}$. This implies that the stratification stabilizes the inflectional instability, as similarly observed in the traditional approximation (Park et al. 2020). While the wavenumber range of the inflectional instability is 0 <  k_x <  1 at $\stackrel{\sim }{f}=0$, it is noticeable that the instability can sustain for k_x >  1 as $\stackrel{\sim }{f}$ increases.

Figure 6 shows growth rate curves for various values of Pe at N = 1, f = 0, k_x = 0, and $\stackrel{\sim }{f}=2$ (i.e., on the equator with 2Ω = 2). Although f = 0 implies that it is inertially stable in the traditional approximation, it is still unstable for finite Pe at $\stackrel{\sim }{f}=2$ and the growth rate increase as Pe increases. At a fixed Pe, the growth rate increases with k_z and the maximum growth rate is reached as k_z → ∞. Furthermore, the growth rate curves converge to a single curve as Pe → 0.

Fig. 6. Inviscid growth rate of the inertial instability for various Péclet numbers at N = 1 and kx = 0 for f = 0 and $\stackrel{\sim }{f}=2$, i.e., (Ω, θ) = (1, 90° ).

Preliminary results presented in this section tell us that nontraditional effects can significantly modify the properties of the inflectional and inertial instabilities. For instance, the unstable regime of the inertial instability is changed as $\stackrel{\sim }{f}$ increases. The maximum growth rate of the inflectional instability increases with $\stackrel{\sim }{f}$ and depends on N. Similarly to the traditional case at $\stackrel{\sim }{f}=0$, thermal diffusion destabilizes the inertial instability when $\stackrel{\sim }{f}>0$. In the following section, we discuss more thoroughly how the inertial instability depends on physical parameters such as f, $\stackrel{\sim }{f}$, N, or Pe using the WKBJ approximation by taking the asymptotic limit k_z → ∞ for the inviscid case at Re = ∞. We derive analytical expressions of the dispersion relations for the inertial instability in the nondiffusive (Pe = ∞) and highly-diffusive (Pe → 0) cases to understand how the inertial instability is modified as the horizontal Coriolis parameter $\stackrel{\sim }{f}$ increases.

4. Asymptotic description of the inertial instability

Our previous work in Park et al. (2020) with the traditional approximation ($\stackrel{\sim }{f}=0$) was successful to perform detailed analyses on the inertial instability employing the WKBJ approximation in the inviscid limit Re = ∞. For instance, asymptotic dispersion relations were explicitly proposed for large k_z in two limits: Pe → 0 and Pe → ∞. From the dispersion relations, we found that the unstable regime of the inertial instability is 0 <  f <  1 and the maximum growth rate ${\sigma }_{max}=\sqrt{f(1-f)}$ is reached as k_z → ∞. In the traditional approximation, σ_max is related to the epicyclic frequency ω_ep, which satisfies ${\sigma }_{max}^2=-{\omega }_{\mathrm{ep}}^2$ and leads to the Solberg-Høiland criterion for stability when ${\omega }_{\mathrm{ep}}^2+N^2>0$ (Solberg 1936; Høiland 1941). Besides, σ_max is analogous to the growth rate suggested in the limit Pr → 0 for the Goldreich-Schubert-Fricke (GSF) instability (Goldreich & Schubert 1967; Fricke 1968) induced by the horizontal and vertical shear (for more details, we refer the reader to Knobloch & Spruit 1982; Maeder et al. 2013; Barker et al. 2019). Besides, the maximum growth rate σ_max is attained independently of the stratification and thermal diffusivity (i.e., N and Pe) while the first-order term of the growth rate strongly depends on them. However, this argument is not applicable without the traditional approximation as demonstrated by numerical results in the previous section. Thus, it is imperative to investigate the inertial instability in the nontraditional case to see how the properties of the instability are changed as the horizontal Coriolis parameter $\stackrel{\sim }{f}$ increases.

In the following subsections, we perform the WKBJ analysis for large k_z at k_x = 0 in the two limits of the Péclet number Pe: the nondiffusive case as Pe → ∞ and the higly-diffusive case as Pe → 0. This analysis is an extension of the work by Park et al. (2020) in the traditional f-plane approximation. For the two cases, we describe below the properties of the inertial instability such as dispersion relations, the maximum growth rate, or regimes of the instability in the parameter space.

4.1. WKBJ formulation in the nondiffusive limit Pe→∞

In this subsection, the thermally nondiffusive case with Pe → ∞ is considered. We verified numerically that the inviscid maximum growth rate of the inertial instability is found as |k_z| → ∞ at k_x = 0, and it has a zero temporal frequency (i.e., σ_i = 0). For simplicity, we focus on this most unstable case at k_x = 0 and consider only k_z >  0 due to the symmetry σ(0, k_z) = σ(0, −k_z) at k_x = 0. The second-order ODE for $\widehat{v}$ in Eq. (30) at k_x = 0 becomes

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\widehat{v}}{\mathrm{d}{y}^2}+\frac{2\mathrm{i}{k}_z\stackrel{\sim }{f}f}{{\sigma }^2+N^2+{\stackrel{\sim }{f}}^2}\frac{\mathrm{d}\widehat{v}}{\mathrm{d}y}+k_z^2\frac{\mathrm{\Gamma }}{{\sigma }^2+N^2+{\stackrel{\sim }{f}}^2}\widehat{v}=0.\end{array}$$(40)

We note that the term $2\mathrm{i}{k}_z\stackrel{\sim }{f}f/({\sigma }^2+N^2+{\stackrel{\sim }{f}}^2)$ multiplied by the first derivative $\mathrm{d}\widehat{v}/\mathrm{d}\mathit{y}$ has the order of k_z and is comparable with all the other terms when the WKBJ approximation is applied for large k_z. For convenience, we introduce a new variable $\widehat{V}$ prior to the WKBJ analysis

$$\begin{array}{c}\hfill \widehat{V}(y)=\widehat{v}(y)exp\left(\frac{\mathrm{i}{k}_z\stackrel{\sim }{f}f}{{\sigma }^2+N^2+{\stackrel{\sim }{f}}^2}y\right).\end{array}$$(41)

Then Eq. (40) becomes

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\widehat{V}}{\mathrm{d}{y}^2}+\frac{k_z^2}{{\sigma }^2+N^2+{\stackrel{\sim }{f}}^2}[\mathrm{\Gamma }+\frac{{\stackrel{\sim }{f}}^2{f}^2}{{\sigma }^2+N^2+{\stackrel{\sim }{f}}^2}]\widehat{V}=0,\end{array}$$(42)

which is now in the form of the Poincaré equation. Applying to Eq. (42) the WKBJ approximation of $\widehat{V}$ for large k_z:

$$\begin{array}{c}\hfill \widehat{V}(y)\sim exp\left[\frac{1}{\delta }{\displaystyle \sum _{l=0}^{\infty }}{\delta }^l{S}_l(y)\right],\end{array}$$(43)

we get

$$\begin{array}{c}\hfill \delta =\frac{1}{k_z},S_0^{\prime 2}=-\frac{\stackrel{\sim }{\mathrm{\Gamma }}}{{\sigma }^2+N^2+{\stackrel{\sim }{f}}^2},S{\prime }_1=-\frac{S{″}_0}{2S{\prime }_0},\end{array}$$(44)

where

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathrm{\Gamma }}\equiv \mathrm{\Gamma }+\frac{{\stackrel{\sim }{f}}^2{f}^2}{{\sigma }^2+N^2+{\stackrel{\sim }{f}}^2}.\end{array}$$(45)

Since ${\sigma }^2+N^2+{\stackrel{\sim }{f}}^2>0$, the exponential behavior of $\widehat{V}$ determined by the sign of S₀ depends on the sign of $\stackrel{\sim }{\mathrm{\Gamma }}$. On the one hand, the solution is evanescent if $\stackrel{\sim }{\mathrm{\Gamma }}<0$:

$$\begin{array}{c}\hfill \begin{array}{cc}\widehat{V}(y)=\hfill & {(-\stackrel{\sim }{\mathrm{\Gamma }})}^{-\frac{1}{4}}\hfill \\ \hfill & [A_{1}exp\left(\stackrel{\sim }{k}{\int }_y\sqrt{-\stackrel{\sim }{\mathrm{\Gamma }}}\mathrm{d}\upsilon \right)+A_{2}exp(-\stackrel{\sim }{k}{\int }_y\sqrt{-\stackrel{\sim }{\mathrm{\Gamma }}}\mathrm{d}\upsilon )],\hfill \end{array}\end{array}$$(46)

where A₁ and A₂ are constants, and $\stackrel{\sim }{k}=k_z/\sqrt{{\sigma }^2+N^2+{\stackrel{\sim }{f}}^2}$. On the other hand, the solution is wavelike if $\stackrel{\sim }{\mathrm{\Gamma }}>0$:

$$\begin{array}{c}\hfill \widehat{V}(y)={\stackrel{\sim }{\mathrm{\Gamma }}}^{-\frac{1}{4}}[B_{1}exp\left(\mathrm{i}\stackrel{\sim }{k}{\displaystyle {\int }_y}\sqrt{\stackrel{\sim }{\mathrm{\Gamma }}}\mathrm{d}\upsilon \right)+B_{2}exp(-\mathrm{i}\stackrel{\sim }{k}{\int }_y\sqrt{\stackrel{\sim }{\mathrm{\Gamma }}}\mathrm{d}\upsilon )],\end{array}$$(47)

where B₁ and B₂ are constants.

The wavelike or evanescent solution behavior changes at turning points y_t where $\stackrel{\sim }{\mathrm{\Gamma }}({\mathit{y}}_t)=0$. To find y_t, it is useful to define the function σ_t(y) such that

$$\begin{array}{c}\hfill \begin{array}{cc}2{\sigma }_t^2(y)\hfill & =f(U\prime -f)-(N^2+{\stackrel{\sim }{f}}^2)\hfill \\ \hfill & +\sqrt{{(f(U\prime -f)+N^2+{\stackrel{\sim }{f}}^2)}^2+4{\stackrel{\sim }{f}}^2{f}^2},\hfill \end{array}\end{array}$$(48)

(see also, Park & Billant 2012, 2013; Park et al. 2017). By finding σ_t(y) = σ at a given σ, we can find turning points where $\stackrel{\sim }{\mathrm{\Gamma }}(\mathit{y})=0$. An example of σ_t(y) is shown in Fig. 7 for f = 0.5, $\stackrel{\sim }{f}=2$, and N = 1. The gray area denotes where the WKBJ solution is wavelike (i.e., $\stackrel{\sim }{\mathrm{\Gamma }}>0$) while the white area denotes where the solution is evanescent (i.e., $\stackrel{\sim }{\mathrm{\Gamma }}<0$). At a given instability growth rate such that |σ| < σ_t, max, there exist two turning points, one for y >  0 and the other for y <  0. For instance, if the growth rate is σ = 0.250, there exist two turning points: one at y_t+ = 1.365 and the other at y_t− = −1.365. In this case, we can construct an eigenfunction such that the solution $\widehat{V}$ decays exponentially as |y| → ∞ while it is wavelike between the two turning points y_t±. If the growth rate is greater than the maximum of σ_t (i.e., |σ| > |σ_t, max|), the solution $\widehat{V}$ is evanescent everywhere with no turning point and we cannot construct an eigenfunction $\widehat{V}$ satisfying the decaying boundary conditions as |y| → ∞. Therefore, to construct the eigenfunction, the growth rate has to lie in the range 0 <  |σ| < σ_t, max where

$$\begin{array}{c}\hfill 2{\sigma }_{\mathrm{t},max}^2=f(1-f)-N^2-{\stackrel{\sim }{f}}^2+\sqrt{{[f(1-f)+N^2+{\stackrel{\sim }{f}}^2]}^2+4{\stackrel{\sim }{f}}^2{f}^2}.\end{array}$$(49)

Fig. 7. Function σt(y) (solid lines) at N = 1 for $(f,\stackrel{\sim }{f})=(0.5,2)$ (i.e., Ω = 1.031 and θ = 76°). Gray and white areas denote the regimes where $\stackrel{\sim }{\mathrm{\Gamma }}$ is positive and negative, respectively. The dashed line represents an example of the growth rate σ = 0.250.

By performing a turning point analysis in the growth rate range 0 <  σ_r <  σ_t, max, we can further derive an asymptotic dispersion relation as done in Park et al. (2020). We consider first the evanescent WKBJ solution outside the turning point y >  y_t+:

$$\begin{array}{c}\hfill \widehat{V}(y)=A_{\infty }{(-\stackrel{\sim }{\mathrm{\Gamma }})}^{-\frac{1}{4}}exp(-\stackrel{\sim }{k}{\displaystyle {\int }_{y_{t+}}^y}\sqrt{-\stackrel{\sim }{\mathrm{\Gamma }}(\upsilon )}\mathrm{d}\upsilon ),\end{array}$$(50)

where A_∞ is a constant. As y approaches the turning point y_t+, the WKBJ solution (50) is not valid and a local solution is needed to find the WKBJ solution below the turning point y <  y_t+. We use a new scaled coordinate $\stackrel{\sim }{y}=(\mathit{y}-{\mathit{y}}_{t+})/ϵ$ where ϵ is a small parameter defined as $ϵ={\left[{\stackrel{\sim }{k}}^2(-\stackrel{\sim }{\mathrm{\Gamma }}{\prime }_{t+})\right]}^{-1/3}$, $\stackrel{\sim }{\mathrm{\Gamma }}{\prime }_{t+}$ is the derivative of $\stackrel{\sim }{\mathrm{\Gamma }}$ at y_t+, and we assume at leading order that $\stackrel{\sim }{\mathrm{\Gamma }}(\mathit{y})\sim \stackrel{\sim }{\mathrm{\Gamma }}{\prime }_{t+}ϵ\stackrel{\sim }{y}$ around the turning point y_t+. The following local equation is obtained

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\widehat{V}}{\mathrm{d}{\stackrel{\sim }{y}}^2}-\stackrel{\sim }{y}\widehat{V}=O(ϵ).\end{array}$$(51)

Solutions of the local equation (51) are the Airy functions: $\widehat{V}(\stackrel{\sim }{y})=a_1\mathrm{Ai}(\stackrel{\sim }{y})+b_1\mathrm{Bi}(\stackrel{\sim }{y})$, where a₁ and b₁ are constants (Abramowitz & Stegun 1972). From the asymptotic behavior of the Airy functions as $\stackrel{\sim }{\mathit{y}}\to \pm \infty$ and of the WKBJ solution as y → y_t+, we obtain the matched WKBJ solution in the region y_t− <  y <  y_t+:

$$\begin{array}{c}\hfill \widehat{V}(y)={\stackrel{\sim }{\mathrm{\Gamma }}}^{-\frac{1}{4}}[C_{+}exp\left(\mathrm{i}\stackrel{\sim }{k}{\displaystyle {\int }_y^{y_{t+}}}\sqrt{\stackrel{\sim }{\mathrm{\Gamma }}}\mathrm{d}\upsilon \right)+C_{-}exp(-\mathrm{i}\stackrel{\sim }{k}{\displaystyle {\int }_y^{y_{t+}}}\sqrt{\stackrel{\sim }{\mathrm{\Gamma }}}\mathrm{d}\upsilon )],\end{array}$$(52)

where constants C_± satisfy

$$\begin{array}{c}\hfill C_{+}=exp(-\mathrm{i}\frac{\pi }{4})A_{\infty }\mathrm{and}{C}_{-}=exp\left(\mathrm{i}\frac{\pi }{4}\right)A_{\infty }.\end{array}$$(53)

Similarly, we consider the WKBJ solution below the lower turning point y <  y_t− decaying exponentially as y → −∞:

$$\begin{array}{c}\hfill \widehat{V}(y)=A_{-\infty }{(-\stackrel{\sim }{\mathrm{\Gamma }})}^{-\frac{1}{4}}exp(-\stackrel{\sim }{k}{\displaystyle {\int }_y^{y_{t-}}}\sqrt{-\stackrel{\sim }{\mathrm{\Gamma }}(\upsilon )}\mathrm{d}\upsilon ).\end{array}$$(54)

After the local analysis around y_t−, we find the matched wavelike solution in the range y_t− <  y <  y_t+:

$$\begin{array}{c}\hfill \widehat{V}(y)={\stackrel{\sim }{\mathrm{\Gamma }}}^{-\frac{1}{4}}[B_{+}exp\left(\mathrm{i}\stackrel{\sim }{k}{\displaystyle {\int }_{y_{t-}}^y}\sqrt{\stackrel{\sim }{\mathrm{\Gamma }}}\mathrm{d}\upsilon \right)+B_{-}exp(-\mathrm{i}\stackrel{\sim }{k}{\displaystyle {\int }_{y_{t-}}^y}\sqrt{\stackrel{\sim }{\mathrm{\Gamma }}}\mathrm{d}\upsilon )],\end{array}$$(55)

where constants B_± satisfy

$$\begin{array}{c}\hfill B_{+}=exp(-\mathrm{i}\frac{\pi }{4})A_{-\infty }\mathrm{and}{B}_{-}=exp\left(\mathrm{i}\frac{\pi }{4}\right)A_{-\infty }.\end{array}$$(56)

Finally, matching the wavelike solutions (52) and (55) between two turning points leads to the following dispersion relation in the quantized form

$$\begin{array}{c}\hfill \stackrel{\sim }{k}{\displaystyle {\int }_{y_{t-}}^{y_{t+}}}\sqrt{\mathrm{\Gamma }+\frac{{\stackrel{\sim }{f}}^2{f}^2}{{\sigma }^2+N^2+{\stackrel{\sim }{f}}^2}}\mathrm{d}y=(m-\frac{1}{2})\pi ,\end{array}$$(57)

where m is the positive integer denoting the mode number. This quantized dispersion relation becomes equivalent to that in the traditional approximation as $\stackrel{\sim }{f}$ becomes zero (Park et al. 2020).

The right-hand side of Eq. (57) is always fixed at a finite m so that the integral on the left-hand side of Eq. (57) should go to zero as $\stackrel{\sim }{k}\to \infty$. This implies that the two turning points converge to zero as $\stackrel{\sim }{k}\to \infty$. Using this property, we can find a more explicit dispersion relation in the Taylor expansion form for the growth rate σ as

$$\begin{array}{c}\hfill \sigma ={\sigma }_0-\frac{{\sigma }_1}{k_z}+O\left(\frac{1}{k_z^2}\right),\end{array}$$(58)

where

$$\begin{array}{c}\hfill 2{\sigma }_0^2=f(1-f)-N^2-{\stackrel{\sim }{f}}^2+\sqrt{{[f(1-f)+N^2+{\stackrel{\sim }{f}}^2]}^2+4{\stackrel{\sim }{f}}^2{f}^2},\end{array}$$(59)

and

$$\begin{array}{c}\hfill {\sigma }_1=\frac{(m-\frac{1}{2})\sqrt{f}\sqrt{{\sigma }_0^2+N^2+{\stackrel{\sim }{f}}^2}}{{\sigma }_0[1+f^2{\stackrel{\sim }{f}}^2/{({\sigma }_0^2+N^2+{\stackrel{\sim }{f}}^2)}^2]}.\end{array}$$(60)

Since σ₁ is positive, the first mode with m = 1 is the most unstable mode. There are also higher-order modes with m >  1 as investigated in Park et al. (2020), but their growth rate is smaller than that of the first mode. Thus, we hereafter consider only the m = 1 mode for asymptotic results. We note that σ = σ₀ is the maximum growth rate σ_max as k_z → ∞, and σ₀ equals to σ_t, max of Eq. (49). While the maximum growth rate in the traditional f-plane approximation is ${\sigma }_{max}=\sqrt{f(1-f)}$, which depends only on f, we see that the maximum growth rate σ₀ now depends on f, $\stackrel{\sim }{f}$, and N altogether. For instance, at a fixed N, the maximum growth rate σ₀ increases with $\stackrel{\sim }{f}$ as shown in Fig. 8a. Compared to the inertially-unstable regime 0 <  f <  1 in the traditional approximation, we take σ_max real and positive and find the following range of the inertial instability for fixed N and $\stackrel{\sim }{f}$:

$$\begin{array}{c}\hfill 0<f<1+\frac{{\stackrel{\sim }{f}}^2}{N^2}.\end{array}$$(61)

Fig. 8. (a) Contours of the maximum growth rate σ0 from Eq. (59) in the parameter space $(f,\stackrel{\sim }{f})$ for N = 2. Black dashed lines represent the upper limit of the instability $f=1+{\stackrel{\sim }{f}}^2/N^2$ from Eq. (61) for different N, and the white dashed line represents fmax from (64) where the maximum growth rate is attained. (b) Contours of σ0 from Eq. (59) in the parameter space $(N,\stackrel{\sim }{f})$ at f = 1. (c) Contours of the first-order term σ1 in the parameter space $(f,\stackrel{\sim }{f})$ for N = 2 and the first branch at m = 1.

If we use the relations f = 2Ω cos θ and $\stackrel{\sim }{f}=2\mathrm{\Omega }sin\theta$, we have the corresponding range

$$\begin{array}{c}\hfill 0<2\mathrm{\Omega }cos\theta <1+\frac{4{\mathrm{\Omega }}^2{sin}^2\theta }{N^2},\end{array}$$(62)

which is equivalent to

$$\begin{array}{c}\hfill {tan}^{-1}\left(\frac{N}{2}\right)<\theta <{90}^{°},\end{array}$$(63)

if 2Ω >  0.

Moreover, at the vertical Coriolis parameter f = f_max

$$\begin{array}{c}\hfill f_{max}=\frac{1}{4}-N^2+\sqrt{{(\frac{1}{4}-N^2)}^2+N^2+{\stackrel{\sim }{f}}^2},\end{array}$$(64)

we can find the maximum growth rate at fixed $\stackrel{\sim }{f}$ and N as

$$\begin{array}{c}\hfill max\left({\sigma }_{max}\right){|}_{\stackrel{\sim }{f},N}={\sigma }_0(f_{max}).\end{array}$$(65)

In Fig. 8b, we see that the maximum growth rate σ₀ at fixed f and $\stackrel{\sim }{f}$ decreases as N increases, which highlights the stabilizing role of the stratification. We also display in Fig. 8c the dependence of the first-order term σ₁ on f and $\stackrel{\sim }{f}$ at N = 2. In this case, σ₁ increases with f at a fixed $\stackrel{\sim }{f}$ and goes to infinity as f reaches its upper limit $f=1+{\stackrel{\sim }{f}}^2/N^2$ since σ₀ is in the denominator of σ₁ and σ₀ → 0 as $f\to 1+{\stackrel{\sim }{f}}^2/N^2$.

Figure 9 shows the growth rate σ versus the vertical wavenumber k_z for different values of $\stackrel{\sim }{f}$ and N at f = 1 and k_x = 0 in the inviscid and nondiffusive limits. The vertical Coriolis parameter f = 1 belongs to the inertially stable regime under the traditional approximation, but we see that it becomes unstable if we consider the nontraditional effects with the horizontal component $\stackrel{\sim }{f}>0$, and the growth rate increases with $\stackrel{\sim }{f}$, as predicted from the WKBJ analysis. We notice that numerical stability results are in good agreement with asymptotic predictions from Eq. (58), especially for large k_z. We also verify that the stratification stabilizes the inertial instability, and the onset of the inertial instability appears at a higher k_z as N increases.

Fig. 9. Numerical results (lines) and asymptotic results (circles) from (58) for various values of $\stackrel{\sim }{f}$ and N at f = 1, kx = 0, Re = ∞, and Pe = ∞.

4.2. The WKBJ formulation in the limit Pe→0

In this subsection, we now turn our attention to the highly-diffusive case in the limit of Pe → 0. This limit is relevant for stellar interiors. To analyze the inertial instability, we consider the case with k_x = 0 and use the fourth-order ODE (34) for the WKBJ analysis. As performed in the previous subsection, we introduce for convenience the variable $\widehat{W}$

$$\begin{array}{c}\hfill \widehat{W}(y)=\widehat{v}(y)exp[-\frac{\mathrm{i}{k}_z\stackrel{\sim }{f}(U-2fy)}{2({\sigma }^2+{\stackrel{\sim }{f}}^2)}],\end{array}$$(66)

to understand more clearly the exponential or wavelike behavior around turning points. The fourth-order ODE (34) can be rewritten in terms of $\widehat{W}$ as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}& \frac{1}{k_z^4}\frac{{\mathrm{d}}^4\widehat{W}}{\mathrm{d}{y}^4}+\frac{\mathrm{i}\stackrel{\sim }{f}(U\prime -2f)}{k_z^3({\sigma }^2+{\stackrel{\sim }{f}}^2)}\frac{{\mathrm{d}}^3\widehat{W}}{\mathrm{d}{y}^3}+\frac{1}{k_z^2}(\frac{\mathrm{\Gamma }}{{\sigma }^2+{\stackrel{\sim }{f}}^2}-1)\frac{{\mathrm{d}}^2\widehat{W}}{\mathrm{d}{y}^2}\hfill \\ \hfill & +[\frac{\mathrm{i}\stackrel{\sim }{f}(U\prime -2f)}{k_z({\sigma }^2+{\stackrel{\sim }{f}}^2)}(\frac{\mathrm{\Gamma }}{{\sigma }^2+{\stackrel{\sim }{f}}^2}+\frac{{\stackrel{\sim }{f}}^2{(U\prime -2f)}^2}{4{({\sigma }^2+{\stackrel{\sim }{f}}^2)}^2})+\frac{U″}{k_z^2({\sigma }^2+{\stackrel{\sim }{f}}^2)}\hfill \\ \hfill & \times (2f+\frac{3(U\prime -2f)}{2({\sigma }^2+{\stackrel{\sim }{f}}^2)})]\frac{\mathrm{d}\widehat{W}}{\mathrm{d}y}+\{[\frac{-\mathrm{\Gamma }}{{\sigma }^2+{\stackrel{\sim }{f}}^2}-\frac{{\stackrel{\sim }{f}}^2{(U\prime -2f)}^2}{4{({\sigma }^2+{\stackrel{\sim }{f}}^2)}^2}]\hfill \\ \hfill & \times [1+\frac{{\stackrel{\sim }{f}}^2{(U\prime -2f)}^2}{4{({\sigma }^2+{\stackrel{\sim }{f}}^2)}^2}]+\frac{1}{k_z}[\frac{\mathrm{i}\stackrel{\sim }{f}U″}{2({\sigma }^2+{\stackrel{\sim }{f}}^2)}(1+\frac{\mathrm{\Gamma }}{{\sigma }^2+{\stackrel{\sim }{f}}^2})\hfill \\ \hfill & +\frac{\mathrm{i}\stackrel{\sim }{f}U″(U\prime -2f)}{{({\sigma }^2+{\stackrel{\sim }{f}}^2)}^2}(f+\frac{3{\stackrel{\sim }{f}}^2(U\prime -2f)}{4({\sigma }^2+{\stackrel{\sim }{f}}^2)})]\}\widehat{W}=O\left(\frac{1}{k_z^2}\right),\hfill \end{array}\end{array}$$(67)