Horizontal shear instabilities in rotating stellar radiation zones: I. Inflectional and inertial instabilities and the effects of thermal diffusion

Rotational mixing, the key process in stellar evolution, transports angular momentum and chemical elements in stellar radiative zones. In the past two decades, an emphasis has been placed on the turbulent transport induced by the vertical shear instability. However, instabilities arising from horizontal shear and the strength of the anisotropic turbulent transport that they may trigger remain relatively unexplored. This paper investigates the combined effects of stable stratification, rotation, and thermal diffusion on the horizontal shear instabilities in the context of stellar radiative zones. The eigenvalue problem describing linear instabilities of a flow with a hyperbolic-tangent horizontal shear profile was solved numerically for a wide range of parameters. As a first step, we consider a polar $f$-plane where the gravity and rotation vector are aligned. Two types of instabilities are identified: the inflectional and inertial instabilities. The inflectional instability that arises from the inflection point is the most unstable when at a zero vertical wavenumber and a finite wavenumber in the streamwise direction along the imposed-flow direction. The three-dimensional inflectional instability is destabilized by stratification, while it is stabilized by thermal diffusion. The inertial instability is rotationally driven, and a WKBJ analysis reveals that its growth rate reaches the maximum $\sqrt{f(1-f)}$ in the inviscid limit as the vertical wavenumber goes to infinity, where $f$ is the dimensionless Coriolis parameter. The inertial instability for a finite vertical wavenumber is stabilized as the stratification increases, whereas it is destabilized by the thermal diffusion. Furthermore, we found a self-similarity in both the inflectional and inertial instabilities based on the rescaled parameter $PeN^2$ with the P\'{e}clet number $Pe$ and the Brunt-V\"{a}is\"{a}l\"{a} frequency $N$.


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.2R (García et al. 2007) and to weak differential rotation in other stars (e.g., Mosser et al. 2012;Deheuvels et al. 2012Deheuvels et al. , 2014Kurtz et al. 2014;Saio et al. 2015;Murphy et al. 2016;Spada et al. 2016;Van Reeth et al. 2016;Aerts et al. 2017;Van Reeth et al. 2018;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 1983Zahn , 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;. 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. , 2018. 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;Ceillier et al. 2013;Marques et al. 2013;Cantiello et al. 2014;Eggenberger et al. 2019). 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 & 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;Maeder & Meynet 1996;Maeder 1997;Kulenthirarajah & Garaud 2018). Their predictions are now being tested using direct numerical simulations (Prat & Lignières 2013Prat 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;). 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 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 (Fig. 1a,b) while taking into account their important thermal diffusion.
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; 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 0 which is the Coriolis parameter defined as f 0 = 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 0 ( f 0 − U ) becomes negative (Fig. 1d). This condition leads to the inertially-unstable range 0 < f 0 < max(U ) for the hyperbolic tangent shear flow whose the maximum growth rate is found in inviscid limit as f 0 (max(U ) − 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: 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: (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 κ 0 and viscosity ν 0 into account: 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 J. Park, V. Prat and S. Mathis: Horizontal shear instabilities in rotating stellar radiation zones 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 stablystratified 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.

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 Ω: where u = (u, v, w) is the velocity, p is the pressure, θ = T − T 0 is the temperature deviation from the reference temperature T 0 , f = (0, 0, 2Ω) is the Coriolis vector which is here antialigned with the gravity g = (0, 0, −g), ρ 0 is the reference density, κ 0 is the reference thermal diffusivity, α T is the thermal expansion coefficient, and ∇ 2 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 where U 0 and L 0 are the reference velocity and length scale, respectively. Such a base flow is balanced with the pressure gradient as whereP is the base pressure profile. We consider a base temperature profileΘ(z) which increases linearly with height z: where ∆Θ 0 is the base-temperature difference along the vertical distance ∆z.

Linearized stability equations
Subject to the base state, we consider perturbationsũ = u − U = (ũ,ṽ,w),p = p −P, andT = θ −Θ. Hereafter, we use dimensionless parameters by converting the equations (4-6) into a set of dimensionless equations with the length scale as L 0 , the velocity scale as U 0 , the time scale as L 0 /U 0 , the pressure scale as ρ 0 U 2 0 , and the temperature scale as (L 0 ∆Θ 0 )/∆z. For infinitesimally small amplitude perturbations, we obtain the following linearized perturbation equations where prime denotes the total derivative with respect to y, N is the dimensionless Brunt-Väisälä frequency where f is the dimensionless Coriolis parameter and Pe is the Péclet number 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 2 is similar to the Richardson number Ri defined as Ri = N 2 v /S 2 v where S v is the vertical shear and N v is the Brunt-Väisälä frequency used for the vertical shear instability study . 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: where i 2 = −1,û = (û,v,ŵ),p, andT 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 equations (10-14) with these normal modes, we obtain the following linear stability equations 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 equations (19-23) into a single ordinary differential equation (ODE) forT . Firstly, from the vertical momentum and temperature equations, we expressŵ andp as a function ofT : From the horizontal momentum equations (20-21), we can expressû andv in terms ofp: Applyingp of (25) toû andv, we have the velocity perturbation u = (û,v,ŵ) expressed byT and the continuity equation (19) becomes the single 4th-order ODE forT : where∇ 4 = (∇ 2 ) 2 , s = −iσ+k x U is the complex Doppler-shifted frequency, and Γ is the function defined as 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.

Numerical method
We solve the set of equations (19-23) numerically by considering an eigenvalue problem in a simplified matrix form: where A and B are operator matrices expressed as where The eigenvalue problem (30) is discretized in the y-direction using the rational Chebyshev functions with the mappingỹ = y/ 1 + y 2 projecting the Chebyshev domainỹ ∈ (−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)

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).
The triangle in (a) and the square in (b) denote the most unstable growth rates of the inflectional and inertial instabilities, respectively.
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 forv from the continuity equation (19) and momentum equations (20-21): (see also Deloncle et al. 2007). It is important to note that the equation (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 (Delon cle 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 instabil-ity 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. 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 σ max = f (1 − f ). In the next section, the reasons why there is an infinite number of branches and why they approach f (1 − f ) as k z → ∞ will be explained by means of the WKBJ approximation. Fig. 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 inertiallyunstable 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). Fig. 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 |v|. For Fig. 4, we showv 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 ofT since this is the quantity we will use for its asymptotic analysis in Sect. 4. The mode shapev of the inflectional instability decreases exponentially as y → ±∞, and if plotted in physical space, we see that the perturbationṽ(x, y) is slightly inclined against the direction of the shear (see also Arobone & Sarkar 2012). On the other hand, the mode shapeT of the inertial instability mode shows there exists a zero crossing for the absolute part ofT at y = 0 while decaying exponentially as y → ±∞. The temperature perturbationT 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.

WKBJ analysis for the inertial instability
The maximum growth rate of the inertial instability is proposed to be σ max = 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 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: In the limit Pe → ∞, we can rearrange equation (35) as 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 equation (36) becomes the second-order ODE. Applying the WKBJ approximation to (36) for large k z : we obtain In this paper, we only consider the inertial instability mode which has the zero frequency (i.e., σ i = 0 and σ 2 + N 2 > 0) thus the sign of Γ determines the behavior of the solutions. We get evanescent solutions if Γ < 0: where A 1 and A 2 are constants, or wavelike solutions if Γ > 0: where B 1 and B 2 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 σ ± : which inform us where Γ = −σ 2 + 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 σ ± = ± √ −Φ 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± .
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 Fig. 4(c,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 1 = 0 and A 2 = 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(σ + ) = 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+ : 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ỹ = (y − y t+ )/ and an approximation Γ(y) ∼ Γ t+ ỹ where Γ t+ is the derivative of Γ in y at y = y t+ , we obtain the following local equation . The solution of this local equation (43) can be expressed in terms of derivatives of the Airy functions:T (ỹ) = a 1 Ai (ỹ) + b 1 Bi (ỹ), where a 1 and b 1 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ỹ → +∞: and whenỹ → −∞: (Abramowitz & Stegun 1972), we obtain the matched WKBJ solution in the region y t− < y < y t+ : where Similarly, the evanescent solution in y < y t− that decays exponentially as y → −∞: matches with the following solution in y t− < y < y t+ after the local solution around y t− is considered: Matching the wavelike solutions (46) and (49): that is we obtain the dispersion relation in the form of a quantization formula: 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 (Fig. 3b,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: we also assume that the growth rate can be expanded as σ = σ 0 − σ 1 k −p z + O(k −2p z ) where p is a positive integer. Applying this growth-rate expansion to the dispersion relation (53), we find that and we obtain the first-order term σ 1 , which leads to a more explicit dispersion relation for very large k z : where Since σ 1 is a positive value, we see that the maximum growth rate σ max = σ 0 = 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 σ 0 . While the maximum growth rate σ max is independent of N, the term σ 1 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.
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.
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 Applying the WKBJ approximation (37), we obtain the following relations at leading order: and at first order: The terms S 0 and S 1 satisfy or The general WKBJ solution with S 2 0 (y) = 1 iŝ but it has no turning point and the two solutions are always exponentially decaying or increasing as y → ±∞. Therefore, we must impose D 3 = D 4 = 0 to construct an eigenfunction. With S 2 0 = −Γ/σ 2 , the WKBJ solution is wavelike if Γ > 0: or evanescent if Γ < 0: 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+ : Around the turning point y t+ , we apply a new scaled coordinatẽ y = (y − y t+ )/ and obtain the following local equation Here, we take = k 2 z (−Γ t+ )/σ 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ỹ → ±∞ (Abramowitz & Stegun 1972) similar to in the previous section, we obtain the WKBJ solution for y t− < y < y t+ T (y) = Γ 1/4 where By matching the solution (68) with the solution for y < y t− : we obtain the following dispersion relation in a quantization form where m 0 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 : where 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.

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 σ 1 in (57); therefore, the ratio γ 1 between the two terms at first order is always larger than unity for positive N if we consider the same branch m = m 0 : 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.

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. 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 2 in Fig. 9b. We see that the rescaled growth rate curves collapse well for high N. This selfsimilarity represented by the rescaled parameter PeN 2 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 Prat & Lignières 2013).
Following the small-Pe approximation (Lignières 1999), we introduce the rescaled temperature deviationθ: In the limit of small Pe, we obtain the set of linear stability equations at leading order: The equations (76-80) can be simplified into an eigenvalue problem forv andθ: where We now see that the dependence of the problem on Pe and N 2 can be studied with the single parameter PeN 2 in the small-Pe limit. As illustrated in Fig. 9b, the growth rate in the limit of small Pe computed from the eigenvalue problem (81) (blue line in Fig. 9b) agrees very well with collapsed growth-rate curves plotted against the parameter PeN 2 for large N.

Effects of Pe on the inertial instability
The WKBJ approximation provides explicit dispersion relations for the inertial instability. They show how the growth rate depends on the stratification in the limit Pe → ∞ and how the growth rate in stratified fluids becomes equivalent to that in unstratified fluids as Pe → 0. Nonetheless, it is imperative to investigate whether this argument is valid for any k x and finite Pe since the WKBJ analysis in this paper is only applied for k x = 0 in the two extreme limits: Pe → ∞ and Pe → 0. In the inertiallyunstable regime f = 0.5 where both the inflectional and inertial instabilities are present (Fig. 8c), we see that the growth rate in the range k z > 0.6 is increased while the growth rate in the range k z < 0.6 is decreased as Pe decreases to zero. For both cases, we clearly see that the growth-rate curves in diffusive fluids (solid lines) asymptote to the curve for unstratified case at N = 0 (dashed line) as Pe → 0. From this asymptotic behavior, we can verify that the growth rate at small wavenumber k z < 0.6 corresponds to the inflectional instability that is stabilized as Pe → 0, while the growth rate at large wavenumber k z > 0.6 corresponds to the inertial instability that is destabilized as Pe → 0. Picking up the growth rate of the inertial instability at (k x , k z ) = (0.445, 1), we show in Fig. 10 the effects of the Péclet number Pe on the inertial instability for different values of the ratio N/ f at f = 0.5. For weakly stratified cases with low N/ f < 1, the growth rate remains constant as σ r 0.239 for a wide range of Pe. On the other hand, the growth rate decreases as Pe increases and asymptotes to σ r 0.190 for high N/ f . Similar to Fig. 9a, the shape of growth-rate curves for high ratio N/ f > 10 resemble with different onsets of stabilization as Pe increases. As the rescaled parameter PeN 2 is applied for growth rate curves (Fig. 10b), we see that the curves for large values of N/ f collapse and match with a stability curve computed from the eigenvalue problem (81) under the small-Pe approximation (blue line in Fig. 10b), similarly to the purely inflectional instability case in Fig. 9b. Fig. 11 shows contours of the growth rate σ r in the parameter space (N/ f, Pe) at k x = 0, k z = 20, and f = 0.5, a typical parameter set for the inertial instability. We see that for a fixed N/ f , the inertial instability destabilizes as the thermal diffusivity increases (i.e., as Pe → 0), similarly to what is observed in Fig. 10. For a fixed Pe, the growth rate decreases as N/ f increases (i.e., the stratification stabilizes the system at a fixed rotation). This can be mathematically predicted from the case for nondiffusive fluids as Pe → ∞ by the asymptotic dispersion relation (56) as the term σ 1 increases with N (i.e., the growth rate decreases with N). We can further derive from (56) the critical value of the ratio N/ f in the limit Pe → ∞ where the growth rate σ becomes zero: We see in Fig. 11 that the asymptotic prediction (87) for the critical value of N/ f lies in the stable regime obtained from numerical results. Also, it is notable from the equation (87) that for a fixed f , the critical ratio N/ f increases with the vertical wavenumber k z . This implies that the characteristic vertical length scale λ z = 2π/k z decreases as the ratio N/ f increases, as observed in other stratified-rotating flows (see e.g., Caton et al. 2000).

Conclusion
This paper investigates the instabilities of horizontal shear flow in stably-stratified, rotating, and thermally-diffusive fluids corresponding to stellar radiative regions. On the one hand, the inflectional shear instability always exists for the horizontal shear flow in a hyperbolic tangent profile whose maximum growth rate σ max = 0.1897 is attained at k x = 0.445 and k z = 0 indepen- Table 1. Summary table for the growth rate and its variation with parameters Pe, N, f , and (k x , k z ) for the three-dimensional inflectional instability (i.e., k z ≥ 0) and the inertial instability. dently of the stratification, rotation, and thermal diffusion. For the three-dimensional inflectional instability for nonzero vertical wavenumber k z > 0, the stable stratification that is inhibited by thermal diffusion has a destabilizing action. This is the opposite of the case of the vertical shear instability, which is inhibited by the stable stratification but favored by thermal diffusion that again diminishes its effects (e.g., Zahn 1983Zahn , 1992. On the other hand, the inertial instability is present in the range of 0 < f < 1 and its maximum growth rate σ max = f (1 − f ) is reached as k z → ∞ in the inviscid limit for both nondiffusive and highdiffusivity fluids. The analysis on the inertial instability for the case k x = 0 and large k z has been elaborated further through the WKBJ approximation in two limits: Pe → ∞ and Pe → 0 (i.e., for low and high thermal diffusivity, respectively), and explicit expressions for asymptotic dispersion relations are provided in both limits. For Pe → ∞, the growth rate decreases with N thus the stratification stabilizes the inertial instability, but the maximum growth rate σ max at infinite k z remains as f (1 − f ) independently of N. In the limit Pe → 0, the growth rate is no longer dependent on the stratification and becomes identical to the inertial instability growth rate for the unstratified case at N = 0. Detailed numerical studies confirm in the general parameter space (k x , k z ) that both the inflectional and inertial instabilities in thermally-diffusive fluids asymptote to those of the unstratified case as the thermal diffusivity increases (i.e., Pe → 0). The selfsimilarity of the growth rate in stratified-rotating fluids is also found such that the instabilities depend on the parameter PeN 2 for small Pe and strong stratification N. As a summary, we describe in Table 1 the growth rate and its variation with parameters Pe, N, f , k x and k z . The present work brings to light two horizontal shear instabilities that probably occur in stellar radiative zones but that had not been considered thoroughly before in stellar astrophysics, especially for the case of inertial instability on a local f -plane with the effects of thermal diffusion. The particularity of the stellar regime is that the high thermal diffusivity can weaken the stabilizing effect of the stratification for the inertial instability thus to enhance its development, while the three-dimensional inflectional instability is suppressed by the high thermal diffusivity and the inflectional instability mode becomes dominantly two-dimensional. We first investigated the linear instabilities in the polar regions but their nonlinear saturation and the resulting anisotropic turbulent transport of angular momentum and chemicals in both the horizontal and vertical directions has to be quantified. To derive the associated prescriptions for stellar evolution models, it is imperative to study the effects of the complete Coriolis acceleration on the instabilities at any co-latitude using the so-called nontraditional f -plane approximation (Gerkema et al. 2008). This will be done in the next article of the series.