Issue 
A&A
Volume 647, March 2021



Article Number  A144  
Number of page(s)  25  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202039148  
Published online  24 March 2021 
The complex interplay between tidal inertial waves and zonal flows in differentially rotating stellar and planetary convective regions
I. Free waves
^{1}
Laboratoire AIM ParisSaclay, CEA/DRF  CNRS  Université Paris Diderot, IRFU/DAp Centre de Saclay,
91191
GifsurYvette, France
email: aurelie.astoul@cea.fr
^{2}
Department of Applied Mathematics, School of Mathematics, University of Leeds,
Leeds,
LS2 9JT, UK
^{3}
Fluid and Complex Systems Research Centre, Coventry University,
Coventry
CV1 5FB, UK
^{4}
IRAP,
Observatoire MidiPyrénées, Université de Toulouse,
14 avenue Edouard Belin,
31400
Toulouse, France
^{5}
Université Grenoble Alpes, CNRS, IPAG,
38000
Grenoble, France
Received:
10
August
2020
Accepted:
21
December
2020
Context. Quantifying tidal interactions in closein twobody systems is of prime interest since they have a crucial impact on the architecture and the rotational history of the bodies. Various studies have shown that the dissipation of tides in either body is very sensitive to its structure and to its dynamics. Furthermore, solarlike stars and giant gaseous planets in our Solar System experience differential rotation in their outer convective envelopes. In this respect, numerical simulations of tidal interactions in these objects have shown that the propagation and dissipation properties of tidally excited inertial waves can be strongly modified in the presence of differential rotation.
Aims. In particular, tidal inertial waves may strongly interact with zonal flows at the socalled corotation resonances, where the wave’s Dopplershifted frequency is cancelled out. The energy dissipation at such resonances could deeply modify the orbital and spin evolutions of tidally interacting systems. In this context, we aim to provide a deep physical understanding of the dynamics of tidal waves at corotation resonances in the presence of differential rotation profiles that are typical of lowmass stars and giant planets.
Methods. In this work, we have developed an analytical local model of an inclined shearing box that describes a small patch of the differentially rotating convective zone of a star or a planet. We investigate the propagation and the transmission of free inertial waves at corotation, and more generally at critical levels, which are singularities in the governing wave differential equation. Through the construction of an invariant called the wave action flux, we identify different regimes of wave transmission at critical levels, which are confirmed with a onedimensional threelayer numerical model.
Results. We find that inertial waves can be fully transmitted, strongly damped, or even amplified after crossing a critical level. The occurrence of these regimes depends on the assumed profile of differential rotation, on the nature as well as the latitude of the critical level, and on wave parameters such as the inertial frequency and the longitudinal and vertical wavenumbers. Waves can thus either deposit their action flux in the fluid when damped at critical levels, or they can extract action flux from the fluid when amplified at critical levels. Both situations can lead to significant angular momentum exchange between the tidally interacting bodies.
Key words: hydrodynamics / waves / planetstar interactions / stars: rotation / planets and satellites: interiors / planets and satellites: dynamical evolution and stability
© A. Astoul 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
Tidal interactions are known to drive the late evolution of shortperiod planetary systems, such as hot Jupiters orbiting around their host star and, in our Solar System, the satellites around Jupiter and Saturn (e.g. Ogilvie 2014; Mathis 2019). In particular, the dissipation of tides in the convective envelope of lowmass host stars and giant planets can modify the spin of the tidally perturbed body, the orbital period, and the spinorbit angle of the perturber (e.g. Hut 1980; Ford & Rasio 2006; Lai 2012; Bolmont & Mathis 2016; Damiani & Mathis 2018). Inertial waves, which are driven by tidal forcing and restored by the Coriolis acceleration, are an important source of tidal dissipation in stellar (Ogilvie & Lin 2007; Barker & Ogilvie 2009; Bolmont & Mathis 2016) and planetary convective zones (Ogilvie & Lin 2004), where the action of turbulent motions on tidal flows is most often modelled as an effective frictional force or a viscous force with an effective viscosity that is much larger than the molecular viscosity (e.g. Zahn 1966, 1977; Duguid et al. 2020). For coplanar and circular systems, inertial waves are excited so long as the companion orbits beyond half its corotation radius (the orbit where the host’s rotation frequency is equal to the mean motion). Lowmass stars from K to F spectral type and giant gaseous planets both harbour a convective envelope surrounding a radiative and a solid (or diluted) core, respectively (e.g. Kippenhahn et al. 2012; Debras & Chabrier 2019). In these objects, inertial waves then propagate in a spherical shell and do not form regular normal modes of oscillation as they do in spherical and ellipsoidal geometries (Greenspan 1969; Bryan 1889, respectively). In contrast, they can focus on limit cycles, also called attractors of characteristics (Maas & Lam 1995), that are confined within the convective envelope (see also Rieutord & Valdettaro 1997). With a nonzero viscosity, attractors take the form of shear layers where the tidal wave’s energy and angular momentum can be deposited via viscous dissipation (Rieutord et al. 2001). Furthermore, viscous dissipation across shear layers can be more important as viscosity is weaker, as demonstrated notably by Ogilvie & Lin (2004) and Auclair Desrotour et al. (2015). In that respect, the tidal dissipation of inertial waves can compete with the dissipation of gravitoinertial waves in the radiative core, or it can be greater by several orders of magnitude than the dissipation of equilibrium tidal flows in the convective zone (i.e. the nonwavelike fluid’s response; see, e.g. Ogilvie & Lin 2007). The dissipation of tidally forced waves can have a great impact on the orbital and rotational evolution of the system (AuclairDesrotour et al. 2014; Bolmont & Mathis 2016; Gallet et al. 2018; Benbakoura et al. 2019). Moreover, the dissipation of the stellar dynamical and equilibrium tides varies significantly along the evolution of the star and is highly dependent on stellar parameters, such as the mass, the angular velocity, and the metallicity of stars (Mathis 2015; Gallet et al. 2017; Bolmont et al. 2017). This makes it desirable to include all stellar processes on tidal interaction, in particular differential rotation.
The frequencyaveraged tidal dissipation is often used to quantify the response of a body subject to tidal perturbations (Ogilvie & Lin 2004; Jackson et al. 2008). Yet, the dissipation of a tidally forced inertial wave is strongly correlated with the presence of an attractor at a specific eigenfrequency of the spherical shell (see Ogilvie 2009; Rieutord & Valdettaro 2010). Tidal dissipation at a given frequency may then alter each orbital and spin element of the twobody systems differently, as postulated, for instance, by Lai (2012) to explain the survival of hot Jupiters with completely damped spinorbit angles; this idea was revisited by Damiani & Mathis (2018) with an improved treatment of dynamical tides in the convective region. Additionally, in the context of the Jupiter and Saturn moon systems, Fuller et al. (2016) and Luan et al. (2018) also investigated the dependence of tidal dissipation on frequency to explain the rapid outward migration of the moons through the resonant locking of tidally forced internal modes in the giant gaseous planets. This concept could, for example, explain the high dissipation observed in Saturn as derived from astrometric measurements at the frequency of Rhea (Lainey et al. 2017) and at the frequency of Titan (Lainey et al. 2020).
Furthermore, the fact that all layers in a star or a planet do not rotate at the same speed (i.e. differential rotation) is rarely taken into account in the determination of tidal dissipation. Yet, differential rotation seems ubiquitous in lowmass stars and giant gaseous planets. The Sun’s surface rotates in ~ 25 days at the equator versus ~35 days near the poles, and a latitudedependent rotational gradient has also been observed in the Sun’s convective envelope thanks to helioseismology (Schou et al. 1998; Thompson et al. 2003). Through asteroseismology, latitudinal shears have been found to be comparable to that of the Sun for Sun analogues (Bazot et al. 2019), and they can be even larger for solarlike stars (Benomar et al. 2018). Essentially, differential rotation in lowmass stars depends on the effective temperature (Barnes et al. 2005, 2017) and seems to be more important when the convective envelope is thinner. Solarlike rotation profiles and antisolarlike rotation profiles (with faster poles and a slower equator) are expected in G and Ktype stars depending on their rotation rates, based on threedimensional numerical simulations (see in particular Brun et al. 2017; Beaudoin et al. 2018), while cylindrical rotation profiles are expected for fast rotators (Gastine et al. 2013). Regarding giant gaseous planets in our Solar System, the extent of zonal winds, which are visible on their surface as bands that run lengthwise, has been recently constrained by the probes Cassini and Juno. They extend to a depth of 3000 km for Jupiter (Kaspi et al. 2017) and penetrate down to 9000 km in Saturn (Galanti et al. 2019). Thus, the outermost molecular convective envelopes (Militzer et al. 2019; Debras & Chabrier 2019) are the seat of cylindrical differential rotation.
The study of the impact of differential rotation on the propagation and dissipation properties of inertial modes of oscillation began with the work of Baruteau & Rieutord (2013). They examined the impact of either a shellular (radial) or a cylindrical rotation profile on free inertial waves in an incompressible background by means of a WentzelKramersBrillouinJeffreys (WKBJ) linear analysis for an inviscid fluid and by solving the linearised hydrodynamics equations for a viscous fluid via a spectral code. Their linear analysis highlighted major differences relative to the case of solidbody rotation. Two regimes of propagation have been found in which inertial modes of oscillation can develop along curved paths of characteristics in the entire convective shell (which the authors named ‘D modes’) or in a restricted region of the convective shell, located between a turning surface and one of the shell’s boundaries (DT modes). Compared to solidbody rotation, the frequency range of the propagation of inertial modes is broader. Baruteau & Rieutord (2013) also pointed out the strong dissipation of wave energy at corotation resonances where the Dopplershifted wave frequency vanishes within the fluid. All these new properties have been retrieved by Guenel et al. (2016a), who in turn examined a conical (latitudinal) rotation profile, which is typical of lowmass (F to Ktype) stars. They also confirmed the existence of unstable inertial modes (i.e. modes with positive growth rates) at corotation resonances, which were only found for shellular rotation in Baruteau & Rieutord (2013). Tidal forcing of inertial waves with conical rotation was introduced by Guenel et al. (2016b) within a linear numerical exploration, which also highlighted the strong dissipation of inertial waves at corotation resonances, particularly at low viscosities. Favier et al. (2014) also studied tidally forced inertial waves, but through nonlinear numerical simulations. Differential rotation was triggered in their simulations by tidal waves depositing energy and angular momentum in an initially uniformly rotating spherical shell. In some cases, they observed hydrodynamical shear instabilities when the Ekman number (the ratio between the viscous and Coriolis accelerations) is sufficiently small.
Understanding how inertial waves interact with corotation resonances is thus a key issue in quantifying tidal dissipation, especially since waves may deeply interact with the background flow at this particular location, which in turn may alter the background flow (as was proposed first by Eliassen & Palm 1961, for terrestrial mountain waves). In binary systems and for latetype stars, Goldreich & Nicholson (1989) showed that the angular momentum transported by gravity waves and exchanged at corotation can lead to the successive synchronisation of the layers, from the base to the top of the radiative envelope. More generally, a body of work in various domains, from astrophysical disks (e.g. Goldreich & Tremaine 1979; Baruteau & Masset 2008; Latter & Balbus 2009; Tsang & Lai 2009) to geophysical fluid dynamics (e.g. Bretherton 1966; Yamanaka & Tanaka 1984), has tried to understand the properties of wave propagation and dissipation around corotation and, more generally, at all special locations in fluids that correspond to singularities in the linear wave propagation equation. We will refer to them as ‘critical levels’ in the following (Maslowe 1986), or as ‘critical layers’ in the case of a viscous medium. This distinction is analogous to the distinction between shear layers and attractors of characteristics that are kinds of singularities for the governing equation of inertial waves in a spherical shell. The aforementioned singularities can act very differently, with either severe absorption at the critical level (as in Booker & Bretherton 1967, for stratified vertical shear flows) or no attenuation if the wave propagates in a peculiar direction (Jones 1967; Acheson 1972; Grimshaw 1975a, for stratifiedvertical shear flows with rotation and magnetism). In other cases, a critical level may even give rise to wave amplification under certain conditions related to the first and second derivatives of the mean flow velocity (Lindzen & Tung 1978; Lindzen & Barker 1985, for barotropic and stratified shear flows, respectively). These studies all used an invariant quantity (the Reynolds stress or thewave action for rotating or magnetic flows) as a diagnostic tool to interpret the role of the critical level in terms of energy transmission and to quantify exchanges between the wave and the mean flow (Eliassen & Palm 1961; Bretherton 1966).
In light of these various studies, it is necessary to carefully consider corotation in differentially rotating convective zones. A local model can notably provide us a detailed understanding of physical processes at critical levels. While the propagation through a critical level of gravitoinertial waves in stratified shear flows and of Rossby waves in baroclinic and barotropic flows has been widely studied in the past decades, the behaviour of inertial waves in a latitudinal sheared flow with critical levels has so far been poorly investigated (e.g. Lindzen 1988, for a review). This is why we develop in this work a local Cartesian shearing box model to understand the complex interplay between tidal waves and zonal flows near critical levels. The concept of a shearing box for tidal flows was introduced by Ogilvie & Lesur (2012) to investigate the interactions between largescale tidal perturbations and convective motions. In our model, we focus on the latitudinal differential rotation of the meanflow, varying the box orientation to model either cylindrical or conical rotation. The behaviour of free inertial waves in this framework is then examined near critical levels using both analytical and numerical approaches.
This paper is organised as follows. In Sect. 2, we describe the local shear model with its main assumptions and the system of governing equations. In Sect. 3, we establish a secondorder ordinary differential equation (ODE) for the latitudinal perturbed velocity, and we derive the propagation properties of inertial waves for an inviscid fluid. This ODE is solved near each critical level for both conical and cylindrical rotation profiles, and we interpretenergy flux exchanges between the waves and the mean flow. We use, in Sect. 4, a threelayer numerical model to test our analytical predictions at critical levels. Frictional damping is included, and nonlinear mean flow profiles are also used. Astrophysical applications with implications for lowmass stars hosting close exoplanets and for giant gaseous planets in our Solar System are discussed in Sect. 5. In Sect. 6, we summarise the main results of the paper and discuss some perspectives and caveats.
Fig. 1 Sketch of the local Cartesian box in the convective region of a lowmass star or giant planet. Global spherical coordinates, including the depth r_{0}, the inclination of the box θ_{0}, and the colatitude θ of a point of interest M inside the box, are shown to facilitate the analogy between the spherical and the Cartesian geometries. 
Correspondence between local and global coordinate systems.
2 Local Cartesian model including differential rotation
2.1 Presentation of the model
The local model takes the form of an inclined sheared box, centred at a point C of a convective shell, as illustrated in Fig. 1. The inclined box model has already been used by Auclair Desrotour et al. (2015) to analytically characterise the properties of tidal gravitoinertial waves in the presence of viscous and thermal diffusion in stably stratified or convective regions, as well as by André et al. (2017) in layered semiconvective regions in giant planet interiors (see also Jouve & Ogilvie 2014, for twodimensional numerical simulations of inertial wave attractors). The local coordinate system (x, y, z) corresponds to the local azimuthal, latitudinal, and radial directions of global spherical coordinates, respectively, as presented in Table 1. The mean flow velocity U is directed along the local azimuthal axis e_{x} (we neglect possible meridional flows), and differential rotation is embodied by a latitudinal shear ∂_{y} U. As the box is tilted by an angle θ_{0} relative to the rotation axis, the rotation vector in the local coordinate system is (1)
where Ω_{0} is the rotation frequency of the star at the pole and and f are the normalised horizontal and vertical Coriolis components, respectively. It should be noted that the inclusion of both of these components means that we have gone beyond the traditional fplane approximation (see also, Gerkema et al. 2008). Furthermore, we made several hypotheses to model wave propagation in a latitudinal shear flow. The buoyancy acceleration was kept in the fluid equations for the background flow. The effective gravity acceleration g also includes the centrifugal acceleration, and the fluid’s angular velocity is assumed to be small compared to the critical angular velocity , where G, M, and R are the gravitational constant, the mass, and the radius of the body, respectively. Thus, the geometry of the body is close to spherical. Furthermore, the vector g is supposed to be uniform and constant in the whole box. This requires that the typical length of the box L satisfies L ≪ H_{p}, where H_{p} = −p(dz∕dp) is the vertical pressure scale height, with p being the pressure. We can assume this because tidally excited waves are expected to have smallscale structures (Ogilvie & Lin 2004; Rieutord & Valdettaro 2010; André et al. 2017). Moreover, the dimensions of the box were chosen to be small compared to the depth of the convective envelope so as to remove curvature effects.
2.2 Mean flow profile
In global spherical geometry, the mean flow based on a conical rotation profile Ω(θ) is written (e.g.in Guenel et al. 2016a) as: (2)
where e_{φ} is the azimuthal unit vector and r and θ are the radius and colatitude, respectively. We introduce u_{0} = rsinθΩ_{0}e_{φ}, the mean flow at a point M inside the box (see Fig. 1), without differential rotation. We also use the shear contrast δΩ = Ω(θ) − Ω_{0} (i.e. the difference between the angular frequency at colatitude θ and at the pole). The shear contrast is positive for the Sun since the equator rotates faster than the pole, and negative for antisolarlike rotating stars. Using the notations of Fig. 1, the centre C of the box is located at a distance r_{0} sinθ_{0} from the rotation axis. Accordingly, the latitudinal coordinate of the point M in the local frame is (3)
It should be noted that the radial coordinate r of the point M in spherical geometry can be written as r = r_{0} + z. Nevertheless, we neglected vertical displacements in the expression of the local shear because we are interested in how the (onedimensional) horizontal shear affects the wave dynamics, contrary to many studies on differential rotation in stars that have focused on the vertical shear (e.g. Mathis et al. 2004, 2018; Decressin et al. 2009; Alvan et al. 2013). Since y∕r and thus θ_{0} − θ are small, we provide the correspondences in terms of mean flows and shears between the two geometries in Table 1.
As an example, the shear contrast from solidbody rotation used by Guenel et al. (2016a) was: (4)
where χ is the magnitude of the shear between the equator and the pole. Performing a secondorder Taylor expansion around a fixed colatitude θ_{0}, such that θ = θ_{0}− y∕r_{0}, and at a specified depth r_{0} inside the convective region, the local mean flow U can be recast as
We point out that the Taylor expansion must be pushed further at the pole θ_{0} = 0 (and at the pole θ_{0} = π with an opposite sign): (6)
Accordingly, we can approximate a conical shear as a linear mean flow at the first order when the box is tilted. We recall that conical shear has been observed in the solar convective zone and is expected in slowly and moderately rotating solarlike stars (we refer the reader to Sect. 5.1 for a detailed discussion; see also Brun et al. 2015; Beaudoin et al. 2018; Benomar et al. 2018; Bazot et al. 2019). When the box is at the pole, y becomes the distance from the rotation axis (hereafter the ‘axial distance’). Thus, the mean flow mimics a cylindrical differential rotation that can be modelled using a cubic yprofile given by Eq. (6). This rotation profile is found in Jupiter and Saturn, as well as in rapidly rotating stars, as demonstrated, for instance, by Gastine et al. (2013) and Brun et al. (2015).
2.3 System of equations
To derive the system of governing equations for tidal waves in the local reference frame, we made several hypotheses. Stratification terms, which usually drive the propagation of internal gravity waves, were kept for the sake of clarity and will be methodically kept or removed after applying the Boussinesq approximation and setting the equations for inertial waves. Moreover, we assume that the action of turbulence can be modelled as a Rayleigh friction term in the momentum equation with an effective frictional damping rate σ_{f}. This simplifies the analytical solution of the fluid equations compared to the usual modelling of turbulence as an effective viscous force (see in particular Ogilvie 2009). The momentum, continuity, and thermodynamic equations for tidal waves in a differentially rotating Cartesian framework are thus:
where u, p, ρ, and f denote the velocity, pressure, density, and volumetric tidal forcing, respectively. We have also introduced c_{s}, the sound speed, and is the total derivative operator.
All variables are then linearised at the first order: Zeroorder terms correspond to background equilibrium quantities, and firstorder terms represent the leading perturbation. The local velocity, density, and pressure are therefore written as: (10)
where u = (u, v, w) in the local Cartesian basis. We have introduced the dimensionless parameter ϵ: (11)
where we use 1∕(2Ω_{0}), a characteristic time scale, and L, a characteristic length scale of the mean flow. These notations are based on those of Grimshaw (1975a) and adapted to our model. In the following, we will work with dimensionless variables using the above scaling, including 2Ω_{0} L to scale velocity and ρ_{T}gL to scale pressure, with ρ_{T} the reference density. The dimensionless momentum equation of the mean flow is: (12)
with n the unit vector parallel to the rotation axis. Projecting Eq. (12) into Cartesian coordinates, one can derive: (13)
At the leading order in ϵ, one can recognise the hydrostatic balance and, at the first order, the geostrophic balance (the set is akin to the thermalwind equilibrium assumption, see e.g. Grimshaw 1975a; Yamanaka & Tanaka 1984). We underline that tending ϵ to zero is similar to assuming the Boussinesq approximation. Indeed, all density variations are neglected, except the ones involved in the buoyancy force. The dimensionless BruntVäisälä frequency is (14)
where we introduce the dimensionless number , which is small when filtering acoustic waves. Consequently, the curl of Eq. (12) gives (15)
where we neglect the secondorder terms in ϵ.
Now, we make several assumptions to treat the propagation of inertial waves. As the convective motions are essentially adiabatic, the convective zone can be assumed to be neutrally stratified to a first approximation. Hence, the BruntVäisälä frequency N is cancelled out in the third density relationship of Eq. (15). Moreover, we make the Boussinesq approximation, which means that we neglect terms in ϵ and F in the final set of perturbed equations. Thus, the dimensionless linearised momentum, continuity, and thermodynamic equations are ultimately:
where v is the latitudinal velocity perturbation. We emphasise that, although vertical stratification has been filtered in the limit where N goes to zero, a horizontal stratification term remains in Eq. (18). As a result, we consider the inertial waves propagating in the inclined shear box where the mean flow is maintained by the thermalwind balance.
2.4 Equilibrium state of the background flow
It is worthwhile discussing the choice of keeping buoyancy forces in the zeroorder momentum equation. Without gravitational forces, the momentum equation for mean dimensional variables is written as a geostrophic balance: (19)
This balance satisfies the TaylorProudman theorem (Rieutord 2015), namely the geostrophic flow is independent of the coordinate parallel to the rotation axis. When taking the xaxis (the only nonzero) projection of the curl of this equation, one obtains the following relationship: (20)
Without the vertical stratification embodied by the BruntVäisälä frequency, nor latitudinal stratification, the equilibrium of a ydependent mean flow is thus not ensured for an incompressible fluid. An alternative to conserve the equilibrium without stratification would be to consider a zdependence of the mean flow. Such a possibility is not considered in this paper since we are mainly interested in latitudinal mean flow profiles. Furthermore, in addition to maintaining differential rotation, the latitudinal stratification can allow for the construction of an invariant that is useful for studying energy transfer at critical levels: the wave action flux. This will be discussed further in Sect. 3. Lastly, since at the poles, the latitudinal stratification term will not appear in the perturbed fluid equations (as we can see from Eq. (18)).
3 Dynamics of inertial waves at critical levels: analytical predictions
In this section, we analytically investigate the behaviour of inertial waves at critical levels in a nondissipative fluid at various colatitudes. For this purpose, we consider perturbations q in the normal mode (21)
where ω is the complex inertial frequency, k_{x} and k_{z} are the real streamwise and vertical wavenumbers, respectively, and c.c. is the complex conjugate.
3.1 Wave propagation equation in the latitudinal direction
Using the modal form Eq. (21) for ρ, p, and u, we solved the set of hydrodynamic equations, Eqs. (16)–(18), for the latitudinal velocity v. Considering free inertial waves (i.e. without forcing terms), the set of perturbation equations can be recast into a single secondorder ODE for v: (22)
where the prime now denotes the derivative according to y, and A, B, and C are the coefficients that can be simplified without friction as follows: (23)
where is the absolute wavenumber in the direction perpendicular to the y direction and σ = ω − k_{x}U is the (dimensionless) Dopplershifted wave frequency. We refer the reader to Appendix A for the detailed ODE derivation with friction and tidal source terms. Equation (22) becomes singular when A = 0 or σ = 0, and these singular points are called critical levels (see e.g. Bretherton 1966; Grimshaw 1975a). The critical level where the Dopplershifted frequency equals zero (i.e. σ = 0) can be met when the mean flow matches the local phase velocity, and this is also known as ‘corotation resonance’ (e.g. in Goldreich & Nicholson 1989; Goldreich & Tremaine 1979; Ogilvie & Lin 2004). When the Coriolis acceleration is not taken into account (i.e. when treating internal gravity waves), the corotation resonance is the unique critical level (see e.g. Booker & Bretherton 1967). At colatitudes other than the poles, the critical levels come in three flavours, the corotation σ = 0 and two other critical levels that are defined, in our model, by (we recall that is the latitudinal component of the rotation vector). These critical levels were similarly reported for vertical shear flows, as in the studies of Jones (1967) and Grimshaw (1975a) for vertical and inclined rotation vectors, respectively. In these works, the Dopplershifted frequency at critical levels other than the corotation resonance equals ± 2Ω_{v}, where Ω_{v} is the vertical component of the rotation vector.
3.2 Propagation properties
3.2.1 Dispersion relation, group, and phase velocities
The threedimensional dispersion relation is a fourthorder equation in Dopplershifted frequency when injecting wavelike solutions into the three directions, x, y, and z, in Eq. (22). In order to understand the main properties of waves at the critical level, we made the shortwavelength approximation in the meridional plane, as in Baruteau & Rieutord (2013) and Guenel et al. (2016a). This involves keeping only the secondorder derivatives in the y and z directions, and it reduces the relation dispersion to a secondorder equation when injecting plane wavelike solutions. In the local meridional plane, the differential equation reduces to a Poincarélike equation: (24)
where we recover the Poincaré equation (for the propagation of inertial waves in the inviscid limit; Cartan 1922) in the meridional plane when there is no shear (U′ = 0) and at the poles ( and f = 1). Moreover, we set v ∝ exp−i(k_{z}z − k_{x}x) so as to write the wave dispersion relation for the Dopplershifted frequency σ: (25)
where is the norm of the wave vector in the meridional plane (e.g. for fixed k_{x}), as in Baruteau & Rieutord (2013). Compared to solidbody rotation (see e.g. Rieutord 2015), an additional term () is present, which accounts for the latitudinal shear. Assuming that σ^{2} takes positive values (as in Baruteau & Rieutord 2013; Guenel et al. 2016a), we therefore introduce (26)
We can then specify the phase velocity in the meridional plane: (27)
In the same way, we can derive the expression for the group velocity in the meridional plane: (28)
We note that without differential rotation, the group velocity reduces to its wellknown expression for solidbody rotation (e.g. see Rieutord 2015): (29)
Moreover, as in solidbody rotation, the group velocity (Eq. (28)) and the phase velocity (Eq. (27)) lie in perpendicular planes: v_{g} ⋅v_{ϕ} = 0.
When the box is located at the north pole (θ_{0} = 0 in Fig. 1), by setting κ = 1 − U′, we recover (30)
as in Latter & Balbus (2009) and Baruteau & Rieutord (2013), where κ can be identified with the epicyclic frequency of Baruteau & Rieutord (2013) and k_{y} corresponds to the cylindrical component of the wavenumber (k_{s}).
3.2.2 Phase and group velocity at singularities
In this section, we derive the conditions required to meet singularities in terms of wavenumbers and shear, and we examine the implications for the phase and group velocities. When the box is inclined, for σ → 0 we must have γ → 0, meaning that v_{ϕ} →0 while v_{g} ⋅e_{y}→∞ and v_{g} ⋅e_{z}→∞ (possibility 1).
Guenel et al. (2016a) found similar results by studying the propagation of free inertial waves in a global frame with conical shear, namely, when their parameter (which is homogeneous to a frequency and equivalent to our γ parameter) goes to zero, the group velocity goes to infinity while the phase velocity is cancelled out. According to their work, an inertial wave may propagate across the corotation.
Now to get , we either need k_{y}→∞ at fixed k_{z}, which implies v_{ϕ} →0 and v_{g} →0 and means that inertial waves cannot get through the critical level (possibility 2), or k_{z}→ 0 at fixed k_{y}, which gives v_{ϕ} ⋅e_{z} → 0 and v_{g} ⋅e_{y} → 0, while and v_{g} ⋅e_{z}→ f∕k_{y}; the wave may then cross the critical level with some preferential direction (possibility 3).
Again, these conditions share some similarities with those observed for corotation in a global spherical geometry. The first of the three possibilities is analogous to the global phase and group velocities tending to zero when k_{s} →∞, with k_{s} the axial wavenumber in cylindrical coordinates (Baruteau & Rieutord 2013; Guenel et al. 2016a). This makes sense since here the axial distance is s = rsinθ and y ~ r_{0}(θ_{0} − θ). However, the second possibility is slightly different from both these previous works, in that k_{z}  → 0 at fixed k_{s}, where k_{z} is the global vertical wavenumber along the rotation axis, unlike our local vertical wavenumber k_{z} which is along the spherical radial coordinate.
We point out that the singularities at arise in our model because the rotation vector is inclined with respect to the local vertical axis of the box. In the global model of Guenel et al. (2016a), three conditions for a wave to meet the corotation exist, and these conditions are actually quite similar to the three above possibilities for waves in our model to interact either with the corotation σ = 0 or with the other critical levels at . Hence, the local critical levels at behave somewhat similarly to the corotation in the global framework, as if we partially broke the degeneracy in the local framework of the origin of the corotation found in the global framework.
When the box is at the north pole, the three possibilities to meet corotation are similar but lead to different relationships for the phase and group velocities: (i) κ → 0 (i.e. U′→ 1), meaning that v_{ϕ} →0 and v_{g} →0: The wave is totally absorbed at corotation. (ii) k_{y}→∞ at fixed k_{z}, which implies v_{ϕ} →0 and v_{g} →0: We can draw the same conclusion as in the previous case and this case is also analogous to the second possibility when the box is tilted. (iii) k_{z}  → 0 at fixed k_{y}, which gives v_{ϕ} →0, v_{g} ⋅ e_{y} → 0 while v_{g} ⋅e_{z}→ κ∕k_{y}: the wave energy does not cross the corotation in the latitudinal direction (equivalent to the vertical paths of characteristic in global cylindrical geometry, as in Baruteau & Rieutord 2013). At the north pole^{1}, we actually have a perfect match with the conditions given by Baruteau & Rieutord (2013) when using a cylindrical rotation profile for the mean flow.
3.2.3 Energetical aspects
In this section, we examine the energetic balance associated with inertial waves in our inclined shear box model, without assuming the shortwavelength approximation. This energetic balance does not include potential energy because of the adiabaticity of the convective region, but two additional terms relative to the solidbody rotation case appear, which come from the differential rotation. We denote the displacements along the vertical and latitudinal directions with η and ζ, respectively. Considering that ω = k_{x}c, where c is the longitudinal phase velocity (e.g. as in Booker & Bretherton 1967), we can use the firstorder definition (31)
This first allows us to express the perturbed density from Eq. (18) as (as a reminder, the symbol ′ has been dropped for perturbed quantities). Then, by multiplying the momentum equation, Eq. (16), by ρ_{0} u, we get the energy balance equation: (32)
where e_{k} = ρ_{0}u^{2}∕2 is the kinetic energy density and pu is the socalled acoustic flux. We now integrate the above energy balance equation over x and z, as well as over one wave period as the perturbed quantities have a wavelike form in these directions. Further assuming that the box is δ thick in the y direction, the energetic balance yields: (33)
where we introduce, from left to right, the power of the external pressures at the boundaries − δ∕2 and δ∕2 on the perturbed latitudinal flow, the work of the shear, the viscous dissipation, and the forcing power, which read, respectively: (34)
where the bar represents the average in the (z, x) plane over one period. We note that the energy density and the acoustic flux in the x and z directions drop out in Eq. (33) when integrating because of the wave periodicity in those directions. The quantity can also be seen as the power transferred from the mean flow to the perturbation (or conversely) by the Reynolds stress: (35)
where we use partial integration and the periodicity of perturbations in the x and z directions. At the pole, , so we recover the definition of the Reynolds stress in Miles (1961), who studied the stability of a twodimensional stratified ysheared flow (i.e. ). This quantity can also be called the latitudinal flux of horizontal momentum (in the (z, x) plane sense) in reference to the vertical flux of horizontal momentum in stratified zsheared flows. Moreover, we emphasise that the latitudinal flux of energy is not conserved, even in the inviscid freewave problem. This is due to the presence of the shear, as already stated for example by Eliassen & Palm (1961), who studied stratified vertically sheared flows. They underline that when the mean flow varies with height, the kinetic energy of the mean motion can be converted into wave energy. Without friction and forcing, the yderivative of the latitudinal flux is: (36)
Using the same method as Broad (1995), we multiplied the xprojection of the inviscid forcefree momentum equation by ζ: (37)
By multiplying by (U − c), the latitudinal flux of energy can thus be written as: (38)
By differentiating this relationship with respect to y, and by equalising with Eq. (36), one can obtain: (39)
that is with τ the Reynolds stress. Equation (39) is naturally satisfied at corotation, where U −c = 0, or if the Reynolds stress is uniform. Booker & Bretherton (1967) have shown that the Reynolds stress is discontinuous at a critical level, highlighting exchanges between wave energy and the mean flow. Compared to the analysis of Broad (1995) for threedimensional stratified shear flows, Eq. (39) is not vectorial, because our base flow is unidirectional.
3.2.4 Polarisation relations
For the forthcoming analysis, it is useful to derive expressions of the perturbed projected velocities and the perturbed reduced pressure^{2} Π = p∕ρ_{0}, namely the polarisation relations (see Appendix B for more details). In the inviscid freewave problem, these perturbed quantities can be written in terms of the latitudinal velocity, its derivative, and the shear: (40)
Without shear and at θ_{0} = 0, we recover the polarisation relations in the solidbody rotation case (see e.g. Rieutord 2015).
3.2.5 Conservation of the wave action flux
While the latitudinal flux of energy is not conserved in the whole domain, there is a conserved quantity, called the wave action flux, as introduced in Grimshaw (1975a): (41)
which is the latitudinal flux averaged over vertical and longitudinal wavelengths divided by the Dopplershifted frequency. A general treatment for the derivation of the wave action as a conserved quantity can similarly be found in Andrews & McIntyre (1978). The wave action flux is related to the Reynolds stress τ as . By using the expression for the perturbed reduced pressure Π derived in the previous section, the wave action flux now reads: (42)
Unlike the latitudinal flux of energy, but similarly to the Reynolds stress, this wave action flux is conserved along the latitudinal direction. One can demonstrate that in the whole domain except at critical levels by using the expression for the reduced pressure in Eq. (40) and the ODE (Eq. (22)). Several works have shown that a properly defined (i.e. conserved) angular momentum transport parameter can be found in zsheared mean flows without rotation (Booker & Bretherton 1967), with rotation under the traditional approximation (Jones 1967), and with rotation under the nontraditional approximation (Grimshaw 1975a). Verifying the conservation in the whole domain except at critical levels is really important because it brings to the fore energy transfers that are due to the critical levels. We specify that is a measure of wave energy through a surface (in the (z, x) plane) since is the energy density transported by the group velocity^{3} V_{g} in the latitudinal direction (e.g. Bretherton & Garrett 1968; Mathis & de Brye 2012). It should be underlined that the wave action flux has been defined in the inviscid limit and is not conserved when the friction is taken into account.
3.3 Inertial waves at critical levels when the box is tilted
In this section, we analytically investigate waves passing through the various critical levels in the tilted box. We examine the behaviour of the waves around the corotation σ = 0 and the critical levels when the box is tilted (for the corotation when the box is at the pole, see Sect. 3.4).
3.3.1 Critical levels at
In this subsection, we treat both singularities, , simultaneously. Although Eq. (22) does not have analytical solutions in general, it is still possible to study the behaviour of an inertial wave close to the critical levels defined by by approximating the ODE through its firstorder Taylor expansion in the vicinity of these singularities, and then by applying the Frobenius method. We introduced y_{±}, the location of the related critical level . For a linear mean flow profile U = Λy, with Λ a constant, y_{±} are given by: (43)
Without any assumption on the mean flow profile, the firstorder Taylor expansion of the ODE (Eq. (22)) near y_{±} is: (44)
where the symbol ± refers to the regular singularities^{4} y_{+} and y_{−} and is U′ evaluated at these singularities. The Frobenius method consists in injecting the power function into Eq. (44), with λ a constant to be determined (see e.g. Morse & Feshbach 1953). The corresponding indicial equation is then (46)
Therefore, the two independent solutions of Eq. (44) can be written as follows: (48)
where a_{n} and b_{n} are complex constants. Both solutions are valid in the vicinity of the critical level around which they are built in the complex plane, up to the next singularity if it exists. The coefficients a_{0} and b_{0} are unconstrained and depend on boundary conditions, unlike the other factors that can be determined by injecting the solutions from Eq. (48) into the linearised ODE (Eq. (22)) around y_{±} at the right order for the desired coefficients.
Near the critical points y_{±}, the total solution v_{±} is well approximated by the lowest orders of v_{1,±} and v_{2,±}: (49)
Owing to the existence of a branch point at y_{±} (since λ_{±} is complex),reconnecting solutions on either side of the critical levels is not straightforward. This requires both physical and mathematical arguments (see in particular Miles 1961; Booker & Bretherton 1967; Ringot 1998). In order to remove the degeneracy of the path from positive to negative y − y_{±} (i.e. choosing either e^{+iπ} or e^{−iπ}), we made use of a complex inertial frequency ω = ω_{R} + iω_{I}, assuming the radiation condition ω_{I} > 0. This condition ensures a nongrowing wave towards infinity. The Taylor expansion of the base flow at the first order in y − y_{±} gives (50)
and, by definition, we have (51)
Consequently, the solution below the critical level is unambiguous in terms of the above solution coefficients and depends on (52)
In other words, when taking y − y_{±} to decrease from positive to negative values, its complex argument changes continuously from 0 to . Thus, the appropriate path for determining the branch of passes under(above) y_{±} as long as () (the same reasoning can be found in Grimshaw 1975a). Therefore, the solution on both sides of the critical level y_{±} is: (53)
The remaining issue is now to know in which direction the wave is propagating. The second part of the solution can be assimilated to a wavelike solution with the varying latitudinal wavenumber . Moreover, according to Eq. (42), the wave action flux on either side of y_{±} is: (54)
The group velocity gives the direction towards which the energy is transported (we recall that , with V_{g} the group velocity and the local energy density). By consequence, for the solution featuring the coefficient b_{0}. If k_{z}f is positive, this wave transports energy downwards (upwards) across the critical level (). If k_{z} f is negative, the wave transports energy upwards (downwards) across the critical level (). In all cases, the action flux of the wave with the amplitude b_{0} will be transmitted (in the direction given by the sign of k_{z}f and the critical level y_{+} or y_{−}) by a factor , where (55)
with α_{k} = k_{z}∕k_{x} and , after passing through the critical level. Such a wave will always be attenuated since . The transmission factors and are displayed in Fig. 2 in terms of the absolute value of the shear Rossby number R_{o}  and the ratio of wavenumbers α_{k}. The lower the amplitude of the Rossby number and the lower the inclination, the more likely the wave is to be strongly attenuated at any α_{k}. We reiterate that a low Rossby number refers either to fastrotating stars or to low differential rotation. At the equator, one should note that f = 0, thus there is no transmission nor exchange of wave action flux near the critical levels y_{±} in the inviscid limit (see Eq. (54)). Results are the same for θ_{0} + kπ∕2 with k ∈ {0, 1, 2, 3} and for negative Rossby numbers. However, it has to be emphasised that the cases where the inclination satisfies θ_{0} = kπ with k ∈ {0, 1} are not well described by the attenuated factor and require a specific treatment, as discussed in Sect. 3.4.
It is important to note that, with fixed parameters {k_{z}, θ_{0}, R_{o}}, the attenuation of the wave action flux is specific to a single direction of wave propagation (i.e. the solution featuring the coefficient b_{0}). The solution of coefficienta_{0} is not affected by the attenuation. This is the socalled valve effect introduced by Acheson (1972) in the context of hydromagnetic waves in a rotating fluid. It was also evidenced by Grimshaw (1975a) and further discussed in Grimshaw (1979) for magnetogravitoinertial waves in an inviscid and compressible zsheared fluid.
Fig. 2 Transmission rate of a wave passing through any of the critical levels defined by as a functionof the absolute value of the shear Rossby number R_{o} = U′ (where U′ is scaled by 2Ω_{0}) and the ratioof the horizontal wavenumbers α_{k} = k_{z}∕k_{x} for a colatitude of the box θ_{0} = 10° (left panel) and θ_{0} = 80° (right panel). 
3.3.2 Inertial wave crossing corotation
We performed the same analysis as in the previous section to treat the corotation point y_{0} where σ = 0 (i.e. U(y_{0}) = U_{0} = ω∕k_{x}). The linearised ODE (Eq. (22)) near the corotation using the Taylor expansion of σ and U at the lowest orders is: (56)
where and are the first and second derivatives of the mean flow profile U evaluated at the critical level y_{0}. The singularity at the corotation is a regular singularity, and we can again use the Frobenius method. The indicial equation has solutions λ = {2, 1}. Since the difference between the two values of the exponent λ is an integer, one expects a second independent solution v_{2} of Eq. (56) that includes a logarithmic part, such as (e.g. Schmid et al. 2002) (57)
with the first solution and a_{n}, b_{n}, and L complex coefficients. However, when injecting v_{1} + v_{2} into Eq. (56), one finds L = 0, meaning that a sole polynomial solution in the form (58)
includes all the solutions of Eq. (56), with c_{0} = b_{0} and a_{1} = b_{1} + a_{0} determined by boundary conditions, and determined by recurrence via the expansion of Eq. (22) around y_{0}. As a result, thewave action flux given by Eq. (42) becomes (59)
just below and above the corotation, and it is continuous there, similarly as in Grimshaw (1975a), but here without being restricted to a linear mean flow profile. Hence, no transfer of wave action flux is expected at corotation in the inviscid limit when the box is inclined relative to the rotation axis (i.e. for conical differential rotation), regardless of the mean flow profile. This result also holds true when the box is located at the equator.
As in the works of Grimshaw (1975a) and Jones (1967), it is tempting to investigate the asymptotic behaviour of a wave when y →∞ in order to better constrain the propagation of waves through one or multiple critical levels. Nevertheless, the term in the ODE (Eq. (22)), which cannot be overlooked as it was in the aforementioned studies since we do not have vertical stratification, makes the singularity y = ∞ an essential (or irregular) singularity, and the Frobenius method cannot be applied. This term also prevents us from applying an analysis such as the WKBJ approximation because, even far from critical levels, the coefficients C of the ODE (in Eq. (22)) still have a strong dependence on the latitudinal coordinate when the box is tilted.
3.4 Inertial waves when the box is at the poles
When the box is located at the north or south pole, and the ODE (Eq. (22)) is greatly simplified. For θ_{0} = 0, the dimensionless wave propagation equation becomes (60)
At the southpole (i.e. θ_{0} = π), the term 1 − U′ in Eq. (60) is replaced by 1 + U′. We note that this equation is reminiscent of the differential equation for Rossby waves in the βplane, that is, and there is a constant df∕dy = β, with f the Coriolis parameter (e.g. Miles 1961; Grimshaw 1975b; Gliatto & Held 2020). However, we cannot make a direct comparison at corotation, because the singularity in the equations for Rossby waves and inertial waves is not of thesame order. We have a secondorder pole around the corotation, while only firstorder poles are found in the aforementioned studies. In fact, Eq. (60) is similar to the wave equation in stratified zsheared flows (e.g. Jones 1968).
In our polar configuration, the ycoordinate is now the axial distance, and this means that the mean flow has a cylindrical profile. Sucha rotation profile is expected in giant planets such as Jupiter and Saturn (Kaspi et al. 2017; Galanti et al. 2019, respectively) as a natural outcome of the ProudmanTaylor theorem for fastrotating bodies. The propagation and dissipation of inertial modes of oscillations in the presence of critical levels for this kind of mean flow have been investigated by Baruteau & Rieutord (2013) in a spherical shell.
3.4.1 Analytical solutions with constant shear
Analytical solutions of the ODE Eq. (60) are difficult to find for general profiles of the mean flow (e.g. a quadratic mean flow profile). A linear mean flow profile, on the other hand, has analytic solutions, which is why we use such a profile in this section (i.e. U = R_{o}y, with R_{o} the shear Rossby number that is taken as constant here). Equation (60) then becomes: (61)
where y_{0} = ω∕(k_{x}R_{o}) and α_{k} = k_{z}∕k_{x} is the vertical to longitudinal wavenumber ratio. When the box is located at the south pole, the lefthand term in the bracket is in the numerator. This equation takes the form of Whittaker’s equation (see Abramowitz & Stegun 1972), and solutions can be written in terms of the Whittaker functions M: (62)
where ỹ = 2k_{⊥}(y − y_{0}), (63)
and A and B are complex constants given by boundary conditions. The Whittaker function M_{0,−μ} allows a quite straightforward analytic continuation: (64)
By consequence, the solution below the critical point y = y_{0} is: (65)
Although the Whittaker functions do not feature precisely as wavelike forms, we can already get a good idea of the attenuation factor thanks to analytic continuation, as will be shown in the following section.
It is important to point out that μ can be real or complex depending on the value of (66)
which will simply be denoted by R in the following. This can drastically change the behaviour of a wave passing through the corotation. A necessary, but not sufficient, condition to find an instability is that R < 1∕4, as we will demonstrate in Sect. 3.4.4. This condition is similar to the MilesHoward theorem for stratified zsheared flow (Miles & Howard 1964; Lindzen 1988). In these studies, the prerequisite for instability is that Ri < 1∕4, where Ri is the Richardson number, namely the squared ratio of the BruntVäisälä frequency to the vertical (or radial in global spherical geometry, Alvan et al. 2013) shear. In our model, unlike cases where the box is tilted, a WKBJ analysis can typically be performed for a linear mean flow, provided that R ≫ 1∕4, in line with the condition of stability derived in the coming sections and detailed in Appendix C.
These various situations regarding the value of at the north and south poles are illustrated in Fig. 3. We stress the particular case where R_{o} = 1 (R_{o} = −1) at the north (south) pole, and where the differential equation and its solutions take a quite simple form: (67)
Solutions are then fully evanescent for such shears. One can notice that Eq. (67) is the same far from corotation for any mean flow.
Finally, it is clear from Fig. 3 that wave propagation is the same at the north and south poles provided a Rossby number of opposite sign. As a result, only the equations at the north pole will be treated in the following, and the word ‘pole’ now refers to the north pole.
3.4.2 Frobenius method at the pole
Though analytic solutions are known, it is still useful to determine Frobenius solutions near corotation for two main reasons. First, these solutions can be derived for any mean flow profile near the corotation. Close to corotation, the mean flow is approximated by a Taylor expansion at the firstorder . Secondly, Frobenius solutions may feature wavelike forms, which is helpful for physical interpretation. Therefore, Eq. (60) can be written near corotation as (68)
where and . The indicial equation gives: (69)
In the two next subsections, we examine the two cases where μ is imaginary or real.
Fig. 3 Diagram showing how compares to 1∕4 for two positions of the box (θ_{0} = 0 and π), as well as for a range of Rossby numbers R_{o} and a range of the vertical to longitudinal wavenumber ratio α_{k}. The solid and dashed black lines mark where at the north and south poles, respectively. The purple domain shows where , and the white region where , regardless of whether the box is at the north or the south pole. In the dark red region, and and vice versa in the orange region. 
3.4.3 Theoretical stable regime (R > 1/4)
We address here the case where R > 1∕4. The same analysis from Sect. 3.3 can be carried out to determine how a wave behaves upon crossing the corotation. The solutions of the indicial equation can be recast as (70)
The firstorder solutions to Eq. (68) in the vicinity of y_{0} are: (71)
for . One can recover the same form of Frobenius solutions as in Alvan et al. (2013), who examined radially stratified mean flows in spherical geometry. As , the wave action flux Eq. (42) reduces to (72)
that is, injecting the solutions on both sides of the critical level: (73)
This formulation is quite similar to the expression of the Reynolds stress (τ) in vertically stratified mean flows, which can be found in Booker & Bretherton (1967) in Cartesian geometry. We recall that, indeed, . Moreover, given Eq. (35), the Reynolds stress in our model reads . Using the polarisation relations for u, we recover the wave action flux in Eq. (72).
The prefactor i in the solutions (Eq. (71)) below the corotation does not affect the energy flow and simply indicates that the wave undergoes a phase shift of π∕2 through the critical level (see also Alvan et al. 2013). Above the critical level, the normalised Dopplershifted frequency satisfies sign (σ) = − ς, as was the case for the corotation in the inclined box. The sign is reversed below the critical level. Thus, the first solution of main amplitude a_{0}  carries its latitudinal flux of energy upwards (downwards) for ς = +1 (ς = −1), while the second solution transfers its energy in the opposite direction in the various cases. Therefore, the energy flux of an upward or downward wave is always attenuated by a factor (74)
This attenuation factor is shown versus R_{o} and α_{k} in Fig. 4. We observe that the wave is largely absorbed at the critical level and thus deposits most, if not all, its energy for most couples (α_{k}, R_{o}).
Fig. 4 Transmission rate T_{F} of the wave action flux across the corotation when the box is located at the pole for R > 1∕4. It is displayed against the Rossby number R_{o} and the ratioof the horizontal wavenumbers α_{k} = k_{z}∕k_{x}. The forbidden region where R < 1∕4 is shown in white, and the dark red cone corresponds to values of T_{F} that are lower than 10^{−9}. 
3.4.4 Possible unstable regime (R < 1/4)
We now deal with the case where R < 1∕4 (i.e. μ is real). Contrary to the situation where R > 1∕4, we can no longer assimilate solutions to wavelike functions. The exponential form of solutions for R < 1∕4 near the critical level reads (75)
and makes this region fully evanescent. Furthermore, the associated wave action flux is (76)
Without knowing the direction of the wave or the energy flux, since , it is difficult to assess the impact of the critical level on wave propagation, that is, whether it will attenuate waves or, on the contrary, amplify them.
Lindzen & Barker (1985) found a way to investigate the behaviour of internal gravity waves in the presence of a vertical shear, passing through a critical level in a regime similar to ours (Ri < 1∕4) where solutions are fully evanescent. Their work, which was carried out in local Cartesian geometry, was taken up by Alvan et al. (2013) in global spherical geometry applied to the radiative zone of solarlike stars and evolved stars. The goal of the method is to determine the reflection and transmission coefficients in a threezone model. The evanescent region where the Richardson number satisfies Ri < 1∕4 and where the critical level is located (zone II) is sandwiched between two propagating wave layers (zones I and III). Using a linear mean flow profile so as to establish solutions inside zone II, Lindzen & Barker (1985) and Alvan et al. (2013) both used continuity relations of the perturbed vertical (or radial) velocity and its derivatives at the interfaces between zones in order to get the transmission and reflection coefficients. The critical level is located in the middle of zone II which has a width δ. By consequence, the reflection and transmission coefficients depend, in their works, on the shear and, more precisely, the Richardson number, as well as on the width δ. They both found that, depending on the Richardson number and δ, the reflection and transmission coefficients can be greater than one, meaning that the wave can be overreflected and/or overtransmitted and thus extract energy and angular momentum fluxes from the mean flow, which can, in turn, lead to potential shear instabilities after successive encounters of the wave with the critical layer. However, this result is conditioned by the geometry of the model. As shown by Lindzen (1988) in his review, and references therein, models with one or even two layers with evanescent waves, and potentially a wavelike region, do not allow such phenomena. A first region that allows the wave propagation is mandatory and is combined with a ‘sink’ that pulls the wave to cross the critical level. According to Lindzen & Barker (1985) and Lindzen (1988), the nature of the sink for wave flux can be either another propagative region or an evanescent region that is, as in zone II, subject to friction processes. Given this peculiar geometry, instabilities can occur under boundary conditions that allow the wave to successively return to the critical level. Many studies have tried to relate overreflection and shear instability for a specific wave geometry (see in particular the reviews of Lindzen 1988; Harnik & Heifetz 2007, for internal gravity waves and Rossby waves).
In the present study, we do not further investigate shear instability by doing, for instance, a temporal analysis to estimate the waves’ growth rate (as in Lindzen & Barker 1985; Watts et al. 2004, who considered an initial value problem). On the contrary, we give arguments of necessary but not sufficient conditions to find instabilities, such as R < 1∕4. It is important to note that R is constant in the whole domain for a linear mean flow profile, and thus one is stuck with either a propagative (stable) or an evanescent regime. Therefore, finding an adequate geometry to allow overreflection and overtransmission requires at minimum that the Rossby number is not the same in the whole domain by using, for instance, a nonlinear mean flow profile. Furthermore, in the particular case where R = 0 (i.e. R_{o} = 1 or α_{k} = 0 in Eq. (60) when the wave with k_{z} = 0 propagates in the (x, y)plane), a necessary condition for instability is given by the inflection point theorem (Schmid et al. 2002). This theorem is particularly used to study barotropic instabilities for Rossby waves (see e.g. Lindzen & Tung 1978). In other words, a necessary condition to have unstable modes for R_{o} = 1 is that U″ is cancelled out in the domain of wave propagation.
We summarise in Table 2 the main analytical results of Sects. 3.3 and 3.4 regarding wave and wave action flux transmission, when the box is tilted relative to the rotation axis with a random angle, at the north pole, and at the equator, in the inviscid limit. We note that when ‘no’ is given in both the attenuation and amplification columns, the wave is fully transmitted across the critical level, regardless of the wavenumbers and of the mean flow profile.
Summary of analytical results at each critical level when the box is inclined or not with respect to the rotation axis, in the [0, π∕2] quadrant, and for positive wavenumbers.
4 A threezone numerical model
In order to test the analytical predictions of the previous section, we built up a threezone numerical model to simulate waves passing through critical levels. A similar model has been used, for instance, by Jones (1967) to explore the behaviour of internal gravity waves passing through critical levels in a fluid with rotation and vertical shear. In our model, we solve the two firstorder ODEs satisfied by v and Π, the combinationof which led to the wave propagation equation, Eq. (22). By imposing boundary conditions such that waves satisfy the dispersion relations (see also Appendix D.1), we examine the dynamics of inertial waves propagating in the shear region. Also, whenever possible, we analytically calculate the wave transmission and reflection coefficients as the wavelike solution crosses the shear region.
Fig. 5 Mean flow profiles used in the threezone numerical model against y. In the noshear regions I (y < 0) and III (y > 1), the mean flow is uniform and set respectively to U = 0 and U = Λ. The shearedregion II (greyshaded) can have a linear, square, or cubic mean flow profile. 
4.1 Description of the model
The mean flow profile that is used in the threezone model is illustrated in Fig. 5. The zone with shear (zone II) is surrounded by two noshear regions, one with no mean flow (zone I) and one with a uniform mean flow (zone III). In the whole domain, the mean flow profile that we adopt is expressed as (77)
where U(y) is continuous at each interface and n is an integer: n = 1 for a linear shear flow, n = 2 for a square shearflow, or n = 3 for a cubic shear flow (see also Fig. 5). In zone I, we assume that there is an incident wave that enters the shear zone as well as a wave that is reflected at the interface between zones I and II or in zone II, that is: (78)
where A_{I} and A_{R} are the amplitudes and k_{I} and k_{R} are the wavenumbers of the incident and reflected waves, respectively. We further imposed as a boundary condition in zone III a transmitted wave that propagates towards positive yvalues: (79)
where A_{T} and k_{T} are the amplitude and wavenumber of the transmitted wave, respectively. We imposed A_{T} = 1 without loss of generality and computed the remaining amplitudes A_{I} and A_{R}. More detailson the solutions and the dispersion relations of the waves in zones I and III can be found in Appendix D.1. We ensured that the transmitted wave carries energy upwards by deriving the wave action flux in zones I and III (see Appendix D.2).
We imposed the continuity of the latitudinal velocity and reduced pressure at the interfaces (at y = 0 and y = 1). By doing so, the wave action flux is continuous at both interfaces. Thus, in the absence of critical points, the wave action flux is conserved in the whole domain, namely , where is the wave action flux of the transmitted wave, and , where and are the incident and reflected wave action fluxes, respectively.
To solve the ODE in the three zones, and in particular near singularities, we used MATLAB’s solver ode15s, which is suitable for solving stiff differential equations (Shampine & Reichelt 1997). To avoid strict singularities at , we added a small friction, σ_{f} = 10^{−8}, to our set of units. Given the boundary conditions, the numerical solver deals with two firstorder ODEs for v and Π, which take the form (80)
and where we recall that s = σ + iσ_{f} is the modified Dopplershifted frequency due to Rayleigh friction. While A_{T} = 1 is imposed by the boundary condition, we computed A_{I} and A_{R} by comparing numerical solutions of the system Eq. (80) at y = 0 with the definition of velocity in zone I (Eq. (78)) and its associated reduced pressure (see Eq. (D.9)).
4.2 Numerical exploration at the pole for a constant shear
4.2.1 Reflection and transmission coefficients
In most cases, for any inclination of the box and any mean flow profile, there is no analytical solutions in zone II. Nevertheless, we have shown in Sect. 3.4.1 that, when the box is at the pole and for a linear mean shear flow, solutions can be found in terms of Whittaker functions. In this section, we will find reflection and transmission coefficients similarly to how it was done in Lindzen & Barker (1985) and Alvan et al. (2013), though there are a few differences. In particular, our present study differs from the latter due to the treatment of inertial waves in convective regions (instead of gravity waves in stably stratified radiative regions in their case) with a latitudinal shear (instead of a vertical or radial shear). Our study, however, uses a local Cartesian model as in Lindzen & Barker (1985). Moreover, our boundary conditions are different, as detailed in Sect. 3.4.4, and the thickness of our shear region is fixed to one in scaled units, while Lindzen & Barker (1985) and Alvan et al. (2013) leave the thickness δ as a control parameter. We also checked for the existence of a critical level in the shear zone and the frequency range that delineates the regimes with and without the critical level.
We considered the perturbed reduced pressure Π and velocity v to be continuous at the interfaces y = 0 and y = 1. In the presence of the critical level y_{0} in zone II, we have a set of four analytical solutions whose values at y = 0 and y = 1 allow us to determine the reflection and transmission coefficients. The solutions to the wave propagation equation in zones I, II (below and after the critical level), and III are: (82)
where A and B are complex coefficients that we will express below. We recall that A_{T} = 1 in our numerical model. In the shear region (zone II), the reduced pressure perturbation is given by Eqs. (40), which, at the pole, can be recast as (83)
with R_{o} = Λ. In the regions with no shear (zones I and III), Π takes a simpler expression, with R_{o} = 0 and σ = ω in zone I, and σ = ω − k_{x}U(1) = ω − k_{x}R_{o} in zone III (noted as σ_{3} in the following). We note that while the reduced pressure is kept continuous to conserve the wave action flux across the interfaces, the first derivative of the latitudinal velocity v′ is not necessarily continuous at the interfaces. To find the transmission and reflection coefficients, we solved the system of equations that consists of matching conditions at interfaces as follows:
 1.
v is continuous at y = 0 : A_{I} + A_{R} = v_{W}(0);
 2.
Π is continuous at y = 0 :
 3.
v is continuous at y = 1 : e^{ikT} = v_{W}(1),
 4.
Π is continuous at y = 1 :
with the Whittaker functions v_{W}. At the interfaces below and above the critical level y_{0} (i.e. y = 0 and y = 1), we have (84)
It should be noted that the first derivative of the Whittaker functions can be computed either numerically or analytically via the relationships in Abramowitz & Stegun (1972).
The equations in the above continuity relationships 1 to 4 are independent two by two (1 & 2 and 3 & 4), and A and B can be found first: (85)
The amplitude of the incident and reflected waves can be written, in terms of A and B, as: (86)
The transmission and reflection coefficients are then: (88)
We emphasise that these factors depend notably on the location of the critical level and on the inertial frequency, which was not the case in Lindzen & Barker (1985) and Alvan et al. (2013).
We display in Figs. 6 and 7 the transmission and reflection coefficients as a function of the Rossby number R_{o} and the normalised inertial frequency ω in the regimes R < 1∕4 and R > 1∕4, respectively. We chose k_{x} = k_{z} = 0.1 in these plots, and R = 1∕4 gives , which delineates the two regimes, as we can see in Fig. 3. Areas that are hatched do not possess the corotation point σ = 0. In addition, the wavenumberk_{T} was chosen with a positive sign in these regions to maintain an upward propagating wave (see Appendices D.1 and D.2). Areas that are not hatched feature a critical level, according to the Table D.1. In the case where R < 1∕4, overreflection and overtransmission are both possible (see Fig. 6). One should notice that for R_{o} > 1, we always have by definition of R (89)
regardless of y, which makes the solutions of Eq. (61) tend towards pure exponential functions (i.e. without any imaginary part). Also, we do not see any overreflection, nor overtransmission, in the hatched areas where there is no corotation point. This highlights the essential role of the critical level in inducing the overreflection or overtransmission of inertial waves crossing the shear region in this regime.
The analysis of the regime where R > 1∕4 (in Fig. 7) is trickier. According to our discussion in Sect. 3.4.3, we expect a strong attenuation of the wave and of the wave action flux, as shown in Fig. 4. From this figure and for α_{k} = 0.1, the damping is very strong for low positive Rossby numbers. This tendency is also found for both transmission and reflection coefficients. Nevertheless, one can also observe an unexpected regime of overtransmission near R_{o} = 0.8 and low frequency ω. Still, we must not forget that solutions in this regime, even near the critical level (see Eq. (71)), are not rigorously equivalent to wavelike functions. In particular, the amplification term that can be found at the first order in the Frobenius solutions becomes more prominent as the thickness of the shear zone is larger. This is especially true for the transmission coefficient. Assuming that Eq. (71) holds throughout zone II and corresponds to upward and downward waves, the transmission coefficient is modulated by 1 − y_{0}^{1∕2}∕0 − y_{0}^{1∕2}, the amplitude ratio between the transmitted and incident waves. This term can be greater than one in the shear region. In particular, it is always greater than one when y_{0} < 0 (i.e. no critical level in the regime R > 1∕4; hatched areas in Fig. 7). On the contrary, this ratio is not present for the reflection coefficient since R is a function of the incident and reflected waves evaluated at y = 0. Though this handwaving explanation does not formally demonstrate the origin of this amplification, it stresses the important role of the shearregion thickness and, more generally, of the geometry of the model.
In order to clarify whether the amplification is due to the geometry or the critical level, we need to investigate how the wave action flux changes before and after the critical level. The wave action flux is indeed the relevant quantity to investigate energy flux exchanges at a critical level.
Fig. 6 Transmission coefficient (T, left panel) and reflection coefficient (R, right panel) when the box is at the pole. The coefficients are plotted in the regime R < 1∕4 (a possibly unstable case, see Sect. 3.4.4) as a function of the Rossby number R_{o} and the inertial frequency ω. The hatched areas do not feature critical points and correspond to regions where ω > k_{x} R_{o} in our peculiar geometry (see Appendix D.2 for this particular matter). Vertical and longitudinal wave numbers are fixed: k_{x} = 0.1 and k_{z} = 0.1. Moreover, the contours that correspond to coefficients R and T that are equal to one are indicated by solid black lines. Crosses mark the set of parameters used in Fig. 9for the analysis of the behaviour of the velocity in the threelayer model. 
Fig. 7 Same as Fig. 6, but for R > 1∕4. When R_{o} < 0 in these panels, hatched areas mark instances where ω > 0 (see Appendix D.2). 
Fig. 8 Ratio of the wave action flux above and below the critical level against the Rossby number and the inertial frequency for R >1∕4 (left panel) and R < 1∕4 (right panel) when the box is at the pole. As in Fig. 7, k_{x} = k_{z} = 0.1 and hatchedzones represent areas without critical levels. Again, crosses indicate the set of parameters chosen to analyse the behaviour of the velocity in the threelayer model. 
4.2.2 Wave action fluxes below and above the shear region
Since v and Π are continuous at the interfaces between the shear and noshear regions, the wave action flux is preserved and continuous in all three zones in the absence of friction and critical levels. However, it is discontinuous at the corotation point, as demonstrated in Sects. 3.4.3 and 3.4.4. Given the amplitude of the incident and reflected waves (Eq. (86)), we can calculate the ratio of the wave action flux below and after the corotation (see Appendix D.2 for the detailed calculation): (90)
The signs + or − can be chosen in regards to the wave action flux of the transmitted wave that can be positive or negative, depending on the presence of the critical level, while the energy flux is always positive in order to have an upward propagating wave in zone III (see Appendix D.2 for a more detailed discussion). This wave action flux ratio is displayed in Fig. 8 in the two regimes, R ≶ 1∕4. As expected, this ratio is equal to one when no critical level is present (hatched areas). Unlike in the previous section, the regime where R > 1∕4 no longer has amplification areas, and everywhere. This supports the idea that the critical level has nothing to do with the amplification phenomenon observed in the left panel of Fig. 7. As already observed in Fig. 4, the damping due to the critical level is strong except for Rossby numbers close to the threshold between the two regimes. Moreover, means that since the minus sign is taken in Eq. (90). Therefore, no overreflection due to the critical level is expected in this regime. The other regime (R < 1∕4, right panel) features areas where the wave is overreflected for (i.e. when ) and areas where the wave is overtransmitted for . For the first inequality (), the threshold between under and overreflection (around R_{o} ≈ 0.9) is the same as for the reflection factor (in the right panel of Fig. 6). For the second one (), the comparison with the transmission factor (in the left panel of Fig. 6) is more questionable. Still, these two points suggest that the critical level can induce the overreflection and overtransmission phenomena in the regime where R < 1∕4.
4.2.3 Numerical solutions
In the previous sections, we examined how the shear parameter and the inertial wave frequency impact the reflection and transmission coefficients as well as the wave action flux. We now study particular cases of wave propagation through the critical level for fixed sets of parameters in both regimes, R ≶ 1∕4. To do that, we numerically solved Eqs. (80) for the threelayer model described in Sect. 4.1, for θ_{0} = 0 and a linear shear flow (n = 1 in zone II). We selected three pairs of values for the inertial frequency and the shear, two in the regime R > 1∕4 and one in the regime R < 1∕4. These values are marked by crosses in Figs. 6–8. In each case, the latitudinal velocity and the wave action flux were successively calculated through the three zones and plotted in Fig. 9. The numerical solution, which was computed by imposing the boundary condition A_{T} = 1 and the continuous interfacial conditions for v and Π at y = {0, 1}, is the sum of the incident and reflected waves in zone I and is equal to a transmitted wave in zone III. The expressions for the incident, reflected, and transmitted waves are given by Eqs. (78) and (79) (see also Appendix D.1). In the shear region (zone II) of the upper panel of Fig. 9, the Whittaker solution has been added and it matches perfectly with the numerical solution below and above the critical level in each case. Moreover, Frobenius approximations for the latitudinal velocity (Eq. (71) when R > 1∕4 and Eq. (75) when R < 1∕4) and for the wave action flux (Eq. (73) when R > 1∕4 and Eq. (76) otherwise) have also been included. The coefficients a_{0} and b_{0} were determined by matching the numerical solution and its derivative to the Frobenius approximation for the velocity close to the critical level. For both the latitudinal velocity and the wave action flux, this firstorder approximation gives satisfactory agreement with numerical solutions, although a slight deviation (regarding velocity) from the numerical solution can be observed as one moves away from the critical level. In particular, it should be mentioned that for the far left panels, corresponding to an upward wave is sufficient to correctly fit the numerical solution, meaning that the counterpropagating wave in zone I is reflected at y = 0. However, for the middle and right panels, it is not clear whether the firstorder Frobenius solutions can be reconnected to the incident and reflected waves at y = 0.
We now examine attenuation or amplification phenomena in each column of panels in Fig. 9. In the left panels, for which R > 1∕4, the latitudinal velocity is strongly attenuated at the critical levels and so is the wave flux action, as we can expect from the left panel of Fig. 8. While the transmitted wave is totally absorbed, the reflected wave remains, which is consistent with analytical values of the transmission and reflection coefficients in Fig. 7 (see white crosses). In the middle panels, where we also have R > 1∕4, the wave is overtransmitted but not overreflected, which is also consistent with Fig. 7 (see black crosses). In view of the wave action flux, the amplification of the transmitted wave does not seem to be related to the critical level because this quantity is greatly reduced after the critical level (see also the white cross in the left panel of Fig. 8). The third column of panels now depicts the regime where R < 1∕4. The wave is overreflected by a factor of ~1.5 and overtransmitted by a factor of ~2, in concordance with the reflection and transmission coefficients plotted in Fig. 6. The wave action flux is negative, and by a factor of three, as observed in Fig. 8. These three case studies reinforce the idea from Booker & Bretherton (1967) in the case of stratified zsheared flows: that the wave energy can be lost to the mean flow or, on the contrary, that the wave can take energy from the mean flow.
Fig. 9 Numerical outputs of the threelayer model. Top: real part of the latitudinal velocity v against y. The quantities v_{num}, v_{I}, v_{R}, v_{T}, v_{F}, and v_{W} are the numerical, incident, reflected, transmitted, firstorder Frobenius, and Whittaker velocities, respectively. Bottom: wave action flux against y. The quantities , , , and are the numerical, transmitted, incident and reflected, and Frobenius wave action fluxes, respectively. For all panels, θ_{0} = 0, the mean flow is linear in the greyshaded shear regions, and the critical level is marked by dashed lines. The horizontal wave numbers are set to k_{x} = k_{z} = 0.1. From left to right: (i) ω = 0.02 and R_{o} = 0.3 (R > 1∕4), (ii) ω = 0.002 and R_{o} = 0.8 (R > 1∕4), and (iii) ω = 0.09 and R_{o} = 1.8 (R < 1∕4). 
4.3 Numerical exploration at constant shear when the box is inclined
We now investigate wave propagation through the different critical levels when the box is inclined with respect to the rotation axis. We still assume that the shear region (zone II) has a linear shear flow profile (n = 1). In contrast to the polar configuration, we do not have analytical solutions to the ordinary differential equation. Instead of going for an extensive numerical investigation of the parameter space, we instead focus on the dynamics of inertial waves in our threelayer model as they cross the critical levels and σ = 0. Our results are presented in Fig. 10 for a box inclination of 10° relative to the rotation axis, a shearfixed to R_{o} = 0.3, and wavenumbers set to k_{x} = k_{z} = 1. The value of the frequency ω determines the existence and the nature of the critical level as detailed in the Table D.1. As in Fig. 9, we plot in each column of Fig. 10 the latitudinal velocity and the wave action flux. From left to right, we illustrate our results for the critical levels , σ = 0, and (i.e. there are two critical levels in the rightmost panels). One can notice that the firstorder Frobenius approximation is no longer in good agreement with the numerical solution in the entire shear region, though it remains a reasonable approximation in the vicinity of a critical level. Unlike the polar case, the discrepancy far outside the critical levels is due to the linear approximation that the governing ODE takes around critical levels (see Eqs. (44) and (56)).
In the left panels of Fig. 10, the reflected and transmitted waves are strongly attenuated at the critical level y = y_{+}. Part of the wave energy is laid down to the mean flow, as corroborated by the drop in the wave action flux. In the middle panels, we do not see any discontinuity at the corotation y = y_{0}, which is in line with the theoretical analysis in Sect. 3.3.2 for a constant shear. However, the wave is overreflected and overtransmitted, possibly due to the polynomial form of the solutions in the Frobenius series around y = y_{0}. In the right panels, the wave encounters successively critical levels at y = y_{0} and at y = y_{−}. Although the wave going through the shear region is not attenuated at the corotation y = y_{0}, it is completely absorbed at the second critical level y = y_{−} where the wave action flux drops to zero. This is consistent with the transmission coefficient in the left panel of Fig. 2. The latitudinal velocities displayed in the topleft and topright panels of Fig. 10 support the concept of a valve effect. Indeed, according to our analysis with the Frobenius method in Sect. 3.3.1, and given the shear and wavenumbers of Fig. 10, the attenuation is strong for a downward wave meeting the critical level y_{+} (first panel), while the attenuation is strong for an upward wave that meets the critical level y_{−} (third panel). Before the critical level y_{−} and after the critical level y_{+}, we observe fast oscillations of shorter and shorter periods close to the critical level, as already evidenced by Booker & Bretherton (1967). The analysis to determine how the wave is reflected in the shear zone can hardly be taken any further because Frobenius solutions are not fully separable into upward and downward waves.
We emphasise that the behaviour of the wave at corotation y = y_{0} when the box is at the pole stands out as being clearly different from the case when the box is inclined for a linear mean flow profile. This is particularly true in terms of the absolute value of the wave action flux that is subject to rise and drop in the polar configuration, whereas it remains conserved when the box is inclined. In this inclined case, the only way for a wave to be attenuated without friction is for it to meet critical levels . Depending on the critical level encountered, an upward or downward wave will not be attenuated in the same way, which is a phenomenon also known as the valve effect. In addition, we no longer observe amplification due to the critical level, but there is still ‘geometric’ amplification (for instance, in the middle panels where we can observe overtransmission and overreflection), which can be explained by the exponential form of the Frobenius series.
Fig. 10 Same quantities as in Fig. 9 (when the box is at the pole), but for a box tilted by 10 degrees relative to the pole and for different values of the shear, wavenumbers, and inertial frequencies. We note that, unlike Fig. 9, we do not have analytical solutions in the shear region. In all panels, the horizontal wavenumbers are set to k_{x} = k_{z} = 1 and the shear is R_{o} = 0.3. From left to right: inertial frequencies set to ω = 0.31, ω = 0.16, and ω = 0.1. Third panels: we indicate y_{0} and y_{−} with dasheddotted and vertical dashed lines, respectively. 
4.4 Numerical results with a nonconstant shear
The choice of a linear mean flow profile allows the resolution of the ordinary differential equation at the pole and a simpler implementation of the Frobenius method. Nevertheless, the correspondence between a global and a local mean flow as presented in Sect. 2.2 involves higherorder terms than a simple linear dependence. Therefore, it is important to examine the effect of different mean flow profiles with nonzero U″. Now assuming n > 1 for the shear flow profile in zone II (see Eq. (77)), the Frobenius method in the inclined and polar cases still holds provided that (91)
Moreover, we specify that the critical level in y < 0 for even values of n is not examined given our numerical setup where the shear region is located in the range y ∈ [0, 1]. The conditions to have critical levels inside the shear region (zone II) are the same as those for a constant shear (see Appendix D.2).
We show in Table 3 numerical values of the reflection and transmission coefficients along with the ratio of the wave action flux below and after the critical level in the six parameter sets illustrated in Figs. 9 and 10 for linear, square, and cubic mean flow profiles. For all cases, R, T, and change quite significantly between n = 1, 2, and 3. At the pole, the (‘geometric’) overtransmission found for R_{o} = 0.8 disappears for n = 2 and n = 3, where the reflected and transmitted waves are strongly attenuated. Similarly, in the case where R_{o} = 1.8, the overtransmission and overreflection disappear when n = 3, whereas , which entails that the transmitted wave has taken wave action flux from the mean flow just after the critical level. We find that the Frobenius method does not give consistent results for the case where R_{o} = 1.8 and n = {2, 3}. When the box is tilted and for the critical level y_{+}, the reflected wave is more attenuated than the transmitted one, which is consistent with the stronger attenuation of the counterpropagating wave discussed in the previous section. At the corotation and when the box is inclined (fifth row of Table 3), the wave action flux remains the same in the whole domain for linear or nonlinear mean flow profiles, in agreement with the wave action flux derived analytically in Eq. (59). We speculate that the overtransmission and overreflection are due to the polynomial form of the solutions in the shear region. These wave amplifications may be related to shear instabilities in this particular threelayer configuration, and they are probably not linked with the presence of a critical level whose implications on the flow are well diagnosed by the wave action flux. In the last case (sixth row), where two critical levels coexist, the transmitted wave is less attenuated at y = y_{−} from n = 1 to n = 3, but it remains more attenuated than the reflected wave, as discussed in the previous section. Again, no jump is found in the wave action flux at the corotation despite analytical predictions.
The same analysis was carried out for different amplitudes of the Rayleigh friction force, up to σ_{f} = 10^{−2}. Of course, the wave action flux is no longer constant, but we observe that all three parameters of Table 3 change very little compared to the case where σ_{f} = 10^{−8}. This result is consistent with Alvan et al. (2013) in the context of gravity waves and vertical shear. We comment that while a low friction is mandatory in the numerical code to solve the ODE at singularities y_{±}, it is not the case at the corotation y = y_{0}.
Reflection R and transmission T coefficients and the ratio of the wave action flux above and below a critical level for a linear, square, and cubic mean flow (n = 1, 2, and 3, respectively).
5 Astrophysical discussion
5.1 Latitudinal differential rotation in stars
In stars, latitudinal differential rotation is often characterised by the difference in rotation frequency between the equator and the pole, that is the quantity ΔΩ = Ω_{eq} − Ω_{0}, where Ω_{eq} is the rotation frequency at the equator (e.g. Barnes et al. 2017). We will now refer to ΔΩ as the shear contrast. Different regimes are distinguished according to the value of ΔΩ: antisolarlike rotation for ΔΩ < 0, cylindrical rotation for ΔΩ≪ 1, and solarlike rotation rotation for (nottoolow) positive ΔΩ (e.g. ΔΩ∕Ω_{0} ≃ 0.3 for the Sun). Several works based on threedimensional numerical simulations have explored the range of physical parameters leading to each aforementioned regime in stars and in giant planets (e.g. Gastine et al. 2013; Varela et al. 2016; Beaudoin et al. 2018) In particular, Brun et al. (2017) derived a criterion based on mixing length theory and calibrated with threedimensional simulations that determines the rotation profile of a solarlike star. This criterion is based on the fluid Rossby number R_{of}, defined as (92)
where R_{of,⊙} = 0.89 is the solar fluid Rossby number and Ω_{*} and M_{*} are the mean rotation and the mass of the star respectively, normalised with their solar values. Brun et al. (2017) highlighted the following three regimes: (i) R_{of} > 1 for antisolarlike rotation, (ii) 0.3 < R_{of} < 0.9 for solarlike rotation, and (iii) R_{of} ≲ 0.3 for cylindrical rotation. Furthermore, they introduced the shear contrast ΔΩ_{S} at the colatitude 30° as (93)
since the rotation frequency is often illdefined at low colatitudes in threedimensional numerical simulations in spherical geometry. From their threedimensional simulations, Brun et al. (2017) obtained the following scaling: (94)
where ΔΩ_{S,⊙}≃ 565 × 10^{−9} s^{−1} is the solar value of ΔΩ_{S} calculated from García et al. (2007). Using the mean flow profile (Eq. (4)), ΔΩ_{S} can be related to our shear contrast ΔΩ via (95)
We show in Fig. 11 the quantity ΔΩ∕Ω_{0}, expressed inEq. (95), versus the age of solarlike stars for K to G spectral types. To compute this quantity, we used grids of the onedimensional stellar evolution code STAREVOL (see Amard et al. 2019, for details of the code). In light of Fig. 11, the stars in the premain sequence (age ≲ 100 Myr) exhibit cylindrical rotation as they are fast rotating. During the main sequence, stars mostly feature solarlike rotation, while antisolarlike rotation is observed at the end of the main sequence from 0.8 to 1.1 M_{⊙}. According to Fig. 11, a limit on the absolute value of the normalised shear contrast can be set to ΔΩ∕Ω_{0} < 0.5. However, as already stressed by Benomar et al. (2018), the latitudinal shear inferred by asteroseismology can be much larger than predicted by numerical simulations. This can actually be inferred by comparing the shear factors in their work (see Table S3 from the supplementary materials in Benomar et al. 2018) with ours given in Fig. 11. Moreover, according to their study, cylindrical and antisolar differential rotation are hardly unambiguously detectable. Finally, one should recall that the Brun et al. (2017) scaling laws given in Eqs. (92) and (94) are derived for K and G spectral type stars only.
Since we now know values of the shear contrast, we can calculate the ‘shear’ Rossby number, R_{o} = U′∕(2Ω), given by thefollowing relationship: (96)
which has been derived from Eq. (5) by keeping only zeroorder terms in y. Taking χ ≃ 0.3 as a representative value of the shear contrast for main sequence G and K stars, we find that the Rossby number is maximal when θ_{0} ≃ 55°, its maximum value being R_{o} ≃−0.17. In particular, R_{o} ≃−0.013 for θ_{0} = 10° and R_{o} ≃−0.076 for θ_{0} = 80°. These values are useful for interpreting wave flux action transmission at critical levels , considering Fig. 2 (we recall that at the corotation the wave action flux is fully transmitted). A downward (upward) propagating wave through () is: (1) totally absorbed provided that k_{z} ≳ 0.1k_{x} (α_{k} ≳ 0.1) at θ_{0} = 10°, and that k_{z} ≳ k_{x} (α_{k} ≳ 1) at θ_{0} = 80°, (2) strongly attenuated for k_{z} ~ 0.1k_{x} (α_{k} ~ 0.1) at θ_{0} = 80°, and (3) fully transmitted given that k_{z} ≪ 10^{−1}k_{x} (α_{k} ≪ 10^{−1}) for both inclinations.
These results also hold for antisolarlike differential rotation since the transmission factor is a function of R_{o}. For larger values of R_{o} (i.e. for larger values of the shear contrast), waves are less damped at critical levels at a given α_{k} = k_{z}∕k_{x}. The connection between this ratio of the vertical and azimuthal wavenumbers in the local model and an equivalent ratio of global wavenumbers in the spherical geometry is not straightforward. A first hint can be to state that k_{z} ~ k_{r}, where k_{r} is the wavenumber in the global radialdirection, while k_{x} ~ m∕(r_{0} sinθ_{0}), where m is the azimuthal order of the considered mode of the tidal potential (when m≠ 0; Zahn et al. 1997). Then, we get α_{k} ≡ k_{r}r_{0} sinθ_{0}∕m ≡ 2π r_{0}∕λ_{r} × sinθ_{0}∕m by introducing λ_{r}, the radial wavelength of the tidal wave. In the case where r_{0} > λ_{r}m, we should thus be in the regime where the tidal wave is attenuated.
Fig. 11 Absolute values of the shear contrast normalised by the rotation at the pole against the age of K and G spectral type stars. Solid lines feature 0.3 < R_{of} < 0.9, dashed lines R_{of} ≲ 0.3, and transparent lines R_{of} > 1. 
5.2 Cylindrical differential rotation in Jupiter and Saturn
Mathis et al. (in prep.) have developed an equatorial model to examine inertial wave properties in the outer convective layers of giant gaseous planets, such as Saturn and Jupiter, which are subject to cylindrical differential rotation. In their model (built in cylindrical coordinates), they derived a Schrödingerlike differential equation for under the anelastic approximation, where ρ is the density, r the axial distance coordinate, and v_{r} the axial velocity. For free inertial waves, their secondorder differential equation is: (97)
where l_{Θ} and m denote the equatorial and azimuthal wavenumbers, respectively, is the (linear) Dopplershifted frequency, and κ_{r} is the ‘axial’ epicyclic frequency defined as . For cylindrical differential equation, the corotation resonance results in critical cylinders (see Baruteau & Rieutord 2013) characterised by a critical axial distance r = r_{c}. The Taylor expansion of Eq. (97) at the first order around r_{c} gives: (98)
by setting the local Rossby number in cylindrical coordinates. By writing , Eq. (98) becomes very similar to our ODE Eq. (22) when the box is located at the south pole and for a constant shear (i.e. with , f = −1, and ). We note that, when the local shear box model is located at the equator, the latitudinal coordinate y is directed along the (vertical) rotation axis, whereas when the box is at the poles, y is the axial distance. That is why a polar configuration of the box best reproduces the ‘equatorial’ model of Mathis et al. (in prep.). Moreover, the convention of a plus sign in the Dopplershifted frequency explains why Eq. (98) is analogous to our wave propagation equation when the box is at the south pole rather than at the north pole.
In Fig. 12, we show Jupiter’s and Saturn’s local Rossby number from Mathis et al. (in prep.). Cylindrical differential rotation extends in the outer layer of the convective envelope of both planets, in agreement with the Juno and Cassini Grand Finale observations(Kaspi et al. 2017; Galanti et al. 2019, respectively). According to Fig. 3, where we notably plotted at the south pole (dark red and purple areas), we can assess the role of the critical level for wave transmission across the corotation, in terms of the wavenumber ratio α_{k}. In Fig. 12, the Rossby number satisfies R_{o} ≲ 0.27 for Saturn and R_{o} ≲ 0.07 for Jupiter. Given this range of values, two regimes can be evidenced for waves and wave action fluxes through the corotation: (1) waves are strongly attenuated for l_{Θ} ≳ m (see also Fig. 4 for R_{o} = −R_{o} as the transmission factor in this figure is plotted for the north pole) and (2) waves can be overreflected and overtransmitted for l_{Θ} ≪ m and can potentially lead to instabilities given specific boundary conditions.
To give an idea of the values these wavenumbers can take, we have listed in Table 4 three typical orbital states, where, in order, the asynchronous, eccentricity, and obliquity tides are supposed to be dominant (Ogilvie 2014). These states are described by the ‘spherical’ quadrupolar components of the dominant terms in the tidal potential, namely the degree l and the order m of the spherical harmonics. The analogy with the equatorial model is then made to get m, and l_{Θ} is chosen toapproximate, as best as possible, the behaviour of the Legendre polynomial around the equator with a simple trigonometric function ℜ{exp[i(l_{Θ}Θ + ϕ)]}, where ϕ is the appropriate phase (Mathis et al., in prep.). To find the associated wave attenuation at corotation for the three main tides, one can use Figs. 3 and 4 for α_{k} ≥ 1 and R_{o} ≲ 0.27 and look at the south pole (as Fig. 4 is plotted at the north pole, one has to take the opposite Rossby number). From Fig. 3, we can assess that waves excited by these tides are always in the socalled stable regime for these ranges of R_{o} and α_{k}, which excludes an amplification of these waves. Moreover, in Fig. 4, we also observe that waves are completely absorbed at corotation for our given ranges of parameters. By consequence, waves excited by the asynchronous, inclined, or eccentric tides in Jupiter and Saturn are expected to transfer all their wave action flux to the mean flow at corotation.
Fig. 12 Rossby numbers for Jupiter and Saturn as a function of the axial distance r. 
6 Conclusion and perspectives
The present study was motivated by the works of Baruteau & Rieutord (2013) and Guenel et al. (2016a,b), who showed that differential rotation can strongly affect the propagation and dissipation properties of (tidal) linear inertial waves. They considered different rotation profiles typical of stellar and planetary interiors and pointed out that tidal waves can deeply interact with zonal flows at corotation resonances, leading to intense wave energy dissipation, along with possible instabilities. In this paper, we have investigated the transmission of free inertial waves with latitudinal stratification and differential rotation at the corotation resonance (characterised by a zero Dopplershifted wave frequency) and, more broadly, at critical levels (any singularities of the governing secondorder wave propagation equation in the inviscid limit). For this purpose, we built a new local Cartesian box model with horizontal shear, modelling a small patch of the convective zone of a lowmass star or a giant planet. By considering the inclination of the local reference frame relative to the rotation axis, we have examined the effect on wave propagation through a critical level of a conical rotation profile at a general colatitude when the box is tilted, or of a cylindrical rotation profile when the box is at the north or south poles. These rotation profiles are inspired by those observed or expected in the Sun, lowmass stars, and the giant gaseous planets in our Solar System. Three critical levels can be identified when the box is inclined relative to the rotation axis: the corotation resonance and two other critical levels that arise from the inclination between the gravity and the rotation vectors, which are defined by a Dopplershifted frequency equal to plus or minus the latitudinal component of the rotation frequency. When the box is at the poles, critical levels are restricted to the corotation.
In order to diagnose the behaviour of a wave passing through a critical level for both aforementioned rotation profiles, we made use of an invariant called the wave action flux, which is independent of the latitudinal coordinate in a nondissipative fluid flow. This invariant was used when the ‘directional’ flux of angular momentum (here latitudinal) cannot be constructed easily from the mean perturbed velocity, as is the case, for example, in Lindzen & Tung (1978) for Rossby waves in planeparallel shear flows. The wave action flux has already been used in vertically stratified shear flows in the presence of rotation or magnetic fields to interpret the role of critical levels (Grimshaw 1975a, 1979; Andrews & McIntyre 1978; Mathis 2009; Mathis & de Brye 2012). Using the condition that this invariant is discontinuous at critical levels, we demonstrated in Sect. 3 that waves can be either fully transmitted, damped, or even amplified after passing through critical levels as a result of wave action flux exchanges. These different regimes of wave transmission are found with both conical and cylindrical rotation profiles; they depend on the critical level encountered, on the wave properties (e.g. the propagation direction and wavenumbers), and on the profile of the mean flow. Table 2 summarises the main analytical results.
We then compared our analytical results with a threelayer numerical model that comprises a shear zone, where the critical level is located, and two surrounding shearfree zones that allow incident, reflected, and transmitted waves. A difference with the analytical model is the introduction of a small dissipative force under the form of a Rayleigh friction (also called frictional force by Ogilvie 2009) to avoid strict singularities. This does not seem to affect the results since analytical and numerical results match quite well when using a powerlaw mean flow profile and varying the friction. This conclusion is also shared by the work of Alvan et al. (2013), who studied corotation resonances for gravity waves propagating in stratified and vertically shear flows.
Based on the analytical results, we discussed possible applications to stellar and planetary interiors in Sect. 5. We have estimated the rate of differential rotation in solarlike stars using the shear contrast (the rotation difference between the pole and theequator) and in giant gaseous planets through the local Rossby number (the ratio between the shear and the rotation frequency incylindrical coordinates). We find that throughout the lifetimes of K and Gtype stars, as well as for Jupiter and Saturn at the present time, a regime where inertial waves are strongly damped is largely preferred in the convective envelope of these objects. Similar conclusions were found by Alvan et al. (2013) for internal gravity waves through critical levels in the core of solarlike stars.
It is interesting to discuss different regimes of the wave transmission in terms of angular momentum transfer for cases of strong damping and wave amplification. First, we have to underline that the theoretical analysis presented in Sect. 3 (using the Frobenius method) does not adequately characterise wave (over)reflection as similarly observed in the numerical section (Sect. 4). What we can access is the wave action flux on either side of the critical level. The analysis of the changes in the wave action flux across the critical level allows us to understand whether energy is deposited into or extracted from the mean flow, in line with the work carried out by Miles (1961), Booker & Bretherton (1967), Grimshaw (1975a), and Lindzen & Barker (1985). We have demonstrated that, in the presence of a locally conical differential rotation, a valve effect can be found for critical levels other than the corotation, analogous to the results of Acheson (1972) and Grimshaw (1975a) for hydromagnetic and gravitoinertial waves with a vertical shear. For these peculiar critical levels, waves can be attenuated when going in one direction, mainly featured by the sign of the rotation components in the box, or fully transmitted when going in the other direction. For cylindrical differential rotation, we have found a criterion analogous to the MilesHoward theorem for stratified shear flows (Miles & Howard 1964), which, for inertial waves such as those studied in this work, can be formulated as: (99)
The above criterion depends on the shear Rossby number R_{o} = U′∕(2Ω) and the vertical (k_{z}) and longitudinal (k_{x}) wavenumbers. This last point is an important difference with the MilesHoward criterion, which does not involve wavenumbers. For this reason, the analogy between Eq. (99) and the MilesHoward stability criterion must be taken with care. We also stress that Eq. (99) is very different from the Rayleigh’s inflection point theorem for Rossby waves, which are a subclass of inertial waves, when neglecting the vertical or the radial perturbed velocity (Bretherton 1966; Lindzen & Tung 1978). The MilesHoward criterion allows us to distinguish between critical levels where strong wave attenuation is expected for Ri > 1∕4 (where Ri is the Richardson number) and those where overreflection and overtransmission can lead to potential shear instabilities for Ri < 1∕4. Lindzen (1988) warns, however, that overreflection and overtransmission are a necessary but not sufficient condition for shear instability. Such amplifications leading to the instability require peculiar conditions in a threelayer model, where the shear zone that features the critical level is surrounded by a region of incoming propagating waves and by a ‘sink’ zone to force waves to cross the evanescent shear zone. Special boundary conditions are necessary for the wave to successively return to the critical level and induce wave amplitude growth. Recent studies (see e.g. Carpenter et al. 2012, for a review) have revisited instabilities in stratified shear flows by studying multiple counterpropagating waves that can interact with one another to grow in amplitude with time (with conditions such as phaselocking). A parallel was drawn between overreflection mechanisms and interacting counterpropagating waves by Harnik & Heifetz (2007) to describe baroclinic instabilities for Rossby waves.
Contrary to what our results predict, Baruteau & Rieutord (2013) did not observe any instabilities of inertial waves when using cylindrical differential rotation. Several reasons can be put forward to explain this discrepancy, such as boundary conditions (as discussed in the previous paragraph) or the values of the shear and horizontal wavenumbers, since the inertial waves may not be in the regime which allows instabilities, according to the criterion in Eq. (99), which needs to be adapted further to the global cylindrical geometry used in the work of Baruteau & Rieutord (2013). We stress that when exploring different powerlaws for the mean flow profiles, overreflection was not retrieved for a nonlinear mean flow in cylindrical differential rotation. Guenel et al. (2016a,b) did observe instabilities with conical differential rotation, but only for sufficiently low viscosities, whereas our study shows little dependence on the friction and rather highlights possible overtransmission for nonlinear flows. Lastly, we underline that a temporal analysis on the growth rate of perturbations should be undertaken to unravel instabilities, which has not been performed in this paper but has been in other separate papers (Park et al. 2020, 2021).
This ab initio analytical study is thus a first step towards understanding how inertial waves interact with a mean flow subject to latitudinal differential rotation at critical levels in the context of tidal dissipation in differentially rotating stars and planets. Possible feedbacks of the perturbed wave on equilibrium quantities and the mean flow are not taken into account in this study, nor are nonlinearities in the perturbed hydrodynamical wave equations. Nonetheless, they should be considered in future studies since Barker & Ogilvie (2010) and Baruteau & Rieutord (2013) suggested important nonlinear effects for inertial waves at corotation. Finally, magnetism may also play an important role in dissipating or redistributing angular momentum at critical levels through magnetic stresses (e.g. Wei 2016, 2018; Lin & Ogilvie 2018; Astoul et al. 2019).
Acknowledgements
We would like to thank the anonymous referee, as well as A. Barker, for the helpful comments and suggestions regarding our work. A.A., J.P., and S.M. acknowledge funding by the European Research Council through the ERC grant SPIRE 647383. The authors acknowledge the PLATO CNES funding at CEA/IRFU/DAp and IRAP. This research hasmade use of NASA’s Astrophysics Data System and of the software MATLAB version R2018a. Finally, we thank Quentin André and Jéremy Ahuir for fruitful discussions on the details of the analytical model.
Appendix A Derivation of the wave propagation equation
In this section, we detail the derivation of the ODE (cf. Eqs. (22) and (23)) in a dissipative medium and for forced inertial waves. Writing the perturbation variables as wavelike functions in x and z directions, the hydrodynamic equations, (Eqs. (16), (17), and (18)), are, respectively:
with s = σ + iσ_{f} a complex frequency that includes the Rayleigh friction frequency σ_{f}. We have removed the symbol ′ in each perturbed quantity, and it now refers to the derivative with respect to y. The projection of the curl of the perturbed momentum Eq. (A.1) on the (e_{x}, e_{y}, e_{z}) basis is:
We note that the term is of order ϵ and thus neglected. This enables us to define and use the quantity Π = p∕ρ_{0}, which we call reduced pressure in the present study. Using the continuity equation, the linear combination ik_{z} (A.4) − ik_{x}(A.6) reads: (A.7)
where . Thus, gives a secondorder ordinary differential equation on the latitudinal velocity v, with a source term S: (A.8)
Coefficients are written as: (A.9)
Appendix B Derivation of the polarisation relations
The NavierStokes equation with Rayleigh friction, without forcing, and projected onto the Cartesian basis reads:
where we use Eq. (A.3) to replace the perturbed density. In order to obtain an equation for the reduced pressure and the latitudinal velocity perturbations only, we applied the following linear combination: (B.4)
which yields, using the continuity Eq. (A.2): (B.5)
To get the perturbed vertical velocity, one can use the linear combination , which yields, after some algebra: (B.6)
Furthermore, gives the perturbed longitudinal velocity: (B.7)
Appendix C The WKBJ approximation
Rapidly oscillating solutions are often studied within the WKBJ approximation (Press 1981). At the pole, the solution of the ODE (60) in this approximation takes the form: (C.1)
where A and B are the complex amplitudes of the wave function Ψ and κ(y) is the complex potential associated with the ODE. We can determine the validity domain of this approximation in the same way Alvan et al. (2013) did. In the WKBJ approximation, Ψ satisfies (C.2)
where for a constant shear. The WKBJ approximation is valid provided that: (i) R_{o} < 1 so that R > 0, and (ii) is negligible in front of . Then, we introduce the Liouville transformation: (C.3)
Equation (C.2) thus becomes: (C.5)
The WKBJ approximation states that Φ≪ 1. Given the definition of f, this leads to the validity condition R≫ 1∕4.
Appendix D Analytical properties in the noshear regions (zones I and III)
D.1 Wavelike solutions
The ordinary differential equation without shear can be written as: (D.1)
For the following, we introduce s_{1} = ω + iσ_{f} in zone I, and s_{3} = ω − k_{x}Λ + iσ_{f} in zone III. In the noflow region I, the analytic solution of Eq. (D.1) can be written in the form: (D.2)
which is the sum of an incident wave and a reflected wave of amplitudes A_{I} and A_{R} and wavenumbers: (D.3)
respectively. These relationships are given by the dispersion relation satisfied by the incident and reflected waves: (D.4)
Similarly, we write the latitudinal velocity in the uniform mean flow region III as: (D.5)
where A_{T} is the amplitude of the transmitted wave and (D.6)
is its wavenumber, also given by the dispersion relation of the transmitted wave. It is important to stress that all varying parameters in the model, such as Λ, k_{x}, k_{z}, and ω, have been chosen such that k_{I}, k_{R}, and k_{T} are not complexwithout friction. This prevents widely diverging waves in zones I or III.
Requirements on ω for the existence of critical points inside the shear region II depending on the sign of the shear Λ.
D.2 Wave action flux
As we did in Sect. 3.2.5, one can also determine the wave action flux (Eq. (41)) in region I, (D.7)
with σ_{3} = ω − k_{x}Λ, and where we use the reduced pressure for all three waves: (D.9)
It is noteworthy that in the noshear regions the latitudinal flux of energy is preserved (see also Eliassen & Palm 1961), and so is the wave action flux. Indeed, the characteristic frequencies (ω and σ_{3}) in both zones (I andIII, respectively) are constants. Furthermore, as we stated in Sect. 3.2.5, the direction of the group velocity is given by . By consequence, the direction of energy propagation is constrained by the sign of the Dopplershifted frequency σ_{3} in zone III, and by the sign of inertial frequency ω in zone I. Considering ω > 0, the incident wave is properly named since the group velocity is positive; likewise, the reflective wave has negative group velocity and is thus moving downwards. It is a little more complicated for the socalled transmitted wave. The sign of σ_{3} is directly related to the presence or the absence of a critical level inside zone II.
Table D.1 summarises which frequency has to be exited so that waves can meet one or several critical points in the shear region of range 0 < y < 1 for a linear mean flow profile.
Thus, if corotation is met in zone II with Λ > 0, we automatically have a negative Dopplershifted frequency in zone III (i.e. σ_{3} < 0). Therefore, one must choose the − sign in the expression of k_{T} (Eq. (D.6)) in order to construct a wave that moves away from the critical level. The same reasoning applies for the critical level . It is trickier for the singularity because , so σ_{3} can be either positive or negative in the interval [0, π∕2] depending on the values of and k_{x}Λ.
References
 Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions (New York: Dover Publications Inc.) [Google Scholar]
 Acheson, D. J. 1972, J. Fluid Mech., 53, 401 [NASA ADS] [CrossRef] [Google Scholar]
 Alvan, L., Mathis, S., & Decressin, T. 2013, A&A, 553, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Amard, L., Palacios, A., Charbonnel, C., et al. 2019, A&A, 631, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 André, Q., Barker, A. J., & Mathis, S. 2017, A&A, 605, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Andrews, D. G., & McIntyre, M. E. 1978, J. Fluid Mech., 89, 647 [NASA ADS] [CrossRef] [Google Scholar]
 Astoul, A., Mathis, S., Baruteau, C., et al. 2019, A&A, 631, A111 [CrossRef] [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Le PoncinLafitte, C., & Mathis, S. 2014, SF2A2014: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. J. Ballet, F. Martins, F. Bournaud, R. Monier, & C. Reylé, 199 [Google Scholar]
 Auclair Desrotour, P., Mathis, S., & Le PoncinLafitte, C. 2015, A&A, 581, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barker, A. J., & Ogilvie, G. I. 2009, MNRAS, 395, 2268 [Google Scholar]
 Barker, A. J., & Ogilvie, G. I. 2010, MNRAS, 404, 1849 [NASA ADS] [Google Scholar]
 Barnes, J. R., Collier Cameron, A., Donati, J. F., et al. 2005, MNRAS, 357, L1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Barnes, J. R., Jeffers, S. V., Haswell, C. A., et al. 2017, MNRAS, 471, 811 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., & Masset, F. 2008, ApJ, 672, 1054 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., & Rieutord, M. 2013, J. Fluid Mech., 719, 47 [Google Scholar]
 Bazot, M., Benomar, O., ChristensenDalsgaard, J., et al. 2019, A&A, 623, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beaudoin, P., Strugarek, A., & Charbonneau, P. 2018, ApJ, 859, 61 [Google Scholar]
 Benbakoura, M., Réville, V., Brun, A. S., Le PoncinLafitte, C., & Mathis, S. 2019, A&A, 621, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Benomar, O., Bazot, M., Nielsen, M. B., et al. 2018, Science, 361, 1231 [Google Scholar]
 Bolmont, E., & Mathis, S. 2016, Celest. Mech. Dyn. Astron., 126, 275 [Google Scholar]
 Bolmont, E., Gallet, F., Mathis, S., et al. 2017, A&A, 604, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Booker, J. R., & Bretherton, F. P. 1967, J. Fluid Mech., 27, 513 [NASA ADS] [CrossRef] [Google Scholar]
 Bretherton, F. P. 1966, Quart. J. R. Meteorol. Soc., 92, 325 [Google Scholar]
 Bretherton, F. P., & Garrett, C. J. R. 1968, Proc. R. Soc. London Ser. A, 302, 529 [Google Scholar]
 Broad, A. S. 1995, Quart. J. R. Meteorol. Soc., 121, 1891 [Google Scholar]
 Brun, A. S., García, R. A., Houdek, G., Nandy, D., & Pinsonneault, M. 2015, Space Sci. Rev., 196, 303 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Strugarek, A., Varela, J., et al. 2017, ApJ, 836, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Bryan, G. H. 1889, Phil. Trans. R. Soc. London Ser. A, 180, 187 [Google Scholar]
 Carpenter, J. R., Tedford, E. W., Heifetz, E., & Lawrence, G. A. 2012, Appl. Mech. Rev., 64, 061001 [Google Scholar]
 Cartan, E. 1922, Bull. Sci. Math., 46, 317 [Google Scholar]
 Damiani, C., & Mathis, S. 2018, A&A, 618, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Debras, F., & Chabrier, G. 2019, ApJ, 872, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Decressin, T., Mathis, S., Palacios, A., et al. 2009, A&A, 495, 271 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duguid, C. D., Barker, A. J., & Jones, C. A. 2020, MNRAS, 491, 923 [Google Scholar]
 Eliassen, A., & Palm, E. 1961, On the Transfer of Energy in Stationary Mountain Waves (Det Norske VidenskapsAkademi i Oslo. Geofysiske publikasjoner) (I kommisjon hos Aschehoug) [Google Scholar]
 Favier, B., Barker, A. J., Baruteau, C., & Ogilvie, G. I. 2014, MNRAS, 439, 845 [NASA ADS] [CrossRef] [Google Scholar]
 Ford, E. B., & Rasio, F. A. 2006, ApJ, 638, L45 [NASA ADS] [CrossRef] [Google Scholar]
 Fuller, J., Luan, J., & Quataert, E. 2016, MNRAS, 458, 3867 [NASA ADS] [CrossRef] [Google Scholar]
 Galanti, E., Kaspi, Y., Miguel, Y., et al. 2019, Geophys. Res. Lett., 46, 616 [NASA ADS] [CrossRef] [Google Scholar]
 Gallet, F., Bolmont, E., Mathis, S., Charbonnel, C., & Amard, L. 2017, A&A, 604, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gallet, F., Bolmont, E., Bouvier, J., Mathis, S., & Charbonnel, C. 2018, A&A, 619, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 García, R. A., TurckChièze, S., JiménezReyes, S. J., et al. 2007, Science, 316, 1591 [Google Scholar]
 Gastine, T., Wicht, J., & Aurnou, J. M. 2013, Icarus, 225, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Gerkema, T., Zimmerman, J. T. F., Maas, L. R. M., & van Haren, H. 2008, Rev. Geophys., 46, RG2004 [Google Scholar]
 Gliatto, M. T., & Held, I. M. 2020, J. Atm. Sci., 77, 859 [Google Scholar]
 Goldreich, P., & Nicholson, P. D. 1989, ApJ, 342, 1075 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857 [Google Scholar]
 Greenspan, H. P. 1969, The theory of rotating fluids (Cambridge University Press) [Google Scholar]
 Grimshaw, R. H. J. 1975a, J. Fluid Mech., 70, 287 [Google Scholar]
 Grimshaw, R. H. J. 1975b, Tellus Ser. A, 27, 351 [Google Scholar]
 Grimshaw, R. 1979, Geophys. Astrophys. Fluid Dyn., 14, 303 [Google Scholar]
 Guenel, M., Baruteau, C., Mathis, S., & Rieutord, M. 2016a, A&A, 589, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Guenel, M., Mathis, S., Baruteau, C., & Rieutord, M. 2016b, ArXiv eprints [arXiv:1612.05071] [Google Scholar]
 Harnik, N., & Heifetz, E. 2007, J. Atm. Sci., 64, 2238 [Google Scholar]
 Hut, P. 1980, A&A, 92, 167 [NASA ADS] [Google Scholar]
 Jackson, B., Greenberg, R., & Barnes, R. 2008, ApJ, 678, 1396 [Google Scholar]
 Jones, W. L. 1967, J. Fluid Mech., 30, 439 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, W. L. 1968, J. Fluid Mech., 34, 609 [Google Scholar]
 Jouve, L., & Ogilvie, G. I. 2014, J. Fluid Mech., 745, 223 [Google Scholar]
 Kaspi, Y., Guillot, T., Galanti, E., et al. 2017, Geophys. Res. Lett., 44, 5960 [Google Scholar]
 Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structure and Evolution (Berlin: Springer) [Google Scholar]
 Lai, D. 2012, MNRAS, 423, 486 [Google Scholar]
 Lainey, V., Jacobson, R. A., Tajeddine, R., et al. 2017, Icarus, 281, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Lainey, V., Casajus, L. G., Fuller, J., et al. 2020, Nat. Astron., 4, 1053 [CrossRef] [Google Scholar]
 Latter, H. N., & Balbus, S. A. 2009, MNRAS, 399, 1058 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, Y., & Ogilvie, G. I. 2018, MNRAS, 474, 1644 [NASA ADS] [CrossRef] [Google Scholar]
 Lindzen, R. S. 1988, Pure Appl. Geophys., 126, 103 [Google Scholar]
 Lindzen, R. S., & Barker, J. W. 1985, J. Fluid Mech., 151, 189 [NASA ADS] [CrossRef] [Google Scholar]
 Lindzen, R. S., & Tung, K. K. 1978, J. Atm. Sci., 35, 1626 [Google Scholar]
 Luan, J., Fuller, J., & Quataert, E. 2018, MNRAS, 473, 5002 [CrossRef] [Google Scholar]
 Maas, L. R. M., & Lam, F. P. A. 1995, J. Fluid Mech., 300, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Maslowe, S. A. 1986, Ann. Rev. Fluid Mech., 18, 405 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Mathis, S. 2009, A&A, 506, 811 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S. 2015, A&A, 580, L3 [Google Scholar]
 Mathis, S. 2019, EAS Publ. Ser., 82, 5 [Google Scholar]
 Mathis, S., & de Brye, N. 2012, A&A, 540, A37 [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., Prat, V., Amard, L., et al. 2018, A&A, 620, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Miles, J. W. 1961, J. Fluid Mech., 10, 496 [NASA ADS] [CrossRef] [Google Scholar]
 Miles, J. W., & Howard, L. N. 1964, J. Fluid Mech., 20, 331 [Google Scholar]
 Militzer, B., Wahl, S., & Hubbard, W. B. 2019, ApJ, 879, 78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morse, P. M., & Feshbach, H. 1953, Methods of Theoretical Physics (Feshbach Publishing) [Google Scholar]
 Ogilvie, G. I. 2009, MNRAS, 396, 794 [NASA ADS] [CrossRef] [Google Scholar]
 Ogilvie, G. I. 2014, ARA&A, 52, 171 [Google Scholar]
 Ogilvie, G. I., & Lesur, G. 2012, MNRAS, 422, 1975 [NASA ADS] [CrossRef] [Google Scholar]
 Ogilvie, G. I., & Lin, D. N. C. 2004, ApJ, 610, 477 [Google Scholar]
 Ogilvie, G. I., & Lin, D. N. C. 2007, ApJ, 661, 1180 [NASA ADS] [CrossRef] [Google Scholar]
 Park, J., Prat, V., & Mathis, S. 2020, A&A, 635, A133 [EDP Sciences] [Google Scholar]
 Park, J., Prat, V., Mathis, S., & Bugnet, L. 2021, A&A, 646, A64 [EDP Sciences] [Google Scholar]
 Press, W. H. 1981, ApJ, 245, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Rieutord, M. 2015, Fluid Dynamics: An Introduction (Berlin: Springer) [Google Scholar]
 Rieutord, M., & Valdettaro, L. 1997, J. Fluid Mech., 341, 77 [Google Scholar]
 Rieutord, M., & Valdettaro, L. 2010, J. Fluid Mech., 643, 363 [NASA ADS] [CrossRef] [Google Scholar]
 Rieutord, M., Georgeot, B., & Valdettaro, L. 2001, J. Fluid Mech., 435, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Ringot, O. 1998, PhD thesis (Paris 7), France [Google Scholar]
 Schmid, P., Henningson, D., & Jankowski, D. 2002, Appl. Mech. Rev., 55, B57 [Google Scholar]
 Schou, J., Antia, H. M., Basu, S., et al. 1998, ApJ, 505, 390 [Google Scholar]
 Shampine, L. F. & Reichelt, M. W. 1997, SIAM J. Sci. Comp., 18, 1 [Google Scholar]
 Thompson, M. J., ChristensenDalsgaard, J., Miesch, M. S., & Toomre, J. 2003, ARA&A, 41, 599 [NASA ADS] [CrossRef] [Google Scholar]
 Tsang, D., & Lai, D. 2009, MNRAS, 400, 470 [NASA ADS] [CrossRef] [Google Scholar]
 Varela, J., Strugarek, A., & Brun, A. S. 2016, Adv. Space Res., 58, 1507 [NASA ADS] [CrossRef] [Google Scholar]
 Watts, A. L., Andersson, N., & Williams, R. L. 2004, MNRAS, 350, 927 [CrossRef] [Google Scholar]
 Wei, X. 2016, ApJ, 828, 30 [Google Scholar]
 Wei, X. 2018, ApJ, 854, 34 [Google Scholar]
 Yamanaka, M., & Tanaka, H. 1984, J. Meteorol. Soc. Jpn., 62, 1 [Google Scholar]
 Zahn, J. P. 1966, Ann. Astrophys., 29, 313 [NASA ADS] [Google Scholar]
 Zahn, J. P. 1977, A&A, 500, 121 [NASA ADS] [Google Scholar]
 Zahn, J. P., Talon, S., & Matias, J. 1997, A&A, 322, 320 [NASA ADS] [Google Scholar]
We note that the velocity of energy density, V_{g}, in the latitudinal direction has been named the ‘group velocity’ for obvious physical reasons, but it differs from the group velocity defined in Sect. 3.2, which, unlike V_{g}, depends on latitudinal and vertical wavenumbers.
All Tables
Summary of analytical results at each critical level when the box is inclined or not with respect to the rotation axis, in the [0, π∕2] quadrant, and for positive wavenumbers.
Reflection R and transmission T coefficients and the ratio of the wave action flux above and below a critical level for a linear, square, and cubic mean flow (n = 1, 2, and 3, respectively).
Requirements on ω for the existence of critical points inside the shear region II depending on the sign of the shear Λ.
All Figures
Fig. 1 Sketch of the local Cartesian box in the convective region of a lowmass star or giant planet. Global spherical coordinates, including the depth r_{0}, the inclination of the box θ_{0}, and the colatitude θ of a point of interest M inside the box, are shown to facilitate the analogy between the spherical and the Cartesian geometries. 

In the text 
Fig. 2 Transmission rate of a wave passing through any of the critical levels defined by as a functionof the absolute value of the shear Rossby number R_{o} = U′ (where U′ is scaled by 2Ω_{0}) and the ratioof the horizontal wavenumbers α_{k} = k_{z}∕k_{x} for a colatitude of the box θ_{0} = 10° (left panel) and θ_{0} = 80° (right panel). 

In the text 
Fig. 3 Diagram showing how compares to 1∕4 for two positions of the box (θ_{0} = 0 and π), as well as for a range of Rossby numbers R_{o} and a range of the vertical to longitudinal wavenumber ratio α_{k}. The solid and dashed black lines mark where at the north and south poles, respectively. The purple domain shows where , and the white region where , regardless of whether the box is at the north or the south pole. In the dark red region, and and vice versa in the orange region. 

In the text 
Fig. 4 Transmission rate T_{F} of the wave action flux across the corotation when the box is located at the pole for R > 1∕4. It is displayed against the Rossby number R_{o} and the ratioof the horizontal wavenumbers α_{k} = k_{z}∕k_{x}. The forbidden region where R < 1∕4 is shown in white, and the dark red cone corresponds to values of T_{F} that are lower than 10^{−9}. 

In the text 
Fig. 5 Mean flow profiles used in the threezone numerical model against y. In the noshear regions I (y < 0) and III (y > 1), the mean flow is uniform and set respectively to U = 0 and U = Λ. The shearedregion II (greyshaded) can have a linear, square, or cubic mean flow profile. 

In the text 
Fig. 6 Transmission coefficient (T, left panel) and reflection coefficient (R, right panel) when the box is at the pole. The coefficients are plotted in the regime R < 1∕4 (a possibly unstable case, see Sect. 3.4.4) as a function of the Rossby number R_{o} and the inertial frequency ω. The hatched areas do not feature critical points and correspond to regions where ω > k_{x} R_{o} in our peculiar geometry (see Appendix D.2 for this particular matter). Vertical and longitudinal wave numbers are fixed: k_{x} = 0.1 and k_{z} = 0.1. Moreover, the contours that correspond to coefficients R and T that are equal to one are indicated by solid black lines. Crosses mark the set of parameters used in Fig. 9for the analysis of the behaviour of the velocity in the threelayer model. 

In the text 
Fig. 7 Same as Fig. 6, but for R > 1∕4. When R_{o} < 0 in these panels, hatched areas mark instances where ω > 0 (see Appendix D.2). 

In the text 
Fig. 8 Ratio of the wave action flux above and below the critical level against the Rossby number and the inertial frequency for R >1∕4 (left panel) and R < 1∕4 (right panel) when the box is at the pole. As in Fig. 7, k_{x} = k_{z} = 0.1 and hatchedzones represent areas without critical levels. Again, crosses indicate the set of parameters chosen to analyse the behaviour of the velocity in the threelayer model. 

In the text 
Fig. 9 Numerical outputs of the threelayer model. Top: real part of the latitudinal velocity v against y. The quantities v_{num}, v_{I}, v_{R}, v_{T}, v_{F}, and v_{W} are the numerical, incident, reflected, transmitted, firstorder Frobenius, and Whittaker velocities, respectively. Bottom: wave action flux against y. The quantities , , , and are the numerical, transmitted, incident and reflected, and Frobenius wave action fluxes, respectively. For all panels, θ_{0} = 0, the mean flow is linear in the greyshaded shear regions, and the critical level is marked by dashed lines. The horizontal wave numbers are set to k_{x} = k_{z} = 0.1. From left to right: (i) ω = 0.02 and R_{o} = 0.3 (R > 1∕4), (ii) ω = 0.002 and R_{o} = 0.8 (R > 1∕4), and (iii) ω = 0.09 and R_{o} = 1.8 (R < 1∕4). 

In the text 
Fig. 10 Same quantities as in Fig. 9 (when the box is at the pole), but for a box tilted by 10 degrees relative to the pole and for different values of the shear, wavenumbers, and inertial frequencies. We note that, unlike Fig. 9, we do not have analytical solutions in the shear region. In all panels, the horizontal wavenumbers are set to k_{x} = k_{z} = 1 and the shear is R_{o} = 0.3. From left to right: inertial frequencies set to ω = 0.31, ω = 0.16, and ω = 0.1. Third panels: we indicate y_{0} and y_{−} with dasheddotted and vertical dashed lines, respectively. 

In the text 
Fig. 11 Absolute values of the shear contrast normalised by the rotation at the pole against the age of K and G spectral type stars. Solid lines feature 0.3 < R_{of} < 0.9, dashed lines R_{of} ≲ 0.3, and transparent lines R_{of} > 1. 

In the text 
Fig. 12 Rossby numbers for Jupiter and Saturn as a function of the axial distance r. 

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.