Issue 
A&A
Volume 646, February 2021



Article Number  A64  
Number of page(s)  19  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202038654  
Published online  08 February 2021 
Horizontal shear instabilities in rotating stellar radiation zones
II. Effects of the full Coriolis acceleration
^{1}
AIM, CEA, CNRS, Université ParisSaclay, Université Paris Diderot, Sorbonne Paris Cité, 91191 GifsurYvette, France
email: junho.park@cea.fr
^{2}
Fluid and Complex Systems Research Centre, Coventry University, Coventry CV1 5FB, UK
Received:
14
June
2020
Accepted:
27
November
2020
Context. Stellar interiors are the seat of efficient transport of angular momentum all along their evolution. In this context, understanding the dependence of the turbulent transport triggered by the instabilities of the vertical and horizontal shears of the differential rotation in stellar radiation zones as a function of their rotation, stratification, and thermal diffusivity is mandatory. Indeed, it constitutes one of the cornerstones of the rotational transport and mixing theory, which is implemented in stellar evolution codes to predict the rotational and chemical evolutions of stars.
Aims. We investigate horizontal shear instabilities in rotating stellar radiation zones by considering the full Coriolis acceleration with both the dimensionless horizontal Coriolis component f̃ and the vertical component f.
Methods. We performed a linear stability analysis using linearized equations derived from the NavierStokes and heat transport equations in the rotating nontraditional fplane. We considered a horizontal shear flow with a hyperbolic tangent profile as the base flow. The linear stability was analyzed numerically in wide ranges of parameters, and we performed an asymptotic analysis for large vertical wavenumbers using the WentzelKramersBrillouinJeffreys (WKBJ) approximation for nondiffusive and highlydiffusive fluids.
Results. As in the traditional fplane approximation, we identify two types of instabilities: the inflectional and inertial instabilities. The inflectional instability is destabilized as f̃ increases and its maximum growth rate increases significantly, while the thermal diffusivity stabilizes the inflectional instability similarly to the traditional case. The inertial instability is also strongly affected; for instance, the inertially unstable regime is also extended in the nondiffusive limit as 0 < f < 1 + f̃ ^{2}/N^{2}, where N is the dimensionless BruntVäisälä frequency. More strikingly, in the high thermal diffusivity limit, it is always inertially unstable at any colatitude θ except at the poles (i.e., 0° < θ < 180°). We also derived the critical Reynolds numbers for the inertial instability using the asymptotic dispersion relations obtained from the WKBJ analysis. Using the asymptotic and numerical results, we propose a prescription for the effective turbulent viscosities induced by the inertial and inflectional instabilities that can be possibly used in stellar evolution models. The characteristic time of this turbulence is short enough so that it is efficient to redistribute angular momentum and to mix chemicals in stellar radiation zones.
Key words: hydrodynamics / turbulence / stars: rotation / stars: evolution
© J. Park et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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., UdDoula 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 abinitio 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 spacebased 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 HertzsprungRussell 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), stochasticallyexcited internal gravity waves (e.g., Talon & Charbonnel 2005; Rogers 2015; Pinçon et al. 2017), and mixed gravitoacoustic 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 stateoftheart 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 nearinertial 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 gravitoinertial waves, inertial waves, and waveinduced 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 premainsequence, lowmass stars or in the radiative envelope of rapidly rotating upper mainsequence stars. These regimes should be treated properly to build robust one or twodimensional (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 NavierStokes equations when the fluid is stably stratified and thermally diffusive in the rotating nontraditional fplane 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 WentzelKramersBrillouingJeffreys (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 highdiffusivity 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 2Dsecular 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. NavierStokes equations and base steady state
We consider the NavierStokes 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:
where u = (u,v,w) is the velocity, p is the pressure, Θ is the temperature deviation from the reference temperature T_{0}, f = (0,f_{h, 0},f_{v, 0}) is the Coriolis vector with the horizontal (latitudinal) Coriolis component f_{h, 0} = 2Ω_{0} sin θ and the vertical Coriolis component f_{v, 0} = 2Ω_{0} cos θ where is the stellar rotation rate, θ is the colatitude, ρ_{0} is the reference density, g = (0, 0, −g) is the gravity, ν_{0} is the reference viscosity, κ_{0} 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 ∇^{2} = ∂^{2}/∂x^{2} + ∂^{2}/∂y^{2} + ∂^{2}/∂z^{2} denotes the Laplacian operator. Figure 1 illustrates the local coordinate system of the horizontal shear flow in the radiative zone when the nontraditional fplane 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:
Fig. 1. (a) Schematic of the radiative and convective zones, colored as yellow and orange, respectively, for the configuration of lowmass stars rotating with angular speed Ω_{0}. (b) Horizontal shear flow U(y) in a local nontraditional fplane at a colatitude θ in the radiative zone; z_{c} denotes the transition altitude between the radiative and convective zones. 
where U_{0} and L_{0} 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 thermalwind balance between U(y) and the base temperature is satisfied as follows:
Such a thermalwind 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 fplane (Wang et al. 2014). In addition to the thermalwind balance, the base temperature is considered to be stably stratified with a linearly increasing profile in z; therefore, it has the form
where is the difference in base temperature along the vertical distance Δz at a given y.
2.2. Linearized equations
We consider the velocity perturbation , the pressure perturbation where P is the base pressure, and the temperature perturbation . We nondimensionalize variables in Eqs. (1)–(3) by considering the length scale L_{0}, the velocity scale U_{0}, the time scale as t_{0} = L_{0}/U_{0}, the pressure scale as , and the temperature scale as . For infinitesimally small perturbations, we have the following nondimensional linearized equations:
where the prime symbol (′) denotes the derivative with respect to y, f = 2Ω cos θ and are the dimensionless vertical and horizontal Coriolis parameters, respectively, Ω is the dimensionless stellar rotation rate, Re is the Reynolds number
N is the dimensionless BruntVäisälä frequency^{1}
and Pe is the Péclet number
In Eq. (11), the term on the lefthand side appears due to the thermalwind balance (5).
To derive linear stability equations, we apply a normal mode expansion to perturbations as
where , and are the mode shapes of velocity, pressure, and temperature perturbations, respectively, i^{2} = −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
where with . By eliminating using the continuity Eq. (16), Eqs. (16)–(20) can be gathered into a matrix form as an eigenvalue problem:
where 𝒜 and ℬ are the operator matrices expressed as
where
We discretize numerically the operators 𝒜 and ℬ in the ydirection using the rational Chebyshev function that maps the Chebyshev domain y_{cheb} ∈ ( − 1, 1) onto the physical space y ∈ ( − ∞, ∞) via the mapping where L_{map} is the mapping factor (Deloncle et al. 2007; Park et al. 2020). To distinguish physical and numericallyconverged 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 highlyoscillatory spurious modes from physical modes that decays smoothly at boundaries as y → ±∞. Furthermore, we chose the number of collocation points in the ydirection 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 .
While the stability of horizontal shear flows in stratifiedrotating fluids in the traditional fplane 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 . Instead, there exist two separate symmetry conditions on σ in terms of k_{x} and k_{z} as follows:
and
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
Equations (16)–(20) can be simplified into ordinary differential equations (ODEs) for in particular cases. For instance, in the inviscid and nondiffusive limits (Re → ∞ and Pe → ∞), we can derive the single secondorder ODE for as follows:
where
s = −iσ + k_{x}U is the Dopplershifted frequency, ,
and
For finite Pe in the inviscid limit, we can find for k_{x} = 0 the single fourthorder ODE as follows:
where
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 (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 gravitoinertial 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 (i.e., Ω = 1.031 and θ = 76°) for (a) (k_{x}, k_{z}) = (0.25, 0), (b) (k_{x}, k_{z}) = (0, 6), and (c) (k_{x}, k_{z}) = (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 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 fplane (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 : (a) the mode shape and (b) velocity v in the physical space (x, y) at z = 0 for (k_{x}, k_{z}) = (0.25, 0), and σ = 0.0712, (c) the mode shape and (d) velocity v in the physical space (y, z) at x = 0 for (k_{x}, k_{z}) = (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 , respectively. 
Figure 2b shows an eigenvalue spectrum that has nearlyneutral 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 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 Dopplershifted 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 followup 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 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 inertiallystable regime at f = 0 according to the criterion in the traditional approximation (Arobone & Sarkar 2012). If we consider (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 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 (k_{x}, k_{z}) (solid lines) for Re = ∞ and N = 1 for (a) , (b) , and (c) . 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 fplane approximation where . 
In Fig. 4b, we now consider a finite Péclet number Pe = 1. At (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 still remains as σ_{max} = 0.1897 at k_{x} = 0.445. As increases, we see that the inflectionallyunstable regime at is even more shrunk than the unstable regime at . 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 threedimensional inflectional instability is now more unstable than the twodimensional inflectional instability at finite Pe and . 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 inertiallystable regime in the traditional approximation) when the horizontal Coriolis component becomes positive and the fluid is thermally diffusive.
We now consider in Fig. 4c the case f = 0.5, which belongs to the inertiallyunstable regime in the traditional fplane. We see a clear difference between growthrate contours at (i.e., on the pole with 2Ω = 0.5) and those at (i.e., 2Ω ≃ 2.062 and θ ≃ 76°). While the growth rate smoothly increases as k_{z} increases in the traditional case at , the regime of the inflectional instability is well separated from the regime of the inertial instability for . The inflectional instability at 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 fplane 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 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 secondorder ODE for
Fig. 5. Inviscid growth rate of the inflectional instability for various parameter sets of and N at f = 0 (i.e., θ = 90°), Pe = ∞, and k_{z} = 0. 
We clearly see that Eq. (39) is still independent from f but it depends on N if . For a fixed N, we see in Fig. 5 that the maximum growth rate increases with , while it decreases with N at a fixed . 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 , it is noticeable that the instability can sustain for k_{x} > 1 as increases.
Figure 6 shows growth rate curves for various values of Pe at N = 1, f = 0, k_{x} = 0, and (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 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 k_{x} = 0 for f = 0 and , 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 increases. The maximum growth rate of the inflectional instability increases with and depends on N. Similarly to the traditional case at , thermal diffusion destabilizes the inertial instability when . In the following section, we discuss more thoroughly how the inertial instability depends on physical parameters such as 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 highlydiffusive (Pe → 0) cases to understand how the inertial instability is modified as the horizontal Coriolis parameter increases.
4. Asymptotic description of the inertial instability
Our previous work in Park et al. (2020) with the traditional approximation () 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 is reached as k_{z} → ∞. In the traditional approximation, σ_{max} is related to the epicyclic frequency ω_{ep}, which satisfies and leads to the SolbergHøiland criterion for stability when (Solberg 1936; Høiland 1941). Besides, σ_{max} is analogous to the growth rate suggested in the limit Pr → 0 for the GoldreichSchubertFricke (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 firstorder 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 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 higlydiffusive case as Pe → 0. This analysis is an extension of the work by Park et al. (2020) in the traditional fplane 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 secondorder ODE for in Eq. (30) at k_{x} = 0 becomes
We note that the term multiplied by the first derivative 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 prior to the WKBJ analysis
Then Eq. (40) becomes
which is now in the form of the Poincaré equation. Applying to Eq. (42) the WKBJ approximation of for large k_{z}:
we get
where
Since , the exponential behavior of determined by the sign of S_{0} depends on the sign of . On the one hand, the solution is evanescent if :
where A_{1} and A_{2} are constants, and . On the other hand, the solution is wavelike if :
where B_{1} and B_{2} are constants.
The wavelike or evanescent solution behavior changes at turning points y_{t} where . To find y_{t}, it is useful to define the function σ_{t}(y) such that
(see also, Park & Billant 2012, 2013; Park et al. 2017). By finding σ_{t}(y) = σ at a given σ, we can find turning points where . An example of σ_{t}(y) is shown in Fig. 7 for f = 0.5, , and N = 1. The gray area denotes where the WKBJ solution is wavelike (i.e., ) while the white area denotes where the solution is evanescent (i.e., ). 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 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 is evanescent everywhere with no turning point and we cannot construct an eigenfunction 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
Fig. 7. Function σ_{t}(y) (solid lines) at N = 1 for (i.e., Ω = 1.031 and θ = 76°). Gray and white areas denote the regimes where 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+}:
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 where ϵ is a small parameter defined as , is the derivative of at y_{t+}, and we assume at leading order that around the turning point y_{t+}. The following local equation is obtained
Solutions of the local equation (51) are the Airy functions: , where a_{1} and b_{1} are constants (Abramowitz & Stegun 1972). From the asymptotic behavior of the Airy functions as and of the WKBJ solution as y → y_{t+}, we obtain the matched WKBJ solution in the region y_{t−} < y < y_{t+}:
where constants C_{±} satisfy
Similarly, we consider the WKBJ solution below the lower turning point y < y_{t−} decaying exponentially as y → −∞:
After the local analysis around y_{t−}, we find the matched wavelike solution in the range y_{t−} < y < y_{t+}:
where constants B_{±} satisfy
Finally, matching the wavelike solutions (52) and (55) between two turning points leads to the following dispersion relation in the quantized form
where m is the positive integer denoting the mode number. This quantized dispersion relation becomes equivalent to that in the traditional approximation as becomes zero (Park et al. 2020).
The righthand side of Eq. (57) is always fixed at a finite m so that the integral on the lefthand side of Eq. (57) should go to zero as . This implies that the two turning points converge to zero as . Using this property, we can find a more explicit dispersion relation in the Taylor expansion form for the growth rate σ as
where
and
Since σ_{1} is positive, the first mode with m = 1 is the most unstable mode. There are also higherorder 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 σ = σ_{0} is the maximum growth rate σ_{max} as k_{z} → ∞, and σ_{0} equals to σ_{t, max} of Eq. (49). While the maximum growth rate in the traditional fplane approximation is , which depends only on f, we see that the maximum growth rate σ_{0} now depends on f, , and N altogether. For instance, at a fixed N, the maximum growth rate σ_{0} increases with as shown in Fig. 8a. Compared to the inertiallyunstable 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 :
Fig. 8. (a) Contours of the maximum growth rate σ_{0} from Eq. (59) in the parameter space for N = 2. Black dashed lines represent the upper limit of the instability from Eq. (61) for different N, and the white dashed line represents f_{max} from (64) where the maximum growth rate is attained. (b) Contours of σ_{0} from Eq. (59) in the parameter space at f = 1. (c) Contours of the firstorder term σ_{1} in the parameter space for N = 2 and the first branch at m = 1. 
If we use the relations f = 2Ω cos θ and , we have the corresponding range
which is equivalent to
if 2Ω > 0.
Moreover, at the vertical Coriolis parameter f = f_{max}
we can find the maximum growth rate at fixed and N as
In Fig. 8b, we see that the maximum growth rate σ_{0} at fixed f and decreases as N increases, which highlights the stabilizing role of the stratification. We also display in Fig. 8c the dependence of the firstorder term σ_{1} on f and at N = 2. In this case, σ_{1} increases with f at a fixed and goes to infinity as f reaches its upper limit since σ_{0} is in the denominator of σ_{1} and σ_{0} → 0 as .
Figure 9 shows the growth rate σ versus the vertical wavenumber k_{z} for different values of 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 , and the growth rate increases with , 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 and N at f = 1, k_{x} = 0, Re = ∞, and Pe = ∞. 
4.2. The WKBJ formulation in the limit Pe→0
In this subsection, we now turn our attention to the highlydiffusive 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 fourthorder ODE (34) for the WKBJ analysis. As performed in the previous subsection, we introduce for convenience the variable
to understand more clearly the exponential or wavelike behavior around turning points. The fourthorder ODE (34) can be rewritten in terms of as follows:
where the righthand side of order will be neglected in the WKBJ analysis for large k_{z}. We apply the WKBJ approximation
to the fourthorder ODE (67), and we find and four solutions for S_{0} where
The WKBJ solution with is
where A_{3} and A_{4} are constants. This solution implies that thus simply increases or decreases exponentially as y → ∞ without turning points. Therefore, an eigenfunction satisfying the decaying boundary conditions as y → ∞ cannot be constructed. Thus, we impose A_{3} = A_{4} = 0. The WKBJ solution with can be either evanescent or wavelike depending on the sign of where
On the one hand, the WKBJ solution is evanescent if :
where C_{1} and C_{2} are constants, and . On the other hand, the WKBJ solution is wavelike if :
where D_{1} and D_{2} are constants. In these expressions, we consider only the leading order term S_{0} and neglect higherorder terms.
As conducted in the previous subsection, to derive asymptotic dispersion relations, we need to perform analyses around turning points where . To construct an eigenfunction, we first note that is negative as y → ∞ since . We thus have the exponentially decaying WKBJ solution for as follows:
Around the turning point where , we expand , where is the derivative of at the turning point, , and is a small parameter satisfying . Using this expansion in Eq. (67), we obtain the local equation around :
whose solutions are the Airy functions , where c_{1} and d_{1} are constants. Matching the asymptotic behavior of the Airy functions as and of the WKBJ solution (74) as , we have the following WKBJ solution in the range between the two turning points:
For y → −∞, we consider the decaying WKBJ solution for as
From the local solution behavior around the turning point , we can match the two WKBJ solutions (76) and (77). It leads to the following dispersion relation in a quantized form:
where m_{0} is the positive integer for the mode number.
Furthermore, we apply the Taylor expansion for σ
to the quantized dispersion relation (78), and we get
and
Since σ_{1, 0} is positive, we consider hereafter the first mode with m_{0} = 1 for the asymptotic results of the most unstable mode. In Fig. 10, we display the growth rate σ as a function of the vertical wavenumber k_{z} for various values of f and at small Pe = 0.1 for k_{x} = 0, N = 1, and Re = ∞, and we compare the numerical results with the asymptotic predictions from Eq. (79). We see that they are in very good agreement, especially as k_{z} → ∞.
Fig. 10. Numerical results (lines) of the growth rate σ as a function of k_{z} for various parameter sets of : (1,1) (black, Ω = 0.707 and θ = 45°), (1,2) (blue, Ω = 1.118 and θ = 63.4°), (1,3) (red, Ω = 1.581 and θ = 71.6°), (−0.5, 1) (green, Ω = 0.559 and θ = 116.6°), (−1, 1) (yellow, Ω = 0.707 and θ = 135°), at Pe = 0.1, k_{x} = 0, N = 1, and Re = ∞, and asymptotic results (circles) from Eq. (79) in the limit Pe → 0. 
Figure 11 shows how the maximum growth rate σ_{0, 0} depends on the Coriolis parameters f and . It is very remarkable that for , σ_{0, 0} is positive for any value of f including negative (f < 0) and large (f > 1) values, which are the two ranges outside the inertially unstable regime in the traditional approximation. The maximum value of σ_{0, 0} at a fixed is always attained at f = f_{max} = 0.5 with the corresponding value σ_{max, 0} = 0.5.
Fig. 11. Contours of σ_{0, 0} in the parameter space of . 
Using the expressions of the Coriolis parameters f = 2Ω cos θ and , we can see how the maximum growth rate σ_{0, 0} depends on the rotation rate Ω and colatitude θ. Figure 12 shows the maximum growth rate σ_{0, 0} plotted against the absolute Coriolis parameter 2Ω for various colatitudes. At the northern pole (θ = 0°), the case considered in the traditional approximation, the inertial instability only exists in the range 0 < 2Ω < 1. On the other hand, in the northern hemisphere 0° < θ ≤ 90°, we see that the inertial instability exists in any values of 2Ω > 0, and the growth rate reaches its maximum σ_{max} = 0.5 at 2Ω_{max} = 0.5/cos θ. After the peak, σ_{max} decreases with 2Ω and asymptotes to a constant value as 2Ω → ∞. By considering the limit 2Ω → ∞, we can find the asymptote of σ_{max} as follows:
Fig. 12. Maximum growth rate σ_{max} = σ_{0, 0} of Eq. (80) as a function of the absolute Coriolis parameter 2Ω at various colatitudes θ. Solid lines and dashed lines denote the results in the northern and southern hemispheres, respectively. 
In the southern hemisphere 90° < θ < 180°, the growth rate always increases with 2Ω and it reaches the maximum sin θ/2 as 2Ω → ∞. At the southern pole (θ = 180°), we find no inertial instability.
As a summary, it is remarkable that the inertial instability for the highdiffusivity case (Pe → 0) exists in the ranges
or
This colatitude range 0° ≤θ < 180° is much wider than the range of Eq. (63), which belongs to the northern hemisphere in the nondiffusive case (Pe = ∞). Since the highthermaldiffusivity regime is relevant for stellar interiors, the inertial instability due to the horizontal shear can be an important source of turbulence in stars at any colatitude except at the southern pole. This unstable regime in the limit Pe → 0 is equivalent to the regime of the GSF instability that occurs when the rotation profile varies along the rotation axis for small Pe^{2}.
5. Detailed parametric investigation
The detailed WKBJ analysis performed in the previous section provided explicit expressions for the dispersion relation of the inertial instability in two cases: one with no thermal diffusion and the other with high thermal diffusivity, without the traditional approximation. However, we do not fully understand how the inertial and inflectional instabilities are modified in other regimes such as finite Pe, finite Re, etc. In this section, we present more broad and detailed numerical results to understand the parametric behavior of the inflectional and inertial instabilities on Pe, Re, N, f, and .
5.1. Maximum growth rate of the inflectional instability: effects of N, and Pe
In the traditional approximation (i.e., ), the inflectional instability always has the maximum growth rate σ_{max} ≃ 0.1897 in the inviscid limit Re = ∞ at k_{x} ≃ 0.445 and k_{z} = 0, independently from the Coriolis parameter f, the BruntVäisälä frequency N, and the Péclet number Pe (Deloncle et al. 2007; Arobone & Sarkar 2012; Park et al. 2020). But as shown in Figs. 4 and 5, the maximum growth rate of the inflectional instability depends on the values of N and Pe if . For instance, we plot in Fig. 13 the maximum growth rate σ_{max} of the inflectional instability over the parameter space (k_{x}, k_{z}) at f = 0 (i.e., on the equator) in the inviscid and nondiffusive limits (i.e., Re = ∞ and Pe = ∞). We see how the maximum growth rate σ_{max} is modified by the stratification for various values of . We found numerically that the maximum growth rate σ_{max} of the inflectional instability is attained for a finite k_{x} at k_{z} = 0 for the nondiffusive case Pe = ∞, and that it is independent of f as explained by Eq. (39). While the maximum growth rate σ_{max} = 0.1897 is invariant for , we found that σ_{max} at a fixed is reduced for weakly stratified fluids with N ≪ 0.7 while σ_{max} surpasses 0.1897 for stratified fluids with N ≳ 0.7. The maximum of σ_{max} over N, namely max(σ_{max}), is thus greater than 0.1897 and it increases with .
Fig. 13. Maximum growth rate σ_{max} of the inflectional instability as a function of the BruntVäisälä frequency N for different values of at Re = Pe = ∞. 
Figure 14 shows the maximum growth rate of the inflectional instability over N as a function of for various values of Pe at k_{z} = 0 and f = 0. We see that at a finite Pe, max(σ_{max}) remains constant around 0.1897 for small and then increases with . At Pe = 1, we see that the maximum remains constant in the range . The onset of the growthrate increase is delayed as Pe decreases; therefore, we can say that the inflectional instability is stabilized as the thermal diffusivity increases (i.e., as Pe decreases). This stabilization of the inflectional instability was also reported in the traditional approximation (see e.g., Park et al. 2020), thus we can say that the inflectional instability is always stabilized by the thermal diffusivity for . And we can expect that the inflectional instability will not sustain in stars due to their high thermal diffusivity. As a summary, we display in Table 1 the parametric dependence for the inflectional instability in the nontraditional case (for other parameters f and N in the traditional case, refer to Deloncle et al. 2007; Arobone & Sarkar 2012; Park et al. 2020).
Fig. 14. Maximum of σ_{max} of the inflectional instability over N as a function of for different values of Pe (solid lines) at k_{z} = 0 and f = 0. 
Summary table for the inflectionalinstability growth rate by its variation with parameters , Pe, and Re
5.2. Maximum growth rate of the inertial instability at finite Pe
Figure 15 shows plots of the growth rate as a function of k_{z} for different values of Pe at N = 1 and k_{x} = 0. Two parameter sets are considered: (i.e., Ω = 1 and θ = 90°) and (i.e., and θ = 135°). These Coriolis parameters belong to the inertially stable regime in the nondiffusive limit Pe → ∞ according to the criterion (61). However, we see that it becomes unstable as Pe becomes finite, and the growth rate increases as Pe decreases. For a given parameter set of , we see that all growthrate curves for any Pe reach the same maximum value σ_{0, 0} given by Eq. (80) derived in the limit Pe → 0 as k_{z} → ∞. We also verified numerically that this asymptotic behavior is the same for other values of N at finite Pe. Therefore, the largest growth rate of the inertial instability is always σ_{0, 0}, and we have analogously the same relation between the maximum growth rate of the inertial instability and the Coriolis parameters at finite Pe.
Fig. 15. Growth rate σ_{r} of the inertial instability as a function of k_{z} for values of Pe at N = 1 and k_{x} = 0 for (i.e., Ω = 1 and θ = 90°; solid lines) and for (i.e., Ω = 0.707 and θ = 135°; dashed lines). Dotted lines indicate asymptotes σ_{0, 0} of Eq. (80) predicted from the WKBJ analysis in the limit Pe → 0. 
In Table 2, we summarize the variation of the growth rate of the inertial instability with parameters and Pe, and the unstable regimes in the two limits Pe → ∞ and Pe → 0. The destabilization of the inertial instability by Pe was already reported previously (see e.g., Park et al. 2020), and we confirm in this paper that this destabilization holds beyond the traditional approximation.
Summary table for the variation of the inertialinstability growth rate and the inertially unstable regime
5.3. Effect of the viscosity at finite Re
Now, we investigate how the viscosity modifies the inflectional and inertial instabilities. Figure 16 shows the growth rate of the inflectional and inertial instabilities for various parameters at different Reynolds numbers Re. For viscous cases, we fix the Prandtl number Pr defined as Pr = Pe/Re. We clearly see that the viscosity has a stabilizing effect on both instabilities. For the inflectional instability in Fig. 16a, the inviscid growth rate curve is increased for (black solid line) compared to that of the traditional case at (gray dashed line), while the viscous growth rate at decreases as Re decreases. It is remarkable that the inflectional instability persists at low Reynolds number Re = 10, and that the curve at has the same order of magnitude as the inviscid growth rate at . Therefore, we can conclude the inflectional instability is stabilized by the viscosity while it is destabilized as increases. This trend is summarized in Table 1.
Fig. 16. (a) Growth rate σ of the inflectional instability as a function of the streamwise wavenumber k_{x} for k_{z} = 0, N = 1, Pr = 1, , f = 0 (i.e., Ω = 1 and θ = 90°) for Re = ∞ (black), Re = 50 (blue), and Re = 10 (red). The gray dashdot line denotes the inviscid growth rate of the inflectional instability under the traditional approximation (i.e., ). (b,c) Growth rate σ of the inertial instability as a function of the vertical wavenumber k_{z} for k_{x} = 0, N = 1, , and f = 0.5 (i.e., Ω = 2.02 and θ = 82.9°) for (b) Pr = 1 and Re = ∞ (black), Re = 1000 (blue), and Re = 500 (red), and (c) Pr = 10^{−6} and Re = 10^{4} (black), Re = 5000 (blue), and Re = 1000 (red). Solid lines are numerical results, and dashed lines in (b) and (c) are asymptotic predictions from Eqs. (86) and (89), respectively. 
Figures 16b and c show the growth rate of the inertial instability at k_{x} = 0, N = 1, f = 0.5, and for two Prandtl numbers: Pr = 1 and Pr = 10^{−6}. For both cases, the maximum growth rate is now attained at finite k_{z} due to a stronger stabilization by the viscosity at large wavenumber k_{z}. For large Reynolds numbers and a unity Prandtl number (Pr = 1), we can apply the multiplescale analysis to perturbation equations and propose the viscous growth rate expressed in terms of the inviscid growth rate as
This expression of the growth rate has already been proposed in other works (see e.g., Arobone & Sarkar 2012; Yim et al. 2016). Using the asymptotic expression of the inviscid growth rate obtained from the WKBJ analysis in Sect. 4, we can express explicitly the viscous growth rate of the inertial instability at k_{x} = 0 as
We see in Fig. 16b that this expression for the viscous growth rate (86) is in good agreement with numerical results at finite Reynolds numbers and Pr = 1.
From (86), we can further derive the optimal vertical wavenumber k_{z, c}, at which the maximum growth rate is attained, by solving the equation ∂σ/∂k_{z} = 0:
Applying the optimal wavenumber k_{z, c} to the viscous growth rate (86), we find the critical Reynolds number Re_{c} above which the inertial instability occurs.
In Fig. 17a, we display the optimal wavenumber k_{z, c} (87) as a function of the Reynolds number Re and the corresponding critical Re_{c} at 2Ω = 1, N = 4, and k_{x} = 0 at various colatitudes θ. In Fig. 17b, we plot the critical Reynolds number Re_{c} (88) as a function of 2Ω for different colatitudes. Asymptotic predictions from Eq. (88) are of the same order of magnitude as numerical results for Re_{c}, and the small differences come from the small values of k_{z, c} predicted from the WKBJ analysis, which works well for large k_{z}. However, the interest of the asymptotic expression of Eq. (88) is that it can provide information about the critical Reynolds number Re_{c} in wide parameter spaces without exhaustive numerical computations.
Fig. 17. (a) Optimal wavenumbers k_{z, c} from (87) (solid lines) and k_{z, c0} from (90) (dashed lines) as a function of Re at 2Ω = 1, N = 4, and k_{x} = 0 for colatitudes θ: 30° (blue), 45° (red), and 60° (green). Filled circles denote the critical Reynolds numbers. (b,c) Critical Reynolds number Re_{c} as a function of the absolute Coriolis parameter 2Ω for N = 4, and (b) Pr = 1, (c) Pr = 10^{−6} at different colatitudes θ: 0° (black), 30° (blue), 45° (red), and 60° (green) from numerical results (crosses) and asymptotic results of Re_{c} (solid lines in b) from (88) and Re_{c, 0} (dashed lines in c) from (91). 
For very small Prandtl numbers, which is the relevant regime for the radiative zone of stars (Lignières 1999), we can similarly propose the viscous growth rate using the asymptotic growth rate (79) in the limit Pe = 0:
In Fig. 16c, we see that this viscous growth rate agrees very well with numerical results at finite Reynolds numbers and Pr = 10^{−6}. Moreover, we can further derive the optimal wavenumber k_{z, c0} and the critical Reynolds number Re_{c, 0} in the limit Pr → 0 as follows:
Figure 17a and c show the asymptotic predictions of k_{z, c0} and Re_{c, 0} and comparisons with numerical results at Pr = 10^{−6}. The critical Reynolds number is roughly within the same order of magnitude as numerical results, but the difference is large since the predicted critical wavenumber k_{z, c0} is very small, of order O(1), which lies much below the validity range of the WKBJ approximation. However, we can see that the critical Reynolds number Re_{c, 0} increases with 2Ω except at the northern pole θ = 0° where it is unstable in the range 0 < 2Ω < 1. We also found both numerically and theoretically that Re_{c, 0} does not change significantly with colatitude θ at a given Ω.
6. Effective horizontal turbulent viscosity
6.1. Turbulent viscosity induced by the inertial instability
As smallamplitude perturbations grow due to the instability mechanism, they first follow the linear instability growth, then reach a saturated equilibrium state and undergo a transition to turbulence as the perturbation amplitude increases and nonlinear effects become crucial. The nonlinear saturation of instabilities has been investigated in different astrophysical contexts and has been used to predict the order of magnitude of angular momentum transport, mixing, or dynamo actions (Spruit 2002; Denissenkov & Pinsonneault 2007; Zahn et al. 2007; Fuller et al. 2019).
In this section, we use our results on horizontal shear instabilities to derive an effective turbulent viscosity. We consider perturbations that are developed by the horizontal shear instabilities and reach a saturated state due to nonlinear interactions, which we model here through a turbulent viscosity. We adopt the approach described by Spruit (2002), Fuller et al. (2019) that allows the saturated state when the turbulent damping rate γ_{turb} equals the growth rate of the instability σ. This leads to the following expression of the local effective turbulent viscosity in the horizontal direction ν_{h, h}:
where is the horizontal wavenumber in the ydirection required for instability. We use here the notation ν_{h, h} first introduced by Mathis et al. (2018) that corresponds to the turbulent transport triggered in the horizontal direction by the instability of a latitudinal horizontal shear. We thus make the choice to consider the horizontal wavenumber in the latitudinal ydirection to characterize the horizontal effective turbulent viscosity.
We first consider the case of the effective turbulent viscosity ν_{h, h} induced by the inertial instability, whose dispersion relation has been derived using the WKBJ approximation in Sect. 4. Since the horizontal shear flow U(y) is inhomogeneous in the ydirection, the associated wavenumber k_{y} in the ydirection is not a constant but also a function of y. Moreover, the local wavenumber k_{y}(y) only exists when the solution is wavelike in the regime between the two turning points. For the nondiffusive case at Pe = ∞, we have the WKBJ expression for as
Since the exponent changes with y, we define the horizontal wavenumber as the average of the absolute part of the exponent over the two turning points as follows:
Similarly, we compute ν_{h, h} using the asymptotic dispersion relation (79) in the highly diffusive limit Pe → 0:
Figure 18 shows an example of as a function of in the two limits Pe = ∞ (solid line) and Pe → 0 (dashed line). The asymptotic growth rates σ from Eq. (58) for Pe = ∞ and from Eq. (79) for Pe → 0 derived with the WKBJ analyses are used to compute ν_{h, h}. We see that the effective turbulent viscosity reaches its maximum at a finite and decreases monotonically with . The monotonic decrease of ν_{h, h} is due to the fact that the growth rate reaches asymptotically the maximum growth rate σ_{max} as k_{z} → ∞ while the horizontal wavenumber , proportional to k_{z}, goes to infinity.
Fig. 18. Effective turbulent viscosity ν_{h} as a function of the horizontal wavenumber for f = 0.5, (i.e., Ω = 0.559 and θ = 63.4°), and k_{x} = 0 for Pe = ∞ and N = 1 (solid line), and as Pe → 0 (dashed line). 
Similarly to Fuller et al. (2019), where the minimum of the wavenumber is considered to get ν_{h, h}, we consider the maximum peak of ν_{h, h} as the representative value of the effective horizontal turbulent viscosity for given values of the physical parameters such as f and . In Fig. 19a and b, we plot contours of max(ν_{h, h}) in the parameter space of (N/2Ω, θ) for two different stratified fluids with N = 1 and N = 2 without thermal diffusion (i.e., Pe = ∞). For both cases, the turbulent viscosity max(ν_{h, h}) reaches its maximum around the colatitude θ ≃ 80° in the northern hemisphere close to the equator around the value N/2Ω ∼ 1. The effective turbulent viscosity in the southern hemisphere θ > 90° is zero since there is no instability for negative f = 2Ω cos θ. This implies that a strong turbulence and a large effective turbulent viscosity are expected near the equator in the northern hemisphere due to strong instability. In Fig. 19c, we plot max(ν_{h, h}) in the parameter space of (2Ω, θ) for highdiffusivity fluids in the limit Pe → 0 using Eq. (95). In this case, the effective turbulent viscosity is not zero in the southern hemisphere due to the presence of the inertial instability for negative f. The effective turbulent viscosity max(ν_{h, h}) has a maximum around θ ≃ 70°. This emphasizes that nontraditional effects with the positive horizontal Coriolis parameter play a crucial role on the effective turbulent viscosity, especially near the equator.
Fig. 19. (a,b) Contours of the maximum of the effective turbulent viscosity max(ν_{h, h}) in the parameter space (N/2Ω, θ) for Pe = ∞ and (a) N = 1, (b) N = 2. (c) Contours of the maximum of the effective turbulent viscosity max(ν_{h, h}) in the parameter space (2Ω, θ) as Pe → 0. 
From ν_{h, h}, the latitudinallyaveraged turbulent viscosity can be computed in two ways. On the one hand, we define the averaged turbulent viscosity following Zahn (1992) as
This average was used to compute the transport of angular momentum. On the other hand, we follow the definition by Mathis et al. (2018):
where is the averaged turbulent viscosity obtained in the way the mean transport of chemicals is computed. In Fig. 20, we plot these averaged viscosities. For the nondiffusive case (Pe = ∞), is larger than and they are reduced as N increases. The maxima of the two viscosities are reached around N/2Ω = 1 for both N = 1 and N = 2. For Pe → 0, the maximum is much lower than that for Pe = ∞, and it is reached around 2Ω = 1.
Fig. 20. Latitudinallyaveraged turbulent viscosities (black) and (blue) computed from contours in Fig. 19. 
6.2. Turbulent viscosity induced by the inflectional instability
The inflectional instability can also induce a turbulent transport in stellar radiation zones. Therefore, we also compute the effective turbulent viscosity based on the growth rate σ of the inflectional instability. The issues with the inflectional instability are that we do not have an analytic expression for the growth rate σ, and it is difficult to define systematically the horizontal wavenumber . One possible candidate for is to choose based on the solution’s asymptotic behavior as y → ∞. From numerical computations of the maximum growth rate of the inflectional instability in the parameter space (k_{x}, k_{z}), we display in Fig. 21 the effective turbulent viscosity in the parameter space of at Pe = ∞ and f = 0. It is important to note that the maximum growth rate σ_{max} of the inflectional instability is independent of f at Pe = ∞. We see that at a fixed N, max(ν_{h, h}) increases with for a weak stratification with N < 1 while it decreases with for a strong stratification with N > 1. The effective turbulent viscosity max(ν_{h, h}) is generally of order of the unity O(1) in the parameter space of and is smaller than the effective turbulent viscosity max(ν_{h, h}) induced by the inertial instability. We also verified numerically that the max(ν_{h, h}) of the inflectional instability decreases as the thermal diffusivity becomes finite and small. Therefore, we can conjecture that the effective turbulent viscosity induced by the inertial instability will be larger than that induced by the inflectional instability.
Fig. 21. Contours of the maximum of the effective turbulent viscosity max(ν_{h, h}) induced by the inflectional instability in the parameter space of for f = 0 (i.e., on the equator at θ = 90°) and Pe = ∞. 
7. Discussion with dimensional parameters
7.1. Inertial instability with latitudinal differential rotation
We discussed in the previous sections how the inertial instability is developed at a given latitude θ when the base shear of Eq. (4) is considered. In our analysis, we investigated the growth rate of the inertial instability in the dimensionless form by taking the velocity scale U_{0}, length scale L_{0}, and time scale t_{0} = L_{0}/U_{0}. We note that the length scale L_{0} is different from that of global largescale shear flows in stellar interiors as we consider smallscale shear flows and associated perturbations on a local plane. An important thing to note is that the base flow with positive shear S_{0} = U_{0}/L_{0} at y = 0 is used in the normalization. This implies that the velocity in the azimuthal direction always decreases as the colatitude θ increases. The latitudinal differential rotation in stars is, however, different from this case. For instance, stars can have a conical differential rotation in the simplest form as
where ϵ > 0 corresponds to the solarlike case, in which the equator rotates faster than the pole (e.g., ϵ ≃ 0.3 for the Sun), while ϵ < 0 corresponds to the antisolar case, where the equator rotates slower than the pole (Guenel et al. 2016). If instead of being constant, ϵ is proportional to r^{2}, one obtains the simplest form of cylindrical differential rotation (Baruteau & Rieutord 2013). At a constant radius, the latitudinal differential rotation follows the same law as in the conical case.
In the local frame rotating with Ω_{0} at the radius r = R and the colatitude θ, the relative azimuthal velocity U_{φ} is
(see also Hypolite et al. 2018). By considering the increment in the latitudinal direction on the local frame:
we obtain the base shear
The latitudinal shear S is negative (resp. positive) in the northern hemisphere and positive (resp. negative) in the southern hemisphere if ϵ > 0 (resp. ϵ < 0). Moreover, the latitudinal shear S is zero at the equator since the differential rotation of Eq. (98) is either at the maximum when ϵ > 0 or the minimum when ϵ < 0.
The growth rate of Eq. (79) that we derived from the WKBJ analysis is local but the growth rate induced by the latitudinal differential rotation of Eq. (98) needs to be obtained by the global stability analysis (see e.g., Guenel et al. 2016). However, we can approximately predict the growth rate σ_{θ} in the stellar radiation zones using the latitudinal shear S and the maximum growth rate of Eq. (80), which can be expressed in a dimensional form as follows:
Equation (102) can further be expanded in terms of the stellar rotation Ω_{0}, the colatitude θ, and ϵ as
Figure 22 shows the growth rate σ_{θ} versus the colatitude θ for various ϵ. The growth rate σ_{θ} increases with ϵ and reaches its maximum around θ ≃ 60° and 120°. It is zero at the poles and the equator. We see that, for the same ϵ, the growth rate of the antisolar case for ϵ < 0 is slightly larger than that of the solarlike case for ϵ > 0.
Fig. 22. Growth rate σ_{θ} as a function of the colatitude θ for various values of ϵ. Solid lines denote σ_{θ} with ϵ > 0 and dashed lines denote σ_{θ} with ϵ < 0. 
Furthermore, the latitudinallyaveraged^{3} growth rate can be computed as
In Fig. 23, we display as a function of ϵ for the solarlike and antisolar cases. In the range ϵ < 1, almost does not depend on the sign of ϵ. The growth rate shows a linear relation with ϵ as
Fig. 23. Latitudinallyaveraged growth rate as a function of ϵ for the solarlike case (ϵ > 0, solid line) and the antisolar case (ϵ < 0, dashed line). 
where a_{ϵ} ≃ 0.16. We can derive a similar scaling law for the growth rate in the limit Pe → ∞, more relevant to the solar tachocline case (Garaud 2020), if we use the growth rate (59) and a proper value of N.
7.2. Turbulent viscosities and characteristic time of turbulent transport
Using the averaged growth rate , we can define the turbulent viscosities induced by the horizontal shear due to the latitudinal differential rotation of Eq. (98). Consistently with the notations from Mathis et al. (2018), we define the turbulent viscosity in the latitudinal (horizontal) direction and the turbulent viscosity in the radial (vertical) direction . Following the assumption proposed by Spruit (2002), Fuller et al. (2019) that the shear instability grows and the momentum balances in nonlinear regime with the turbulent Reynolds stress, we find
where ∥ denotes the direction parallel to the stratification, ⊥ denotes the direction perpendicular to the stratification, k_{⊥} and k_{∥} are the characteristic wavenumbers in the horizontal and vertical directions, and l_{⊥} and l_{∥} are the length scales in the horizontal and vertical directions, respectively. We note here that we consider small length scales of the unstable modes that trigger turbulent transport, the length scales different from those of global largescale shear flows. The eddy viscosities are directly proportional to the horizontal shear as in the prescription derived by Mathis et al. (2004) (see their Eq. (18)). We note that it would also be possible to define diffusion coefficients for chemicals and ; for this one should compute .
While the turbulent viscosities depend on the direction with the scaling
(see also, Mathis et al. 2018), the dynamical time scale τ that characterizes the turbulence induced by the horizontal shear, and the transport of momentum and chemicals it triggers both in the vertical and latitudinal directions, is
which is simply the inverse of the growth rate .
Using the relation of Eq. (105), we can express the characteristic time τ in terms of Ω_{0} as
where τ_{0} denotes the rotation period of the star at the pole. We note that the characteristic time scale τ is similar to that proposed by Mathis et al. (2018) in the form τ = 1/S where S(r, θ) = r sin θ∂_{r}Ω characterizes the radial (vertical) shear.
According to the scaling (110), we estimate the transport time τ for the Sun at the level of the tachocline as 44.5 days when we use Ω_{0} ≃ 433 nHz and ϵ = 0.3 (García et al. 2007). For solarlike stars within the mass range 0.9 − 1.1 M_{⊙} (M_{⊙} is the Solar mass), if we make the rough assumption that the latitudinal differential rotation with ϵ = 0.3 is maintained during the evolution^{4}, we can roughly estimate that the transport time τ for the solarlike stars in the case of a slow initial rotation (Gallet & Bouvier 2015) is about 13.3 days, 9.7 days, and 16.8 days at the ages of 10 Myr, 100 Myr, and 1000 Myr, respectively. These time scales will be eventually used to compute the eddy viscosities in (109) upon the choice of the length scales.
For stars with a convective core, numerical simulations showed that the differential rotation in the core is mostly cylindrical (Browning et al. 2004; Augustson et al. 2016). For a 2 M_{⊙} star, Browning et al. (2004) found a rotation contrast within the core of between 30 and 60%. At the boundary between the convective core and the surrounding radiative zone, this is equivalent to ϵ ≃ 0.3 − 0.6. Using this estimate as well as rotation rates typical of Atype stars computed using observed surface velocities (see e.g., Zorec & Royer 2012) and MESA stellar structure and evolution models of a 2 M_{⊙} star with a Solar metallicity (Paxton et al. 2011) that provide us the stellar radius, the characteristic time τ ranges between 0.4 and 1 days during most of the main sequence and can go up to 1.5 days the end of the main sequence (around 1 Gyr).
If the turbulence triggered by the horizontal shear acts to damp its source, the horizontal differential rotation, as proposed by Zahn (1992) we thus predict a very efficient transport of angular momentum that can lead to the socalled shellular rotation where horizontal gradients of the angular velocity are weak. This also may be the origin of the observed very small radial extent of the Solar tachocline (Spiegel & Zahn 1992) and of an efficient mixing of chemicals in this region (Brun et al. 1999).
The behavior of the turbulence generated by horizontal shears has been recently explored in the nonrotating case by Cope et al. (2020) and Garaud (2020). They find when exploring the parameter space using Direct Numerical Simulations (DNS) a turbulent stratified regime where the turbulent transport can be modeled by an eddy diffusivity. This opens an interesting path to verify our predictions when such DNS will take rotation into account.
8. Conclusion
In this paper, we studied horizontal shear flow instabilities in stablystratified and thermallydiffusive fluids on a rotating plane where the full Coriolis acceleration with both vertical and horizontal components of the rotation vector is taken into account. For the canonical shear flow in the hyperbolic tangent profile, there exist two types of shear instabilities: the inflectional instability due to the presence of an inflection point and the inertial instability due to an inertial imbalance in the presence of the Coriolis acceleration. In the presence of positive horizontal Coriolis parameter , we found that both the inflectional and inertial instabilities are strongly affected. For instance, the inflectional instability, whose maximum growth rate is known to be independent of N, f, and Pe in the traditional approximation (Deloncle et al. 2007; Park et al. 2020), has now a maximum growth rate that depends on N and Pe. The horizontal Coriolis parameter destabilizes the inflectional instability for strong stratification while it stabilizes the instability at a small N. The thermal diffusivity at finite Pe also plays a stabilizing role in the inflectional instability. On the other hand, the inertial instability is destabilized by the thermal diffusivity, and the unstable regime is widely extended. For instance, in the nondiffusive limit Pe = ∞, the inertially unstable regime is found to be (i.e., tan^{−1}(N/2) < θ < 90°). More strikingly, it is inertially unstable for any f for high diffusivity fluids as Pe → 0 when (i.e., the inertial instability is active at any colatitude 0 < θ < 180° for the absolute Coriolis parameter 2Ω > 0). These unstable regimes as well as the dispersion relations for the inertial instability are derived analytically using the WKBJ approximation for large vertical wavenumber k_{z}, and this asymptotic analysis demonstrates a very good agreement with numerical results in the inviscid limit. Using the asymptotic expressions of the growth rate for the inertial instability, we also predicted the critical Reynolds number above which the growth rate becomes positive. Finally, we proposed prescriptions for the effective horizontal and vertical turbulent viscosities induced by the inertial and inflectional instabilities and found that the inertial instability plays an important role in the turbulent transport near the equator.
Observational and numerical studies suggest that stellar radiative zones have a mostly uniform rotation, whereas stellar convective zones are differentially rotating, but neutrally stratified. Therefore, one expects the present study to be relevant near the boundary between the convective and radiative zones. Because of the small values predicted for the time that characterizes the turbulence triggered by horizontal shear flows in such regions, this turbulent transport can have a strong impact on the structure and the evolution of stars, for instance by interacting with overshooting flows and by extracting angular momentum and chemical elements from the core to the envelope. In particular, the horizontal shear could be a crucial physical ingredient needed to explain the observed structure of the solar and stellar tachoclines (Spiegel & Zahn 1992; Brun et al. 1999). Other ingredients that are missing in the present work, such as magnetic fields (Gough & McIntyre 1998; Strugarek et al. 2011; AcevedoArreguin et al. 2013; Barnabé et al. 2017), may also play an important role.
Our predictions concerning the turbulent viscosities need to be confirmed by fully turbulent numerical simulations. In particular, global numerical simulations would allow us to validate the local approach used in this work, and the relevance of the latitudinally averaged quantities derived for stellar evolution codes. Besides, the local results could be used to build subgrid models for largeeddy simulations and stellarevolution calculations to better capture smallscale transport processes.
We here use the latitudinal average defined by Zahn (1992) for quantities related to angular momentum. For instance, he defined the shellular rotation as .
Using scaling laws computed in Brun et al. (2017) and grids of stellar models that take rotation into account using the STAREVOL code (Amard et al. 2019), Astoul et al. (2020) demonstrates that 0.15 < ϵ < 0.4 that will not predict orders of magnitude for τ
Acknowledgments
The authors acknowledge support from the European Research Council through ERC grant SPIRE 647383 and from GOLF and PLATO CNES grants at the Department of Astrophysics at CEA ParisSaclay. We thank the referee, Prof. A. J. Barker, for his constructive comments that allow us to improve the article.
References
 Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (Dover Publications) [Google Scholar]
 AcevedoArreguin, L. A., Garaud, P., & Wood, T. S. 2013, MNRAS, 434, 720 [Google Scholar]
 Aerts, C., Mathis, S., & Rogers, T. M. 2019, ARA&A, 57, 35 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Amard, L., Palacios, A., Charbonnel, C., et al. 2019, A&A, 631, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Antkowiak, A. 2005, Ph.D. Thesis, Université Paul Sabatier de Toulouse [Google Scholar]
 Arobone, E., & Sarkar, S. 2012, J. Fluid Mech., 703, 29 [Google Scholar]
 Astoul, A., Park, J., Mathis, S., Baruteau, C., & Gallet, F. 2020, A&A, submitted [Google Scholar]
 Augustson, K. C., Brun, A. S., & Toomre, J. 2016, ApJ, 829, 92 [Google Scholar]
 Barker, A. J., Jones, C. A., & Tobias, S. M. 2019, MNRAS, 487, 1777 [Google Scholar]
 Barker, A. J., Jones, C. A., & Tobias, S. M. 2020, MNRAS, 495, 1468 [Google Scholar]
 Barnabé, R., Strugarek, A., Charbonneau, P., Brun, A. S., & Zahn, J.P. 2017, A&A, 601, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baruteau, C., & Rieutord, M. 2013, J. Fluid Mech., 719, 47 [Google Scholar]
 Belkacem, K., Marques, J. P., Goupil, M. J., et al. 2015a, A&A, 579, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Marques, J. P., Goupil, M. J., et al. 2015b, A&A, 579, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Browning, M. K., Brun, A. S., & Toomre, J. 2004, ApJ, 601, 512 [Google Scholar]
 Brun, A. S., & Browning, M. K. 2017, Liv. Rev. Sol. Phys., 14, 4 [Google Scholar]
 Brun, A. S., Strugarek, A., Varela, J., et al. 2017, ApJ, 836, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., TurckChièze, S., & Zahn, J. P. 1999, ApJ, 525, 1032 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, P., & MacGregor, K. B. 1993, ApJ, 417, 762 [Google Scholar]
 Cope, L., Garaud, P., & Caulfield, C. P. 2020, J. Fluid Mech., 903, A1 [Google Scholar]
 Deloncle, A., Chomaz, J.M., & Billant, P. 2007, J. Fluid Mech., 570, 297 [Google Scholar]
 Denissenkov, P. A., & Pinsonneault, M. 2007, ApJ, 655, 1157 [Google Scholar]
 Eckart, C. 1960, Hydrodynamics of Oceans and Atmospheres (Elsevier) [Google Scholar]
 Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fabre, D., & Jacquin, L. 2004, J. Fluid Mech., 500, 239 [Google Scholar]
 Fricke, K. 1968, Z. Astrophys., 68, 317 [NASA ADS] [Google Scholar]
 Fuller, J., Piro, A. L., & Jermyn, A. S. 2019, MNRAS, 485, 3661 [NASA ADS] [Google Scholar]
 Gagnier, D., & Garaud, P. 2018, ApJ, 862, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Gagnier, D., Rieutord, M., Charbonnel, C., Putigny, B., & Espinosa Lara, F. 2019, A&A, 625, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gallet, F., & Bouvier, J. 2015, A&A, 577, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gallet, F., Bolmont, E., Mathis, S., Charbonnel, C., & Amard, L. 2017, A&A, 604, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Garaud, P. 2020, ApJ, 901, 146 [Google Scholar]
 Garaud, P., Gagnier, D., & Verhoeven, J. 2017, ApJ, 837, 133 [NASA ADS] [CrossRef] [Google Scholar]
 García, R. A., TurckChièze, S., JiménezReyes, S. J., et al. 2007, Science, 316, 1591 [Google Scholar]
 Gerkema, T., & Shrira, V. I. 2005, J. Fluid Mech., 529, 195 [Google Scholar]
 Gerkema, T., Zimmerman, J. T. F., Maas, L. R. M., & van Haren, H. 2008, Rev. Geophys., 46, RG2004 [Google Scholar]
 Goldreich, P., & Schubert, G. 1967, ApJ, 150, 571 [Google Scholar]
 Gough, D. O., & McIntyre, M. E. 1998, Nature, 394, 755 [Google Scholar]
 Griffiths, S. D. 2008, J. Fluid Mech., 605, 115 [NASA ADS] [Google Scholar]
 Guenel, M., Baruteau, C., Mathis, S., & Rieutord, M. 2016, A&A, 589, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hirschi, R., Meynet, G., & Maeder, A. 2004, A&A, 425, 649 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Høiland, E. 1941, in Avhandgliger Norske VidenskapsAkademi i Oslo, I, Math.Naturv. Klasse, 11, 1 [Google Scholar]
 Hypolite, D., Mathis, S., & Rieutord, M. 2018, A&A, 610, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Knobloch, E., & Spruit, H. C. 1982, A&A, 113, 261 [NASA ADS] [Google Scholar]
 Kulenthirarajah, L., & Garaud, P. 2018, ApJ, 864, 107 [Google Scholar]
 Lignières, F. 1999, A&A, 348, 933 [NASA ADS] [Google Scholar]
 Maeder, A. 2003, A&A, 399, 263 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars. [Google Scholar]
 Maeder, A., & Zahn, J.P. 1998, A&A, 334, 1000 [NASA ADS] [Google Scholar]
 Maeder, A., Meynet, G., Lagarde, N., & Charbonnel, C. 2013, A&A, 553, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marques, J. P., Goupil, M. J., Lebreton, Y., et al. 2013, A&A, 549, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., & Zahn, J. P. 2004, A&A, 425, 229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., Palacios, A., & Zahn, J.P. 2004, A&A, 425, 243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., Neiner, C., & Tran Minh, N. 2014, A&A, 565, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., Prat, V., Amard, L., et al. 2018, A&A, 620, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matt, S. P., Brun, A. S., Baraffe, I., Bouvier, J., & Chabrier, G. 2015, ApJ, 799, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Moss, D. 1992, MNRAS, 257, 593 [NASA ADS] [Google Scholar]
 Park, J. 2012, Ph.D. Thesis, Ecole Polytechnique [Google Scholar]
 Park, J., & Billant, P. 2012, J. Fluid Mech., 707, 381 [Google Scholar]
 Park, J., & Billant, P. 2013, J. Fluid Mech., 725, 262 [Google Scholar]
 Park, J., Billant, P., & Baik, J.J. 2017, J. Fluid Mech., 822, 80 [Google Scholar]
 Park, J., Prat, V., & Mathis, S. 2020, A&A, 635, A133 [EDP Sciences] [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 1 [Google Scholar]
 Pinçon, C., Belkacem, K., Goupil, M. J., & Marques, J. P. 2017, A&A, 605, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prat, V., & Lignières, F. 2013, A&A, 551, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prat, V., & Lignières, F. 2014, A&A, 566, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prat, V., Guilet, J., Viallet, M., & Müller, E. 2016, A&A, 592, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Richard, D., & Zahn, J.P. 1999, A&A, 347, 734 [Google Scholar]
 Rogers, T. M. 2015, ApJ, 815, L30 [Google Scholar]
 Schmid, P., & Henningson, D. S. 2001, Stability and Transition in Shear Flows (New York: SpringerVerlag) [Google Scholar]
 Solberg, H. 1936, in Procès Verbaux Ass. Météor., UGGI, 6e Assemblée Générale, Edinburgh, Mém. et Disc., 2, 66 [Google Scholar]
 Spiegel, E. A., & Zahn, J. P. 1992, A&A, 265, 106 [NASA ADS] [Google Scholar]
 Spruit, H. C. 1999, A&A, 349, 189 [NASA ADS] [Google Scholar]
 Spruit, H. C. 2002, A&A, 381, 923 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Strugarek, A., Brun, A. S., & Zahn, J. P. 2011, A&A, 532, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Strugarek, A., Bolmont, E., Mathis, S., et al. 2017, ApJ, 847, L16 [Google Scholar]
 Talon, S., & Charbonnel, C. 2005, A&A, 440, 981 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Talon, S., & Zahn, J. P. 1997, A&A, 317, 749 [NASA ADS] [Google Scholar]
 UdDoula, A., Owocki, S. P., & Townsend, R. H. D. 2009, MNRAS, 392, 1022 [Google Scholar]
 Wang, P., McWilliams, J. C., & Ménesguen, C. 2014, J. Fluid Mech., 755, 397 [Google Scholar]
 Yim, E., Billant, P., & Ménesguen, C. 2016, J. Fluid Mech., 801, 508 [Google Scholar]
 Zahn, J. 1983, in SaasFee Advanced Course 13: Astrophysical Processes in Upper Main Sequence Stars, eds. A. N. Cox, S. Vauclair, & J. P. Zahn, 253 [Google Scholar]
 Zahn, J.P. 1992, A&A, 265, 115 [NASA ADS] [Google Scholar]
 Zahn, J.P., Brun, A. S., & Mathis, S. 2007, A&A, 474, 145 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zeitlin, V. 2018, Phys. Fluids, 30, 061701 [Google Scholar]
 Zorec, J., & Royer, F. 2012, A&A, 537, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
Summary table for the inflectionalinstability growth rate by its variation with parameters , Pe, and Re
Summary table for the variation of the inertialinstability growth rate and the inertially unstable regime
All Figures
Fig. 1. (a) Schematic of the radiative and convective zones, colored as yellow and orange, respectively, for the configuration of lowmass stars rotating with angular speed Ω_{0}. (b) Horizontal shear flow U(y) in a local nontraditional fplane at a colatitude θ in the radiative zone; z_{c} denotes the transition altitude between the radiative and convective zones. 

In the text 
Fig. 2. Eigenvalue spectra in the space (σ_{i}, σ_{r}) at Re = ∞, N = 1, Pe = 1, f = 0.5, and (i.e., Ω = 1.031 and θ = 76°) for (a) (k_{x}, k_{z}) = (0.25, 0), (b) (k_{x}, k_{z}) = (0, 6), and (c) (k_{x}, k_{z}) = (0.25, 6). Triangles denote the maximum growth rates. 

In the text 
Fig. 3. Eigenfunctions of the most unstable modes in Fig. 2 at Re = ∞, N = 1, Pe = 1, f = 0.5, and : (a) the mode shape and (b) velocity v in the physical space (x, y) at z = 0 for (k_{x}, k_{z}) = (0.25, 0), and σ = 0.0712, (c) the mode shape and (d) velocity v in the physical space (y, z) at x = 0 for (k_{x}, k_{z}) = (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 , respectively. 

In the text 
Fig. 4. Contours of the maximum growth rate max(σ_{r}) in the parameter space of (k_{x}, k_{z}) (solid lines) for Re = ∞ and N = 1 for (a) , (b) , and (c) . 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 fplane approximation where . 

In the text 
Fig. 5. Inviscid growth rate of the inflectional instability for various parameter sets of and N at f = 0 (i.e., θ = 90°), Pe = ∞, and k_{z} = 0. 

In the text 
Fig. 6. Inviscid growth rate of the inertial instability for various Péclet numbers at N = 1 and k_{x} = 0 for f = 0 and , i.e., (Ω, θ) = (1, 90° ). 

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

In the text 
Fig. 8. (a) Contours of the maximum growth rate σ_{0} from Eq. (59) in the parameter space for N = 2. Black dashed lines represent the upper limit of the instability from Eq. (61) for different N, and the white dashed line represents f_{max} from (64) where the maximum growth rate is attained. (b) Contours of σ_{0} from Eq. (59) in the parameter space at f = 1. (c) Contours of the firstorder term σ_{1} in the parameter space for N = 2 and the first branch at m = 1. 

In the text 
Fig. 9. Numerical results (lines) and asymptotic results (circles) from (58) for various values of and N at f = 1, k_{x} = 0, Re = ∞, and Pe = ∞. 

In the text 
Fig. 10. Numerical results (lines) of the growth rate σ as a function of k_{z} for various parameter sets of : (1,1) (black, Ω = 0.707 and θ = 45°), (1,2) (blue, Ω = 1.118 and θ = 63.4°), (1,3) (red, Ω = 1.581 and θ = 71.6°), (−0.5, 1) (green, Ω = 0.559 and θ = 116.6°), (−1, 1) (yellow, Ω = 0.707 and θ = 135°), at Pe = 0.1, k_{x} = 0, N = 1, and Re = ∞, and asymptotic results (circles) from Eq. (79) in the limit Pe → 0. 

In the text 
Fig. 11. Contours of σ_{0, 0} in the parameter space of . 

In the text 
Fig. 12. Maximum growth rate σ_{max} = σ_{0, 0} of Eq. (80) as a function of the absolute Coriolis parameter 2Ω at various colatitudes θ. Solid lines and dashed lines denote the results in the northern and southern hemispheres, respectively. 

In the text 
Fig. 13. Maximum growth rate σ_{max} of the inflectional instability as a function of the BruntVäisälä frequency N for different values of at Re = Pe = ∞. 

In the text 
Fig. 14. Maximum of σ_{max} of the inflectional instability over N as a function of for different values of Pe (solid lines) at k_{z} = 0 and f = 0. 

In the text 
Fig. 15. Growth rate σ_{r} of the inertial instability as a function of k_{z} for values of Pe at N = 1 and k_{x} = 0 for (i.e., Ω = 1 and θ = 90°; solid lines) and for (i.e., Ω = 0.707 and θ = 135°; dashed lines). Dotted lines indicate asymptotes σ_{0, 0} of Eq. (80) predicted from the WKBJ analysis in the limit Pe → 0. 

In the text 
Fig. 16. (a) Growth rate σ of the inflectional instability as a function of the streamwise wavenumber k_{x} for k_{z} = 0, N = 1, Pr = 1, , f = 0 (i.e., Ω = 1 and θ = 90°) for Re = ∞ (black), Re = 50 (blue), and Re = 10 (red). The gray dashdot line denotes the inviscid growth rate of the inflectional instability under the traditional approximation (i.e., ). (b,c) Growth rate σ of the inertial instability as a function of the vertical wavenumber k_{z} for k_{x} = 0, N = 1, , and f = 0.5 (i.e., Ω = 2.02 and θ = 82.9°) for (b) Pr = 1 and Re = ∞ (black), Re = 1000 (blue), and Re = 500 (red), and (c) Pr = 10^{−6} and Re = 10^{4} (black), Re = 5000 (blue), and Re = 1000 (red). Solid lines are numerical results, and dashed lines in (b) and (c) are asymptotic predictions from Eqs. (86) and (89), respectively. 

In the text 
Fig. 17. (a) Optimal wavenumbers k_{z, c} from (87) (solid lines) and k_{z, c0} from (90) (dashed lines) as a function of Re at 2Ω = 1, N = 4, and k_{x} = 0 for colatitudes θ: 30° (blue), 45° (red), and 60° (green). Filled circles denote the critical Reynolds numbers. (b,c) Critical Reynolds number Re_{c} as a function of the absolute Coriolis parameter 2Ω for N = 4, and (b) Pr = 1, (c) Pr = 10^{−6} at different colatitudes θ: 0° (black), 30° (blue), 45° (red), and 60° (green) from numerical results (crosses) and asymptotic results of Re_{c} (solid lines in b) from (88) and Re_{c, 0} (dashed lines in c) from (91). 

In the text 
Fig. 18. Effective turbulent viscosity ν_{h} as a function of the horizontal wavenumber for f = 0.5, (i.e., Ω = 0.559 and θ = 63.4°), and k_{x} = 0 for Pe = ∞ and N = 1 (solid line), and as Pe → 0 (dashed line). 

In the text 
Fig. 19. (a,b) Contours of the maximum of the effective turbulent viscosity max(ν_{h, h}) in the parameter space (N/2Ω, θ) for Pe = ∞ and (a) N = 1, (b) N = 2. (c) Contours of the maximum of the effective turbulent viscosity max(ν_{h, h}) in the parameter space (2Ω, θ) as Pe → 0. 

In the text 
Fig. 20. Latitudinallyaveraged turbulent viscosities (black) and (blue) computed from contours in Fig. 19. 

In the text 
Fig. 21. Contours of the maximum of the effective turbulent viscosity max(ν_{h, h}) induced by the inflectional instability in the parameter space of for f = 0 (i.e., on the equator at θ = 90°) and Pe = ∞. 

In the text 
Fig. 22. Growth rate σ_{θ} as a function of the colatitude θ for various values of ϵ. Solid lines denote σ_{θ} with ϵ > 0 and dashed lines denote σ_{θ} with ϵ < 0. 

In the text 
Fig. 23. Latitudinallyaveraged growth rate as a function of ϵ for the solarlike case (ϵ > 0, solid line) and the antisolar case (ϵ < 0, dashed line). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.