© J. Park et al. 2020

1. Introduction

The combination of space-based helio- and asteroseismology has demonstrated that stably stratified rotating stellar radiation zones are the seats of efficient transport of angular momentum throughout the evolution of stars. This strong transport leads to a uniform rotation in the case of the Sun down to 0.2 R_⊙ (García et al. 2007) and to weak differential rotation in other stars (e.g., Mosser et al. 2012; Deheuvels et al. 2012, 2014; Kurtz et al. 2014; Saio et al. 2015; Murphy et al. 2016; Spada et al. 2016; Van Reeth et al. 2016, 2018; Aerts et al. 2017; Gehan et al. 2018; Ouazzani et al. 2019). Four main mechanisms that transport angular momentum and mix chemicals are present in stellar radiation zones (e.g., Maeder 2009; Mathis 2013; Aerts et al. 2019, and references therein): instabilities of the differential rotation (e.g., Zahn 1983, 1992), stable and unstable magnetic fields (e.g., Spruit 1999; Fuller et al. 2019), internal gravity waves (e.g., Zahn et al. 1997; Talon & Charbonnel 2005; Pinçon et al. 2017), and large-scale meridional circulations (e.g., Zahn 1992; Maeder & Zahn 1998; Mathis & Zahn 2004). Important progress has been made in their modeling during the last two decades. However, many simplifying assumptions are still employed in the treatment of these four mechanisms because of their complexity and of the broad range of spatial and temporal scales they involve. For instance, the complex interplay between rotation and stratification for the study of vertical and horizontal shear instabilities is partially treated. On the one hand, the action of the Coriolis acceleration is not taken into account in the state-of-the-art modeling of the vertical turbulent transport due to the instabilities by radial differential rotation. On the other hand, the important action of thermal diffusion has completely been ignored in the studies of the horizontal turbulent transport induced by horizontal and vertical shear instabilities (see e.g., Mathis et al. 2004, 2018, respectively). As a consequence, a lot of work is still required to obtain robust abinitio prescriptions for transport processes in stars and to reconcile stellar models and the observations (e.g., Eggenberger et al. 2012, 2019; Ceillier et al. 2013; Marques et al. 2013; Cantiello et al. 2014). In this work, our aim is to improve our understanding of horizontal shear instabilities in rotating stellar radiation zones.

In his seminal article, Zahn (1992) built the theoretical framework to study the turbulent transport induced by vertical and horizontal shears in convectively stable rotating stellar radiation zones. In such regions, turbulence can be anisotropic because of the combined action of buoyancy and the Coriolis acceleration, which control the turbulent motions along the vertical and the horizontal directions, respectively (e.g., Billant & Chomaz 2001; Davidson 2013; Mathis et al. 2018). As pointed out in Mathis et al. (2018), many theoretical works have been devoted to providing robust prescriptions for the vertical eddy diffusion associated with the instability of vertical shear (e.g., Zahn 1992; Maeder 1995, 1997; Maeder & Meynet 1996; Talon & Zahn 1997; Kulenthirarajah & Garaud 2018). Their predictions are now being tested using direct numerical simulations (Prat & Lignières 2013, 2014; Prat et al. 2016; Garaud et al. 2017; Gagnier & Garaud 2018). However, very few studies have examined the instabilities arising from horizontal shear (e.g., Zahn 1992; Maeder 2003; Mathis et al. 2004). The resulting prescriptions are based mostly on phenomenological approaches and the results of unstratified Taylor–Couette flow experiments (Richard & Zahn 1999) that do not take into account the complex interplay of stratification, rotation, and thermal diffusion, all of which play important roles in stellar radiation zones. For instance, heat diffusion inhibits the effects of the stable stratification for the vertical shear instability (Townsend 1958; Zahn 1983). Mathis et al. (2018) studies the combined action of stratification and rotation on the horizontal turbulent transport. Their approaches to the turbulent transport have been successfully implemented in several stellar evolution codes (Talon et al. 1997; Meynet & Maeder 2000; Palacios et al. 2006; Marques et al. 2013). However, the neglect of the nonadiabatic aspects of the problem and the focus only on the latitudinal turbulent transport induced by the three-dimensional motions resulting from the vertical shear instability constitute weak points of the theory for the rotational mixing. Indeed, in the formalism of Mathis et al. (2018), a stronger turbulent transport is assumed along the horizontal direction than along the vertical direction. This allows us to assume that the horizontal gradients of the angular velocity (and of the fluctuations of temperature and chemical composition) are smoothed out, leading to the so-called “shellular” rotation, that only varies with the radius. Such a radial variation of rotation is pertinent to the vertical shear instability, which has been mostly studied on the local f-plane with vertical stratification (see e.g., Wang et al. 2014). The horizontal shear instability with stellar differential rotation in latitudinal direction has also been studied (see e.g.,Watson 1980; Garaud 2001; Kitchatinov & Rüdiger 2009), but the combined impacts of the rotation, stratification, and thermal diffusion on horizontal shear instability are not yet fully resolved. It is thus mandatory to study the instabilities of horizontal shear in rotating and stably stratified stellar regions (Figs. 1a and b) while taking into account their important thermal diffusion.

Fig. 1. Panels a and b: horizontal shear flow U(y) on a local polar plane rotating with angular speed Ω in the radiative zone of a rotating star. The radiative and convective zones, colored as yellow and orange with the transition altitude zc, are the configuration of low-mass stars; it should be inverted for early-type stars. Panel c: horizontal shear flow profile U(y) with an inflection point U″ = 0 for the inflectional instability. Panel d: Rayleigh discriminant Φ(y) = f0(f0−U′) for the inertial instability configuration.

In classical fluid dynamics, theories of shear flow stability in homogeneous fluids are well known (Schmid & Henningson 2001). According to Rayleigh’s inflection point criterion, an inviscid unstable shear flow should possess an inflectional point where the second derivative of the zonal flow velocity U vanishes (i.e., U″(y) = 0, where y is the local latitudinal coordinate and prime denotes the derivative with respect to y; see Fig. 1c). For instance, a shear flow with a hyperbolic tangent velocity profile can undergo this inflection point instability (Michalke 1964). In geophysical fluid dynamics, this instability is also referred to the barotropic instability (Kundu & Cohen 2001) as it can be triggered by a two-dimensional horizontal perturbation with a finite streamwise zonal wavenumber k_x. However, general three-dimensional perturbations with both k_x and a vertical wavenumber k_z can also trigger the inflectional instability of the horizontal shear flow. So, we use the terminology inflectional instability hereafter. The stability analysis of such shear flow in stratified fluids by Deloncle et al. (2007) revealed that the unstable regime in the parameter space of zonal and vertical wavenumbers (k_x, k_z) is widely broadened for a horizontal shear flow U(y) in strongly stratified fluids. This implies that perturbations with a smaller vertical scale (i.e., large k_z) can still be unstable with a large Brunt–Väisälä frequency N. A selfsimilarity of the three-dimensional inflectional instability is found for strong stratification with a rescaled parameter Nk_z (Deloncle et al. 2007). However, they found that the most unstable perturbation is still two-dimensional with a finite k_x at k_z = 0 independent of the stratification.

The effects of both the stratification and rotation on horizontal shear instability were studied by Arobone & Sarkar (2012). They explored how the instability growth rate in the parameter space (k_x, k_z) is modified as the stratification and rotation change. And they found that the most unstable mode of the inflectional instability is always found at finite k_x and k_z = 0 independent of the stratification and rotation. Moreover, they also underlined that there exists an inertial instability in a finite range of f₀ which is the Coriolis parameter defined as f₀ = 2Ω cos θ_s where Ω is the rotation of a star and θ_s is the colatitude. The origin of the inertial instability is different from the inflectional instability: the horizontal flow can become inertially unstable only in rotating fluids if the Rayleigh discriminant Φ(y) = f₀(f₀ − U′) becomes negative (Fig. 1d). This condition leads to the inertially-unstable range 0 <  f₀ <  max(U′) for the hyperbolic tangent shear flow whose the maximum growth rate is found in inviscid limit as $\sqrt{f_0(max(U\prime )-f_0)}$. Such a mechanism is essentially equivalent to that of the centrifugal instability for rotating flows with cylindrical geometry. And it is centrifugally stable if the cylindrical Rayleigh discriminant Φ_r is always positive:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{r}}=\frac{1}{r^3}\frac{\partial {(r{u}_{\theta })}^2}{\partial r}>0,\end{array}$$(1)

where u_θ is the azimuthal velocity and r is the cylindrical radial coordinate (Kloosterziel & van Heijst 1991; Park & Billant 2013). The stability condition can be extended by taking the stratification of the entropy S into account. This results in the Solberg–Høiland conditions:

$$\begin{array}{cc}& {\mathrm{\Phi }}_{\mathrm{r}}-\frac{1}{C_{\mathrm{p}}\rho }\mathrm{\nabla }p·\mathrm{\nabla }S>0,\hfill \\ \hfill & \frac{\partial p}{\partial z}(\frac{\partial {(r{u}_{\theta })}^2}{\partial r}\frac{\partial S}{\partial z}-\frac{\partial S}{\partial r}\frac{\partial {(r{u}_{\theta })}^2}{\partial z})<0,\hfill \end{array}$$(2)

(see also, Rüdiger et al. 2002). Moreover, the Goldreich–Schubert–Fricke (GSF) criterion proposes the stability condition which takes effects of the thermal diffusivity κ₀ and viscosity ν₀ into account:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{r}}+\frac{{\nu }_0}{{\kappa }_0}{N}^2>0,\end{array}$$(3)

where N is the Brunt–Väisälä frequency (Goldreich & Schubert 1967; Fricke 1968; Maeder et al. 2013). While these studies have proposed stability conditions in the presence of the stratification, rotation, and thermal diffusion, it is still not fully understood how horizontal shear instabilities are modified by the thermal diffusion, which has an essential importance for the dynamics of stellar radiative zones.

In this paper, we investigate the linear stability of horizontal shear flows in stably-stratified, rotating, and thermally-diffusive fluids in the context of stellar radiation zones. In Sect. 2, equations for the linear stability analysis are formulated. In Sect. 3, we compute numerical results of this analysis for both the inflectional and inertial instabilities to characterize their main properties. In Sect. 4, these are compared with results from asymptotic analyses by means of the Wentzel–Kramers–Brillouin–Jeffreys (WKBJ) approximation for the inertial instability in order to understand the important role of thermal diffusion in stably-stratified rotating fluids. In Sect. 5, key scaling laws are derived for both the inflectional and inertial instabilities as a function of the stratification and thermal diffusivity. Finally, conclusions and perspectives of this work are presented in Sect. 6.

2. Problem formulation

2.1. Governing equations and base equilibrium state

We consider the Euler equations under the Boussinesq approximation and the heat transport equation in the Cartesian coordinate (x, y, z) in a local frame rotating with angular velocity Ω:

$$\begin{array}{cc}& \mathrm{\nabla }·\mathit{u}=0,\hfill \end{array}$$(4)

$$\begin{array}{cc}& \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}}\theta \mathit{g},\hfill \end{array}$$(5)

$$\begin{array}{cc}& \frac{\partial \theta }{\partial t}+\mathit{u}·\mathrm{\nabla }\theta ={\kappa }_0{\mathrm{\nabla }}^2\theta ,\hfill \end{array}$$(6)

where u = (u,v,w) is the velocity, p is the pressure, θ = T − T₀ is the temperature deviation from the reference temperature T₀, f = (0,0,2Ω) is the Coriolis vector which is here antialigned with the gravity g = (0, 0, −g), ρ₀ is the reference density, κ₀ is the reference thermal diffusivity, α_T is the thermal expansion coefficient, and ∇² denotes the Laplacian operator. We thus consider as a first step the traditional f-plane case corresponding to the horizontal shear flow on the pole (Fig. 1) to compare with previous literature. In this setup, the action of the latitudinal component of the rotation vector is filtered out (Gerkema & Shrira 2005; Gerkema et al. 2008). The current Cartesian coordinates (x, y, z) on the traditional f-plane is associated with the spherical coordinates (r, θ_s, φ) where r is the radial coordinate and φ is the longitude. For instance, x is the longitudinal coordinate with e_x = e_φ where e denotes the unit vector of the coordinate systems, y is the latitudinal coordinate with e_y = −e_θs, and z is the vertical coordinate with e_z = e_r.

To perform the linear stability analysis, we consider a steady base velocity U = (U(y),0,0) in a hyperbolic tangent form

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

where U₀ and L₀ are the reference velocity and length scale, respectively. Such a base flow is balanced with the pressure gradient as

$$\begin{array}{c}\hfill fU=-\frac{\partial \overline{P}}{\partial y},\end{array}$$(8)

where $\overline{P}$ is the base pressure profile. We consider a base temperature profile $\overline{\mathrm{\Theta }}(z)$ which increases linearly with height z:

$$\begin{array}{c}\hfill \overline{\mathrm{\Theta }}=\frac{\mathrm{\Delta }{\mathrm{\Theta }}_0}{\mathrm{\Delta }z}z,\end{array}$$(9)

where ΔΘ₀ is the base-temperature difference along the vertical distance Δz.

2.2. Linearized stability equations

Subject to the base state, we consider perturbations $\stackrel{\sim }{\mathit{u}}=\mathit{u}-\mathit{U}=(\stackrel{\sim }{u},\stackrel{\sim }{\mathit{v}},\stackrel{\sim }{\mathit{w}})$, $\stackrel{\sim }{p}=p-\overline{P}$, and $\stackrel{\sim }{T}=\theta -\overline{\mathrm{\Theta }}$. Hereafter, we use dimensionless parameters by converting the Eqs. (4)–(6) into a set of dimensionless equations with the length scale as L₀, the velocity scale as U₀, the time scale as L₀/U₀, the pressure scale as ${\rho }_0{U}_0^2$, and the temperature scale as (L₀ΔΘ₀)/Δz. For infinitesimally small amplitude perturbations, we obtain the following linearized perturbation equations

$$\begin{array}{cc}& \frac{\partial \stackrel{\sim }{u}}{\partial x}+\frac{\partial \stackrel{\sim }{v}}{\partial y}+\frac{\partial \stackrel{\sim }{w}}{\partial z}=0,\hfill \end{array}$$(10)

$$\begin{array}{cc}& \frac{\partial \stackrel{\sim }{u}}{\partial t}+U\frac{\partial \stackrel{\sim }{u}}{\partial x}+(U^{\prime }-f)\stackrel{\sim }{v}=-\frac{\partial \stackrel{\sim }{p}}{\partial x},\hfill \end{array}$$(11)

$$\begin{array}{cc}& \frac{\partial \stackrel{\sim }{v}}{\partial t}+U\frac{\partial \stackrel{\sim }{v}}{\partial x}+f\stackrel{\sim }{u}=-\frac{\partial \stackrel{\sim }{p}}{\partial y},\hfill \end{array}$$(12)

$$\begin{array}{cc}& \frac{\partial \stackrel{\sim }{w}}{\partial t}+U\frac{\partial \stackrel{\sim }{w}}{\partial x}=-\frac{\partial \stackrel{\sim }{p}}{\partial z}+N^2\stackrel{\sim }{T},\hfill \end{array}$$(13)

$$\begin{array}{cc}& \frac{\partial \stackrel{\sim }{T}}{\partial t}+U\frac{\partial \stackrel{\sim }{T}}{\partial x}+\stackrel{\sim }{w}=\frac{1}{\mathrm{Pe}}{\mathrm{\nabla }}^2\stackrel{\sim }{T},\hfill \end{array}$$(14)

where prime denotes the total derivative with respect to y, N is the dimensionless Brunt–Väisälä frequency where

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

f is the dimensionless Coriolis parameter

$$\begin{array}{c}\hfill f=\frac{2\mathrm{\Omega }{L}_0}{U_0},\end{array}$$(16)

and Pe is the Péclet number

$$\begin{array}{c}\hfill \mathrm{Pe}=\frac{U_0{L}_0}{{\kappa }_0}·\end{array}$$(17)

We note that N is equivalent to the inverse of the horizontal Froude number F_h for the stratified horizontal shear flow in Deloncle et al. (2007). Also, N² is similar to the Richardson number Ri defined as $Ri=N_{\mathit{v}}^2/S_{\mathit{v}}^2$ where S_v is the vertical shear and N_v is the Brunt–Väisälä frequency used for the vertical shear instability study (Lignières et al. 1999). The nondimensional f is equal to the inverse of the Rossby number Ro = 1/f.

To perform a linear stability analysis, we consider a normal mode representation for the perturbation variables:

$$\begin{array}{c}\hfill (\stackrel{\sim }{\mathit{u}},\stackrel{\sim }{p},\stackrel{\sim }{T})=\mathfrak{R}[(\widehat{\mathit{u}}(y),\widehat{p}(y),\widehat{T}(y))exp[i(k_{x}x+k_{z}z)+\sigma t]],\end{array}$$(18)

where i² = −1, $\widehat{\mathit{u}}=(\widehat{u},\widehat{\mathit{v}},\widehat{\mathit{w}})$, $\widehat{p}$, and $\widehat{T}$ are the latitudinal mode shapes for velocity, pressure, and temperature perturbation, respectively, k_x is the horizontal wavenumber in streamwise direction, k_z is the vertical wavenumber, and σ = σ_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. Expanding the Eqs. (10)–(14) with these normal modes, we obtain the following linear stability equations

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

$$\begin{array}{cc}& (\sigma +i{k}_{x}U)\widehat{u}+(U^{\prime }-f)\widehat{v}=-i{k}_x\widehat{p},\hfill \end{array}$$(20)

$$\begin{array}{cc}& (\sigma +i{k}_{x}U)\widehat{v}+f\widehat{u}=-\frac{\mathrm{d}\widehat{p}}{\mathrm{d}y},\hfill \end{array}$$(21)

$$\begin{array}{cc}& (\sigma +i{k}_{x}U)\widehat{w}=-i{k}_z\widehat{p}+N^2\widehat{T},\hfill \end{array}$$(22)

$$\begin{array}{cc}& (\sigma +i{k}_{x}U)\widehat{T}+\widehat{w}=\frac{1}{\mathrm{Pe}}{\widehat{\mathrm{\nabla }}}^2\widehat{T},\hfill \end{array}$$(23)

where ${\widehat{\mathrm{\nabla }}}^2={\mathrm{d}}^2/\mathrm{d}{\mathit{y}}^2-k^2$ with $k^2=k_x^2+k_z^2$. Due to the symmetry σ(k_x, k_z) = σ(k_x, −k_z) = σ(−k_x, k_z)^* = σ(−k_x, −k_z)^* where the asterisk denotes the complex conjugate; we thus consider only the positive wavenumbers k_x and k_z in this paper. For convenience and mathematical simplicity, we can simplify the set of Eqs. (19)–(23) into a single ordinary differential equation (ODE) for $\widehat{T}$. Firstly, from the vertical momentum and temperature equations, we express $\widehat{\mathit{w}}$ and $\widehat{p}$ as a function of $\widehat{T}$:

$$\begin{array}{cc}& \widehat{w}=[-(\sigma +i{k}_{x}U)+\frac{1}{\mathrm{Pe}}{\widehat{\mathrm{\nabla }}}^2]\widehat{T},\hfill \end{array}$$(24)

$$\begin{array}{cc}& \widehat{p}=\frac{i}{k_z}[-{(\sigma +i{k}_{x}U)}^2+\frac{\sigma +i{k}_{x}U}{\mathrm{Pe}}{\widehat{\mathrm{\nabla }}}^2-N^2]\widehat{T}.\hfill \end{array}$$(25)

From the horizontal momentum Eqs. (20) and (21), we can express $\widehat{u}$ and $\widehat{\mathit{v}}$ in terms of $\widehat{p}$:

$$\begin{array}{cc}& \widehat{u}=\frac{(f-U^{\prime })\frac{\mathrm{d}\widehat{p}}{\mathrm{d}y}+i{k}_x(\sigma +i{k}_{x}U)\widehat{p}}{f(U^{\prime }-f)-{(\sigma +i{k}_{x}U)}^2},\hfill \end{array}$$(26)

$$\begin{array}{cc}& \widehat{v}=\frac{(\sigma +i{k}_{x}U)\frac{\mathrm{d}\widehat{p}}{\mathrm{d}y}-i{k}_{x}f\widehat{p}}{f(U^{\prime }-f)-{(\sigma +i{k}_{x}U)}^2}·\hfill \end{array}$$(27)

Applying $\widehat{p}$ of (25) to $\widehat{u}$ and $\widehat{\mathit{v}}$, we have the velocity perturbation $\widehat{u}=(\widehat{u},\widehat{\mathit{v}},\widehat{\mathit{w}})$ expressed by $\widehat{T}$ and the continuity Eq. (19) becomes the single 4th-order ODE for $\widehat{T}$:

$$\begin{array}{cc}& \frac{{\mathrm{d}}^2\widehat{T}}{\mathrm{d}{y}^2}+(\frac{4s{s}^{\prime }}{s^2-N^2}-\frac{\mathrm{\Gamma }{}^{\prime }}{\mathrm{\Gamma }})\frac{\mathrm{d}\widehat{T}}{\mathrm{d}y}\hfill \\ \hfill & +[k_z^2\frac{\mathrm{\Gamma }}{N^2-s^2}-k_x^2+f\frac{k_x\mathrm{\Gamma }{}^{\prime }}{s\mathrm{\Gamma }}+\frac{2s{s}^{\prime }}{s^2-N^2}(\frac{U^{″}}{U^{\prime }}+\frac{s^{\prime }}{s}-\frac{\mathrm{\Gamma }{}^{\prime }}{\mathrm{\Gamma }})]\widehat{T}\hfill \\ \hfill & +\frac{1}{\mathrm{Pe}}\left(\frac{\mathit{is}}{s^2-N^2}\right)[{\widehat{\mathrm{\nabla }}}^4\widehat{T}+(\frac{s^{\prime }}{s}-\frac{\mathrm{\Gamma }{}^{\prime }}{\mathrm{\Gamma }}){\widehat{\mathrm{\nabla }}}^2\frac{\mathrm{d}\widehat{T}}{\mathrm{d}y}\hfill \\ \hfill & +\{\frac{s^{\prime }}{s}(\frac{U^{″}}{U^{\prime }}-\frac{\mathrm{\Gamma }{}^{\prime }}{\mathrm{\Gamma }})+f(\frac{k_x\mathrm{\Gamma }{}^{\prime }}{s\mathrm{\Gamma }}-\frac{k_z^2}{s^2}(U^{\prime }-f))\}{\widehat{\mathrm{\nabla }}}^2\widehat{T}]=0,\hfill \end{array}$$(28)

where ${\widehat{\mathrm{\nabla }}}^4={({\widehat{\mathrm{\nabla }}}^2)}^2$, s = −iσ + k_xU is the complex Doppler-shifted frequency, and Γ is the function defined as

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

We note that Eq. (28) becomes the 2nd-order ODE in the nondiffusive limit Pe → ∞ while it becomes independent of N for the high-diffusivity case as Pe → 0.

2.3. Numerical method

We solve the set of Eqs. (19)–(23) numerically by considering an eigenvalue problem in a simplified matrix form:

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

where 𝒜 and ℬ are operator matrices expressed as

$$\begin{array}{cc}& \mathcal{A}=\left[\begin{array}{ccc}{\mathcal{A}}_{11}& {\mathcal{A}}_{12}& 0\\ i{k}^{2}f& {\mathcal{A}}_{22}& N^2{k}_z\frac{\mathrm{d}}{\mathrm{d}y}\\ i{k}_x& \frac{\mathrm{d}}{\mathrm{d}y}& k_x{k}_{z}U+\frac{i{k}_z}{\mathrm{Pe}}{\widehat{\mathrm{\nabla }}}^2\end{array}\right],\hfill \end{array}$$(31)

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

where

$$\begin{array}{cc}& {\mathcal{A}}_{11}=i{k}_x(U^{\prime }+U\frac{\mathrm{d}}{\mathrm{d}y})-i{k}_{x}f,\hfill \\ \hfill & {\mathcal{A}}_{12}=k_x^{2}U+(U^{\prime }-f)\frac{\mathrm{d}}{\mathrm{d}y}+U^{″},\hfill \\ \hfill & {\mathcal{A}}_{22}=k_{x}U{\widehat{\mathrm{\nabla }}}^2+f{k}_x\frac{\mathrm{d}}{\mathrm{d}y}-k_x{U}^{″}.\hfill \end{array}$$(33)

The eigenvalue problem (30) is discretized in the y-direction using the rational Chebyshev functions with the mapping $\stackrel{\sim }{\mathit{y}}=\mathit{y}/\sqrt{1+{\mathit{y}}^2}$ projecting the Chebyshev domain $\stackrel{\sim }{\mathit{y}}\in (-1,1)$ onto the physical space y ∈ ( − ∞, ∞) (Deloncle et al. 2007; Park 2012). Vanishing boundary conditions are imposed for y → ±∞ by suppressing terms in the first and last rows of the operator matrices (Antkowiak 2005). We note that the eigenvalue problem (30) is identical to that of Arobone & Sarkar (2012) in the inviscid limit and when the density perturbation analogous to $-\widehat{T}$ is considered. The numerical results are validated with those from Deloncle et al. (2007) and Arobone & Sarkar (2012) in stratified and rotating fluids.

3. General results

Horizontal shear flows in stably stratified and rotating fluids are prone to two types of destabilizing mechanisms: the inflectional and inertial instabilities. The inflectional instability occurs when there exists an inflection point y_i where U″(y_i) = 0. On the other hand, the inertial instability can occur without an inflection point when there is an imbalance between the pressure gradient and the inertial force, the mechanism equivalent for the centrifugal instability in cylindrical geometry (Kloosterziel & van Heijst 1991; Billant & Gallaire 2005; Wang et al. 2014). For the hyperbolic tangent shear flow (7) in stratified-rotating fluids, the inflectional instability is always present while the inertial instability exists only in the range 0 <  f <  1 (Arobone & Sarkar 2012).

Within this range, we display in Fig. 2 examples of eigenvalue spectra in the space (σ_i, σ_r) at f = 0.5, N = 1, and Pe = 0.1 for two sets of (k_x, k_z): (a) (0.445, 0) and (b) (0, 5). In Fig. 2a, the most unstable mode denoted by the triangle represents an inflectional instability mode. This mode has the growth rate σ_r = 0.1897 and zero frequency σ_i = 0 as equivalently reported in Deloncle et al. (2007) for the nondiffusive inflectional instability mode. At k_z = 0, we can derive a single 2nd-order ODE for $\widehat{\mathit{v}}$ from the continuity Eq. (19) and momentum Eqs. (20) and (21):

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\widehat{v}}{\mathrm{d}{y}^2}-(k_x^2+\frac{k_x{U}^{″}}{s})\widehat{v}=0,\end{array}$$(34)

Fig. 2. Eigenvalue spectra at (a) (kx, kz) = (0.445, 0) and (b) (kx, kz) = (0, 5) for f = 0.5, N = 1, and Pe = 0.1. The triangle in (a) and the square in (b) denote the most unstable growth rates of the inflectional and inertial instabilities, respectively.

(see also, Deloncle et al. 2007). It is important to note that the Eq. (34) does not depend on neither stratification, rotation, nor thermal diffusion, and the two-dimensional inflectional instability is consequently found to be independent of N (Deloncle et al. 2007) as well as f and Pe. We note that there also exist neutral waves (i.e., σ_r = 0) with nonzero frequency lying collectively in the range |σ_i|≲k_x = 0.445 as well as stable modes with σ_r <  0. In this study, however, we will consider only the unstable modes with zero frequency σ_i = 0. The eigenvalues in Fig. 2b at (k_x, k_z) = (0, 5) show different spectral behaviors from those in Fig. 2a: eigenvalues are distributed symmetrically with respect to σ_r = 0 for unstable and stable modes while neutral waves lie in the frequency range |σ_i|≲f = 0.5. The most unstable mode denoted by the square corresponds to the inertial instability mode with the growth rate σ_r = 0.4378 and zero frequency σ_i = 0.

In Fig. 3a, we display contours of the most unstable growth rate with thermal diffusion at Pe = 0.1 to locate the inflectional and inertial instabilities in the parameter space (k_x, k_z) for f = 0.5 and N = 1. The growth-rate contours are qualitatively similar to those in Arobone & Sarkar (2012), but we focus on the small-Pe regime. The frequency σ_i of the most unstable mode is zero thus not plotted as contours in the parameter space (k_x, k_z). Furthermore, we numerically verified that the inflectional instability is the most unstable in the two-dimensional case at (k_x, k_z) = (0.445, 0) for any values of N and Pe in the inertially-stable range (f ≥ 1 or f ≤ 0). On the other hand, in the inertially-unstable range 0 <  f <  1, the growth rate increases as k_z increases at k_x = 0.445 due to the inertial instability, and the most unstable growth rate of the inertial instability is found as k_z → ∞ at k_x = 0.

Fig. 3. Panel a: contours of the maximum growth rate σmax in the parameter space of (kx, kz) for f = 0.5, N = 1, Pe = 0.1. Panels b and c: corresponding growth rates σr versus vertical wavenumber kz at panel b: kx = 0 and panel c: kx = 0.445. For clarity, only the first four eigenvalue branches are plotted. The dashed line in panel c denotes the growth rate for rotating case f = 0. Symbols denote the parameters for eigenvalues in Fig. 2 and eigenfunctions displayed in Fig. 4.

Figure 3b shows an example of the growth rate σ_r as a function of k_z at k_x = 0. There is not only one unstable branch but a countless number of growth-rate branches (only the first four branches are shown in Fig. 3b for clarity). We see that the first branch is the most unstable and all the branches asymptote to certain values as k_z increases. Arobone & Sarkar (2012) argues that the inertial instability growth rate approaches ${\sigma }_{\max }=\sqrt{f(1-f)}$. In the next section, the reasons why there is an infinite number of branches and why they approach $\sqrt{f(1-f)}$ as k_z → ∞ will be explained by means of the WKBJ approximation. Figure 3c shows σ_r versus k_z at k_x = 0.445. We see that the first branch starts from σ_r = 0.1897 at k_z = 0 due to the inflectional instability, while other branches start at σ_r = 0 and increase with k_z. The increase of the growth rate σ_r with k_z occurs only in the inertially-unstable regime (0 <  f <  1) since the inflectional instability is stabilized as k_z increases in the inertially-stable regime (see the dashed line in Fig. 3c for f = 0).

Figure 4 shows examples of modes for the inflectional instability at (k_x, k_z) = (0.445, 0) in panels a and b, and the inertial instability at (k_x, k_z) = (0, 5) in panels c and d for the same parameters used in Fig. 3. For both instability modes, the mode shapes are normalized by the maximum of $|\widehat{\mathit{v}}|$. For Fig. 4, we show $\widehat{\mathit{v}}$ as a representative for the inflectional instability to be able to compare with the previous literature (Deloncle et al. 2007). For the inertial instability, we illustrate the behavior of $\widehat{T}$ since this is the quantity we will use for its asymptotic analysis in Sect. 4. The mode shape $\widehat{\mathit{v}}$ of the inflectional instability decreases exponentially as y → ±∞, and if plotted in physical space, we see that the perturbation $\stackrel{\sim }{\mathit{v}}(x,\mathit{y})$ is slightly inclined against the direction of the shear (see also, Arobone & Sarkar 2012). On the other hand, the mode shape $\widehat{T}$ of the inertial instability mode shows there exists a zero crossing for the absolute part of $\widehat{T}$ at y = 0 while decaying exponentially as y → ±∞. The temperature perturbation $\stackrel{\sim }{T}$ in the space (y, z) shows that the inertial instability mode has a wave pattern with a zero-crossing at y = 0. One zero-crossing corresponds to the first branch, which is the most unstable branch for given k_x and k_z, and higher branches have multiple zero-crossings accordingly.

Fig. 4. Examples of eigenmodes at (a, b) (kx, kz) = (0.445, 0) for the inflectional instability (triangle in Fig. 3) and (c, d) (kx, kz) = (0, 5) for the inertial instability (square in Fig. 3). Panels a and c: mode shapes of $\widehat{\mathit{v}}(\mathit{y})$ and $\widehat{T}(\mathit{y})$ (blue: real part, red: imaginary part, black: absolute value), respectively. Panels b and d: perturbations $\stackrel{\sim }{\mathit{v}}(x,\mathit{y})$ and $\stackrel{\sim }{T}(\mathit{y},z)$ in the physical space, respectively.

4. WKBJ analysis for the inertial instability

The maximum growth rate of the inertial instability is proposed to be ${\sigma }_{\max }=\sqrt{f(1-f)}$ (Arobone & Sarkar 2012). We also verified numerically that the most unstable growth rate is reached as k_z → ∞ at k_x = 0. But it is still not clear why this maximum growth rate σ_max is attained as k_z → ∞, and when the inertial instability occurs in parameter ranges of k_z, N, f, and Pe. In this section, we detail mathematical and physical interpretations of the inertial instability using the WKBJ approximation in the limit of large k_z at k_x = 0. Similar asymptotic analyses have been performed to understand the instability mechanisms in rotating shear flows (Park & Billant 2012; Park et al. 2017). To understand the effect of thermal diffusion, we perform an asymptotic analysis at two limits: Pe → ∞ (i.e., no thermal diffusion) and Pe → 0 (i.e., high thermal diffusivity). It is noticeable that the diffusion is often neglected by taking Pe → ∞ to understand geophysical flows in the atmosphere and oceans of the Earth (Yavneh et al. 2001; Park et al. 2018), while the high thermal diffusivity case with Pe → 0 is studied in astrophysical context for shear instability and mixing in stellar interiors (Lignières et al. 1999; Prat & Lignières 2014). In the following subsections, we provide explicit expressions of asymptotic dispersion relations for the inertial instability in stratified and rotating fluids in both limits of Pe.

4.1. The weak diffusion limit: Pe →∞

We first consider the 4th-order ODE (28) for k_x = 0:

$$\begin{array}{cc}& \frac{{\mathrm{d}}^4\widehat{T}}{\mathrm{d}{y}^4}-\frac{\mathrm{\Gamma }{}^{\prime }}{\mathrm{\Gamma }}\frac{{\mathrm{d}}^3\widehat{T}}{\mathrm{d}{y}^3}-[k_z^2(1-\frac{\mathrm{\Gamma }}{{\sigma }^2})+\mathrm{Pe}\left(\frac{N^2+{\sigma }^2}{\sigma }\right)]\frac{{\mathrm{d}}^2\widehat{T}}{\mathrm{d}{y}^2}\hfill \\ \hfill & +[k_z^2\frac{\mathrm{\Gamma }{}^{\prime }}{\mathrm{\Gamma }}+\mathrm{Pe}\left(\frac{N^2+{\sigma }^2}{\sigma }\right)\frac{\mathrm{\Gamma }{}^{\prime }}{\mathrm{\Gamma }}]\frac{\mathrm{d}\widehat{T}}{\mathrm{d}y}-[k_z^4\frac{\mathrm{\Gamma }}{{\sigma }^2}+k_z^2\mathrm{Pe}\left(\frac{\mathrm{\Gamma }}{\sigma }\right)]\widehat{T}=0.\hfill \end{array}$$(35)

In the limit Pe → ∞, we can rearrange Eq. (35) as

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\widehat{T}}{\mathrm{d}{y}^2}-\frac{\mathrm{\Gamma }{}^{\prime }}{\mathrm{\Gamma }}\frac{\mathrm{d}\widehat{T}}{\mathrm{d}y}+k_z^2\frac{\mathrm{\Gamma }}{N^2+{\sigma }^2}\widehat{T}=O\left(\frac{1}{\mathrm{Pe}}\right),\end{array}$$(36)

where higher-order derivatives in the 3rd and 4th orders are also of the order O(1/Pe). Provided that d/dy is of order unity, we can neglect the term on the right-hand side as Pe → ∞, and the Eq. (36) becomes the second-order ODE.

Applying the WKBJ approximation to (36) for large k_z:

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

we obtain

$$\begin{array}{c}\hfill \delta =\frac{1}{k_z},S_0^{{}^{\prime }2}=-\frac{\mathrm{\Gamma }}{N^2+{\sigma }^2},S_1^{\prime }=\frac{\mathrm{\Gamma }{}^{\prime }}{4\mathrm{\Gamma }}·\end{array}$$(38)

In this paper, we only consider the inertial instability mode which has the zero frequency (i.e., σ_i = 0 and σ² + N² >  0) thus the sign of Γ determines the behavior of the solutions. We get evanescent solutions if Γ <  0:

$$\begin{array}{cc}\hfill \widehat{T}(y)=& {(-\mathrm{\Gamma })}^{1/4}[A_{1}exp\left(k_z{\displaystyle {\int }_y\sqrt{\frac{-\mathrm{\Gamma }}{{\sigma }^2+N^2}}\mathrm{d}y}\right)\hfill \\ \hfill & +A_{2}exp(-k_z{\displaystyle {\int }_y\sqrt{\frac{-\mathrm{\Gamma }}{{\sigma }^2+N^2}}\mathrm{d}y})],\hfill \end{array}$$(39)

where A₁ and A₂ are constants, or wavelike solutions if Γ >  0:

$$\begin{array}{cc}\hfill \widehat{T}(y)=& {\mathrm{\Gamma }}^{1/4}[B_{1}exp\left(i{k}_z{\displaystyle {\int }_y\sqrt{\frac{\mathrm{\Gamma }}{{\sigma }^2+N^2}}\mathrm{d}y}\right)\hfill \\ \hfill & +B_{2}exp(-i{k}_z{\displaystyle {\int }_y\sqrt{\frac{\mathrm{\Gamma }}{{\sigma }^2+N^2}}\mathrm{d}y})],\hfill \end{array}$$(40)

where B₁ and B₂ are constants. These exponential behaviors depending on y change at turning points y_t where Γ(y_t) = 0. It is also convenient to introduce the turning growth rates σ_±:

$$\begin{array}{c}\hfill {\sigma }_{\pm }=\pm \mathfrak{R}\left[\sqrt{f(U^{\prime }-f)}\right],\end{array}$$(41)

which inform us where Γ = −σ² + f(U′−f) becomes zero. For instance, if the growth rate is σ = 0.348 that is a growth rate obtained numerically for f = 0.5, N = 1, Pe = ∞, k_x = 0, and k_z = 5, there exist two turning points: one at y = y_t− = −0.56 and the other at y = y_t+ = 0.56 (see Fig. 5). This implies that the solution in the regions y >  y_t+ and y <  y_t− is evanescent while the solution in the region y_t− <  y <  y_t+ is wavelike. We note that ${\sigma }_{\pm }=\pm \Re \left(\sqrt{-\mathrm{\Phi }}\right)$ so we can verify that the region where the Rayleigh’s discriminant Φ(y) = f(f − U′) is negative lies between the two turning points y_t±.

Fig. 5. Turning growth rates σ±(y) for f = 0.5 (black solid lines) and an example of the growth rate σ = 0.348 at (kx, kz) = (0, 5) for N = 1 and Pe = ∞ (dashed line). White and gray areas represent the regions where the solution is exponential and wavelike, respectively.

When the growth rate lies between the two turning growth rates σ₋ <  σ <  σ₊, this implies that the solutions are wavelike with Γ >  0 (gray area in Fig. 5), while the solutions are evanescent outside this range (white area in Fig. 5). For other values of σ, we can decide whether we can construct eigenfunctions. For example, if 0 <  σ <  max(σ₊) or min(σ₋) < σ <  0 like the case σ = 0.4232 <  max(σ₊) = 0.5 in Figs. 4c and d, the solution is exponential outside the turning points y_t± and wavelike in between. In this case, we can construct an eigenfunction which decays exponentially as y → ±∞ due to the presence of two turning points. On the other hand, if the growth rate is either σ >  max(σ₊) or σ <  min(σ₋), then there is no turning point and solutions are always evanescent for all y. Therefore, we must impose A₁ = 0 and A₂ = 0 to make the solution decaying with y → ∞ and y → −∞, respectively, and no eigenfunction can be constructed in this case. This also implies that the growth rate does not surpass $max({\sigma }_{+})=\sqrt{f(1-f)}$ to construct the inertial instability mode, which verifies the conjecture of Arobone & Sarkar (2012).

For an unstable mode that has the growth rate in the range 0 <  σ <  max(σ₊), we can further derive expressions for the asymptotic dispersion relation by performing a turning point analysis. We first consider the evanescent WKBJ solution that decays exponentially for y >  y_t+:

$$\begin{array}{c}\hfill \widehat{T}(y)=A_{\infty }{(-\mathrm{\Gamma })}^{1/4}exp(-k_z{\displaystyle {\int }_{y_{t+}}^y\sqrt{\frac{-\mathrm{\Gamma }}{{\sigma }^2+N^2}}\mathrm{d}y}),\end{array}$$(42)

where A_∞ is a constant. Around the turning point y_t+, the WKBJ solution (42) is no longer valid and we need to find a local solution that matches with the WKBJ solution in the range y_t− <  y <  y_t+. By considering a new scaling $\stackrel{\sim }{\mathit{y}}=(\mathit{y}-{\mathit{y}}_{t+})/ϵ$ and an approximation $\mathrm{\Gamma }(\mathit{y})\sim \mathrm{\Gamma }{\prime }_{t+}ϵ\stackrel{\sim }{\mathit{y}}$ where ${\mathrm{\Gamma }}_{t+}^{\prime }$ is the derivative of Γ in y at y = y_t+, we obtain the following local equation

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\widehat{T}}{\mathrm{d}{\stackrel{\sim }{y}}^2}-\frac{1}{\stackrel{\sim }{y}}\frac{\mathrm{d}\widehat{T}}{\mathrm{d}\stackrel{\sim }{y}}-\stackrel{\sim }{y}\widehat{T}=O(ϵ),\end{array}$$(43)

where $ϵ={[k_z^2(-\mathrm{\Gamma }{\prime }_{t+})/(N^2+{\sigma }^2)]}^{-1/3}$. The solution of this local Eq. (43) can be expressed in terms of derivatives of the Airy functions: $\widehat{T}(\stackrel{\sim }{\mathit{y}})=a_1\mathrm{Ai}\prime (\stackrel{\sim }{\mathit{y}})+b_1\mathrm{Bi}\prime (\stackrel{\sim }{\mathit{y}})$, where a₁ and b₁ are constants to be matched from the asymptotic behavior of the WKBJ solution (42) as y → y_t+. From the asymptotic behaviors of the Airy functions when $\stackrel{\sim }{\mathit{y}}\to +\infty$:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathrm{Ai}}^{\prime }(\stackrel{\sim }{y})\to -\frac{{\stackrel{\sim }{y}}^{1/4}}{2\sqrt{\pi }}exp(-\frac{2}{3}{\stackrel{\sim }{y}}^{3/2}),\hfill \\ \hfill & {\mathrm{Bi}}^{\prime }(\stackrel{\sim }{y})\to \frac{{\stackrel{\sim }{y}}^{1/4}}{\sqrt{\pi }}exp\left(\frac{2}{3}{\stackrel{\sim }{y}}^{3/2}\right),\hfill \end{array}\end{array}$$(44)

and when $\stackrel{\sim }{\mathit{y}}\to -\infty$:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathrm{Ai}}^{\prime }(\stackrel{\sim }{y})\to -\frac{{(-\stackrel{\sim }{y})}^{1/4}}{\sqrt{\pi }}cos(\frac{2}{3}{(-\stackrel{\sim }{y})}^{3/2}+\frac{\pi }{4}),\hfill \\ \hfill & {\mathrm{Bi}}^{\prime }(\stackrel{\sim }{y})\to \frac{{(-\stackrel{\sim }{y})}^{1/4}}{\sqrt{\pi }}sin(\frac{2}{3}{(-\stackrel{\sim }{y})}^{3/2}+\frac{\pi }{4}),\hfill \end{array}\end{array}$$(45)

(Abramowitz & Stegun 1972), we obtain the matched WKBJ solution in the region y_t− <  y <  y_t+:

$$\begin{array}{cc}\hfill \widehat{T}(y)=& {\mathrm{\Gamma }}^{1/4}[C_{+}exp\left(i{k}_z{\displaystyle {\int }_{y_{t+}}^y\sqrt{\frac{\mathrm{\Gamma }}{{\sigma }^2+N^2}}\mathrm{d}y}\right)\hfill \\ \hfill & +C_{-}exp(-i{k}_z{\displaystyle {\int }_{y_{t+}}^y\sqrt{\frac{\mathrm{\Gamma }}{{\sigma }^2+N^2}}\mathrm{d}y})],\hfill \end{array}$$(46)

where

$$\begin{array}{c}\hfill C_{+}=exp(-i\frac{\pi }{4})A_{\infty },C_{-}=exp\left(i\frac{\pi }{4}\right)A_{\infty }.\end{array}$$(47)

Similarly, the evanescent solution in y <  y_t− that decays exponentially as y → −∞:

$$\begin{array}{c}\hfill \widehat{T}(y)=A_{-\infty }{(-\mathrm{\Gamma })}^{1/4}exp\left(k_z{\displaystyle {\int }_{y_{t-}}^y\sqrt{\frac{-\mathrm{\Gamma }}{{\sigma }^2+N^2}}\mathrm{d}y}\right),\end{array}$$(48)

matches with the following solution in y_t− <  y <  y_t+ after the local solution around y_t− is considered:

$$\begin{array}{cc}\hfill \widehat{T}(y)=& {\mathrm{\Gamma }}^{1/4}[B_{+}exp\left(i{k}_z{\displaystyle {\int }_{y_{t-}}^y\sqrt{\frac{\mathrm{\Gamma }}{{\sigma }^2+N^2}}\mathrm{d}y}\right)\hfill \\ \hfill & +B_{-}exp(-i{k}_z{\displaystyle {\int }_{y_{t-}}^y\sqrt{\frac{\mathrm{\Gamma }}{{\sigma }^2+N^2}}\mathrm{d}y})],\hfill \end{array}$$(49)

where

$$\begin{array}{c}\hfill B_{+}=exp\left(i\frac{\pi }{4}\right)A_{-\infty },B_{-}=exp(-i\frac{\pi }{4})A_{-\infty }.\end{array}$$(50)

Matching the wavelike solutions (46) and (49):

$$\begin{array}{c}\hfill \frac{B_{-}}{B_{+}}exp(-2i{k}_z{\displaystyle {\int }_{y_{t-}}^{y_{t+}}\sqrt{\frac{\mathrm{\Gamma }}{{\sigma }^2+N^2}}\mathrm{d}y})=\frac{C_{-}}{C_{+}},\end{array}$$(51)

that is

$$\begin{array}{c}\hfill exp\left(2i{k}_z{\displaystyle {\int }_{y_{t-}}^{y_{t+}}\sqrt{\frac{\mathrm{\Gamma }}{{\sigma }^2+N^2}}\mathrm{d}y}\right)=exp(-i\pi ),\end{array}$$(52)

we obtain the dispersion relation in the form of a quantization formula:

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

where m is the branch number with a positive integer. This quantized dispersion relation implies that there exist an infinite number of discrete growth-rate branches for the inertial instability, as verified in our numerical results (Figs. 3b and c).

The integral on the left-hand side of (53) should become zero as k_z → ∞ since the right-hand side term is fixed with a finite m. This implies that the two turning points should approach each other as k_z goes to infinity while Γ → 0 (i.e., y_t− → y_t+ → 0). In addition to this property and the Taylor expansion of Γ around y = 0:

$$\begin{array}{c}\hfill \mathrm{\Gamma }\simeq -{\sigma }^2-f^2+f{U}^{\prime }{|}_{y=0}+\frac{f{U}^{{}^{″\prime }}{|}_{y=0}}{2}{y}^2,\end{array}$$(54)

we also assume that the growth rate can be expanded as $\sigma ={\sigma }_0-{\sigma }_1{k}_z^{-p}+O(k_z^{-2p})$ where p is a positive integer. Applying this growth-rate expansion to the dispersion relation (53), we find that

$$\begin{array}{c}\hfill p=1,y_{t\pm }^2=\frac{2{\sigma }_0{\sigma }_1}{f{k}_z},{\sigma }_0=\sqrt{f(1-f)},\end{array}$$(55)

and we obtain the first-order term σ₁, which leads to a more explicit dispersion relation for very large k_z:

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

where

$$\begin{array}{c}\hfill {\sigma }_1=\frac{\sqrt{f}}{{\sigma }_0}(m-\frac{1}{2})\sqrt{f(1-f)+N^2}.\end{array}$$(57)

Since σ₁ is a positive value, we see that the maximum growth rate ${\sigma }_{max}={\sigma }_0=\sqrt{f(1-f)}$ is achieved at k_z → ∞ for any values of m. The dispersion relation (56) also implies that the inertial instability only exists in the range 0 <  f <  1 to have a positive real σ₀. While the maximum growth rate σ_max is independent of N, the term σ₁ depends on the stratification and increases as N increases; the growth rate thus decreases with N. We see in Fig. 6 that the asymptotic dispersion relation (56) matches well with the numerical results as k_z increases.

Fig. 6. Growth rates of the first branch for Pe = ∞ (solid line) and Pe = 0.01 (dashed line) at f = 0.5 and N = 1 with predictions from the WKBJ dispersion relations (56) for Pe = ∞ (solid with filled circles) and (72) for Pe → 0 (dashed with empty circles).

The asymptotic dispersion relation (56) is a rich source of information regarding the inertial instability that covers a wide range of parameters f and N. In Fig. 7, we see how the growth rate contours change in the parameter space (f, N) for Pe = ∞, k_x = 0, and k_z = 20. At a fixed f, the growth rate decreases with N. The contours of the asymptotic growth rate match very well with our numerical results especially for higher values of the growth rate. The better agreement at higher growth rates follows from our expansion of (53). Specifically, the turning points are assumed to be close to each other, which is equivalent to assuming that the growth rate is close to its maximum value.

Fig. 7. Growth rate contours in the parameter space (f, N) for Pe = ∞, kx = 0, kz = 20 (solid lines), and predictions from the WKBJ dispersion relation (56) (dashed lines).

4.2. The strong diffusion limit: Pe → 0

Now we investigate the limit Pe → 0 for high thermal diffusivity. The 4th-order ODE (35) can be expressed as

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^4\widehat{T}}{\mathrm{d}{y}^4}-\frac{\mathrm{\Gamma }{}^{\prime }}{\mathrm{\Gamma }}\frac{{\mathrm{d}}^3\widehat{T}}{\mathrm{d}{y}^3}-k_z^2(1-\frac{\mathrm{\Gamma }}{{\sigma }^2})\frac{{\mathrm{d}}^2\widehat{T}}{\mathrm{d}{y}^2}+k_z^2\frac{\mathrm{\Gamma }{}^{\prime }}{\mathrm{\Gamma }}\frac{\mathrm{d}\widehat{T}}{\mathrm{d}y}-k_z^4\frac{\mathrm{\Gamma }}{{\sigma }^2}\widehat{T}=O(\mathrm{Pe}).\end{array}$$(58)

Applying the WKBJ approximation (37), we obtain the following relations at leading order:

$$\begin{array}{c}\hfill \delta =\frac{1}{k_z},S_0^{{}^{\prime }4}-(1-\frac{\mathrm{\Gamma }}{{\sigma }^2})S_0^{{}^{\prime }2}-\frac{\mathrm{\Gamma }}{{\sigma }^2}=0,\end{array}$$(59)

and at first order:

$$\begin{array}{c}\hfill S_1^{\prime }[4S_0^{{}^{\prime }2}-2(1-\frac{\mathrm{\Gamma }}{{\sigma }^2})]+\frac{S_0^{{}^{″}}}{S_0^{\prime }}(6S_0^{{}^{\prime }2}-1+\frac{\mathrm{\Gamma }}{{\sigma }^2})+\frac{\mathrm{\Gamma }{}^{\prime }}{\mathrm{\Gamma }}(1-S_0^{{}^{\prime }2})=0.\end{array}$$(60)

The terms S₀ and S₁ satisfy

$$\begin{array}{c}\hfill S_0^{{}^{\prime }2}=-\frac{\mathrm{\Gamma }}{{\sigma }^2},S_1^{{}^{\prime }}=\frac{\mathrm{\Gamma }{}^{\prime }}{4\mathrm{\Gamma }}-\frac{\mathrm{\Gamma }{}^{\prime }}{\mathrm{\Gamma }+{\sigma }^2},\end{array}$$(61)

or

$$\begin{array}{c}\hfill S_0^{{}^{\prime }2}=1,S_1^{{}^{\prime }}=0.\end{array}$$(62)

The general WKBJ solution with $S_0^{\prime 2}(\mathit{y})=1$ is

$$\begin{array}{c}\hfill \widehat{T}(y)=D_{3}exp(k_{z}y)+D_{4}exp(-k_{z}y),\end{array}$$(63)

but it has no turning point and the two solutions are always exponentially decaying or increasing as y → ±∞. Therefore, we must impose D₃ = D₄ = 0 to construct an eigenfunction. With $S_0^{\prime 2}=-\mathrm{\Gamma }/{\sigma }^2$, the WKBJ solution is wavelike if Γ >  0:

$$\begin{array}{cc}\hfill \widehat{T}(y)=& \frac{{\mathrm{\Gamma }}^{1/4}}{|{\sigma }^2+\mathrm{\Gamma }|}[D_{1}exp\left(i{k}_z{\displaystyle {\int }_y\frac{\sqrt{\mathrm{\Gamma }}}{\sigma }\mathrm{d}y}\right)\hfill \\ \hfill & +D_{2}exp(-i{k}_z{\displaystyle {\int }_y\frac{\sqrt{\mathrm{\Gamma }}}{\sigma }\mathrm{d}y})],\hfill \end{array}$$(64)

or evanescent if Γ <  0:

$$\begin{array}{cc}\hfill \widehat{T}(y)=& \frac{{(-\mathrm{\Gamma })}^{1/4}}{|{\sigma }^2+\mathrm{\Gamma }|}[E_{1}exp\left(k_z{\displaystyle {\int }_y\frac{\sqrt{-\mathrm{\Gamma }}}{\sigma }\mathrm{d}y}\right)\hfill \\ \hfill & +E_{2}exp(-k_z{\displaystyle {\int }_y\frac{\sqrt{-\mathrm{\Gamma }}}{\sigma }\mathrm{d}y})].\hfill \end{array}$$(65)

These solutions in (64) and (65) now have turning points where Γ = 0, which is the same as in the case without thermal diffusion at Pe = ∞. Therefore, the turning growth rates σ_±(y) of the solutions (64)–(65) are essentially the same as the expression (41), and we can construct an eigenfunction if the growth rate σ lies in the ranges min(σ₋) < σ <  0 or 0 <  σ <  max(σ₊).

For Pe → 0, we can also perform a turning point analysis in order to obtain the asymptotic dispersion relation. We consider the WKBJ solution that decays exponentially for y >  y_t+:

$$\begin{array}{c}\hfill \widehat{T}(y)=E_{\infty }\frac{{(-\mathrm{\Gamma })}^{1/4}}{|{\sigma }^2+\mathrm{\Gamma }|}exp(-k_z{\displaystyle {\int }_{y_{t+}}^y\frac{\sqrt{-\mathrm{\Gamma }}}{\sigma }\mathrm{d}y}).\end{array}$$(66)

Around the turning point y_t+, we apply a new scaled coordinate $\stackrel{\sim }{\mathit{y}}=(\mathit{y}-{\mathit{y}}_{t+})/ϵ$ and obtain the following local equation

$$\begin{array}{c}\hfill ϵ\frac{{\mathrm{d}}^4\widehat{T}}{\mathrm{d}{\stackrel{\sim }{y}}^4}-\frac{ϵ}{\stackrel{\sim }{y}}\frac{{\mathrm{d}}^3\widehat{T}}{\mathrm{d}{\stackrel{\sim }{y}}^3}-k_z^2{ϵ}^3\frac{{\mathrm{d}}^2\widehat{T}}{\mathrm{d}{\stackrel{\sim }{y}}^2}+\frac{k_z^2{ϵ}^3}{\stackrel{\sim }{y}}\frac{\mathrm{d}\widehat{T}}{\mathrm{d}\stackrel{\sim }{y}}+\frac{k_z^4{ϵ}^6(-\mathrm{\Gamma }{{}^{\prime }}_{t+})}{{\sigma }^2}\widehat{T}=0.\end{array}$$(67)

Here, we take $ϵ={[k_z^2(-\mathrm{\Gamma }{\prime }_{t+})/{\sigma }^2]}^{-1/3}$ to balance the equation, and the 3rd- and 4th-order derivatives become negligible since they are of order O(ϵ). The final local equation becomes the same as (43) thus the local solution is the sum of derivatives of the Airy functions. By matching the solution around y_t+ using the asymptotic behaviors of the Airy functions as $\stackrel{\sim }{\mathit{y}}\to \pm \infty$ (Abramowitz & Stegun 1972) similar to in the previous section, we obtain the WKBJ solution for y_t− <  y <  y_t+

$$\begin{array}{cc}\hfill \widehat{T}(y)=& \frac{{\mathrm{\Gamma }}^{1/4}}{|{\sigma }^2+\mathrm{\Gamma }|}[F_{+}exp\left(i{k}_z{\displaystyle {\int }_{y_{t+}}^y\frac{\sqrt{\mathrm{\Gamma }}}{\sigma }\mathrm{d}y}\right)\hfill \\ \hfill & +F_{-}exp(-i{k}_z{\displaystyle {\int }_{y_{t+}}^y\frac{\sqrt{\mathrm{\Gamma }}}{\sigma }\mathrm{d}y})],\hfill \end{array}$$(68)

where

$$\begin{array}{c}\hfill F_{+}=exp(-i\frac{\pi }{4})E_{\infty },F_{-}=exp\left(i\frac{\pi }{4}\right)E_{\infty }.\end{array}$$(69)

By matching the solution (68) with the solution for y <  y_t−:

$$\begin{array}{c}\hfill \widehat{T}(y)=E_{-\infty }\frac{{(-\mathrm{\Gamma })}^{1/4}}{|{\sigma }^2+\mathrm{\Gamma }|}exp\left(k_z{\displaystyle {\int }_{y_{t-}}^y\frac{\sqrt{-\mathrm{\Gamma }}}{\sigma }\mathrm{d}y}\right),\end{array}$$(70)

we obtain the following dispersion relation in a quantization form

$$\begin{array}{c}\hfill k_z{\displaystyle {\int }_{y_{t-}}^{y_{t+}}\frac{\sqrt{\mathrm{\Gamma }}}{\sigma }\mathrm{d}y=(m_0-\frac{1}{2})\pi ,}\end{array}$$(71)

where m₀ is the positive integer indicating the branch number. Similarly to the previous section, we Taylor-expand the growth rate σ and we obtain an explicit asymptotic dispersion relation for σ as a function of k_z:

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

where

$$\begin{array}{c}\hfill {\sigma }_{0,0}=\sqrt{f(1-f)},{\sigma }_{1,0}=(m_0-\frac{1}{2})\sqrt{f}.\end{array}$$(73)

We see in Fig. 6 that a numerical result at small Pe = 0.01 matches very well with the asymptotic growth rate (72) in the limit Pe → 0 as k_z increases. It is important to note that the expression of the maximum growth rate is the same as (56) in nondiffusive fluids, and the inertial instability in stratified-rotating fluids with high thermal diffusivity also occurs in the regime 0 <  f <  1. In contrast to the nondiffusive case, the term σ_1, 0 is independent of the stratification N.

4.3. Comparison of the growth rates

It is noteworthy that the asymptotic dispersion relation (72) in fluids with high thermal diffusivity is independent of the stratification, and the expression (56) without thermal diffusion becomes identical to (72) as N → 0. This implies that high thermal diffusivity suppresses the effect of the stable stratification. Thus, stratified fluids with high thermal diffusivity behave like unstratified fluids. The term σ_1, 0 at the first order is always smaller than σ₁ in (57); therefore, the ratio γ₁ between the two terms at first order is always larger than unity for positive N if we consider the same branch m = m₀:

$$\begin{array}{c}\hfill {\gamma }_1=\frac{{\sigma }_1}{{\sigma }_{1,0}}=\sqrt{\frac{f(1-f)+N^2}{f(1-f)}}>1.\end{array}$$(74)

The growth rate of inertial instability in thermally-diffusive fluids is therefore always larger than that in nondiffusive fluids. The effect of high thermal diffusivity has already been reported in previous literature for the vertical shear instability (Lignières 1999; Prat & Lignières 2013) but the above expression shows explicitly and quantitatively how much the inertial instability is destabilized by a high thermal diffusivity as shown in Fig. 6.

5. Parametric study

5.1. Effects of Pe on the inflectional instability

While the inertial instability is accessible by the WKBJ approximation, the inflectional instability should be investigated numerically since it occurs at low k_x and k_z. In Fig. 8a, we see the growth rate σ_r of the inflectional instability as a function of k_z for different values of the Péclet number Pe at k_x = 0.445, N = 2, and f = 0, the regime where only the inflectional instability exists. While the two-dimensional growth rate at k_z = 0 remains the same regardless of Pe, the growth rate σ_r decreases with k_z for any values of Pe, and the stabilization occurs faster as Pe decreases to zero (i.e., as the thermal diffusivity increases). We see that growth-rate curves (solid lines) approach the growth-rate curve for the unstratified case at N = 0 (dashed line) as Pe goes to zero. This implies that high thermal diffusivity suppresses the effect of stratification on the inflectional instability. The suppression of the three-dimensional inflectional instability by the thermal diffusion is also found in fast rotating regime at N/f = 0.1 (Fig. 8b). We see that the instability is sustained at smaller k_z for the fast rotation with small N/f, while the stabilization with k_z is faster as the thermal diffusivity increases.

Fig. 8. Growth rate σr versus vertical wavenumber kz at kx = 0.445 for different Péclet numbers Pe: Pe = ∞ (black solid line), Pe = {100,10,5,1,0.5,0.1,0.01} (gray solid lines) and other parameters: panel a: N = 2, N/f = ∞, panel b: N = 2, N/f = 0.1, panel c: N = 0.5, N/f = 1. Dashed lines denote the growth rate for unstratified case N = 0 at panel a: f = 0, panel b: f = 20, panel c: f = 0.5. Arrows indicate the direction of decreasing Pe.

In Fig. 9, we investigate the effects of both Pe and N on the three-dimensional inflectional instability in the inertially-stable regime at f = 0, k_x = 0.445, and k_z = 1. We see in Fig. 9a that the growth rate increases as the stratification becomes strong while the strong thermal diffusion with small Pe suppresses the instability. For strong stratification N ≥ 10, it is also observed that the shape of growth-rate curves is similar and they just have different onsets of stabilization at lower Péclet numbers (i.e., higher diffusivity) as N increases. Due to this resemblance of growth-rate curves for strong stratification, we plot again the growth rate σ_r versus the rescaled parameter PeN² in Fig. 9b. We see that the rescaled growth rate curves collapse well for high N. This selfsimilarity represented by the rescaled parameter PeN² is reminiscent of the Richardson–Péclet number RiPe in the small-Péclet-number approximation used for vertically sheared flows in stratified and thermally-diffusive fluids (Lignières et al. 1999; Prat & Lignières 2013).

Fig. 9. Growth rate of the inflectional instability versus panel a: Pe and panel b: PeN2 for different values of N at kx = 0.445, kz = 1, and f = 0. Blue line represents the growth rate computed from the eigenvalue problem (81) under the small-Pe approximation.

Following the small-Pe approximation (Lignières 1999), we introduce the rescaled temperature deviation $\widehat{\theta }$:

$$\begin{array}{c}\hfill \widehat{\theta }=\frac{\widehat{T}}{\mathrm{Pe}}·\end{array}$$(75)

In the limit of small Pe, we obtain the set of linear stability equations at leading order:

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

$$\begin{array}{cc}& (\sigma +i{k}_{x}U)\widehat{u}+(U^{\prime }-f)\widehat{v}=-i{k}_x\widehat{p},\hfill \end{array}$$(77)

$$\begin{array}{cc}& (\sigma +i{k}_{x}U)\widehat{v}+f\widehat{u}=-\frac{\mathrm{d}\widehat{p}}{\mathrm{d}y},\hfill \end{array}$$(78)

$$\begin{array}{cc}& (\sigma +i{k}_{x}U)\widehat{w}=-i{k}_z\widehat{p}+\mathrm{Pe}{N}^2\widehat{\theta },\hfill \end{array}$$(79)

$$\begin{array}{cc}& \widehat{w}={\widehat{\mathrm{\nabla }}}^2\widehat{\theta }.\hfill \end{array}$$(80)

Equations (76)–(80) can be simplified into an eigenvalue problem for $\widehat{\mathit{v}}$ and $\widehat{\theta }$: