The complex interplay between tidal inertial waves and zonal flows in differentially rotating stellar and planetary convective regions I. Free waves

Quantifying tidal interactions in close-in two-body systems is of prime interest since they have a crucial impact on the architecture and on 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, like differential rotation which exists in the outer convective enveloppe of solar-like stars and giant gaseous planets. In particular, tidal waves may strongly interact with zonal flows at the so-called corotation resonances, where the wave's Doppler-shifted frequency cancels out. We aim to provide a deep physical understanding of the dynamics of tidal inertial waves at corotation resonances, in the presence of differential rotation profiles typical of low-mass stars and giant planets. By developping an inclined shearing box, 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 one-dimensional three-layer numerical model. We find that inertial waves can be either 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 to the fluid when damped at critical levels, or they can extract action flux to the fluid when amplified at critical levels. Both situations could lead to significant angular momentum exchange between the tidally interacting bodies.


Introduction
Tidal interactions are known to drive the late evolution of short-period 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 low-mass host stars and giant planets can modify the spin of the tidally perturbed body, the orbital period, and the spin-orbit 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 1966Zahn , 1977;;Duguid et al. 2020).For coplanar and circular systems, inertial waves are excited so long as the companion orbits beyond half its co-rotation radius (the orbit where the host's rotation frequency is equal to the mean motion).Low-mass 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, A&A 647, A144 (2021) also called attractors of characteristics (Maas & Lam 1995), that are confined within the convective envelope (see also Rieutord & Valdettaro 1997).With a non-zero 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 gravito-inertial 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 non-wave-like 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 (Auclair-Desrotour 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 frequency-averaged 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 two-body systems differently, as postulated, for instance, by Lai (2012) to explain the survival of hot Jupiters with completely damped spin-orbit 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 low-mass stars and giant gaseous planets.The Sun's surface rotates in ∼25 days at the equator versus ∼35 days near the poles, and a latitude-dependent 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 solar-like stars (Benomar et al. 2018).Essentially, differential rotation in low-mass stars depends on the effective temperature (Barnes et al. 2005(Barnes et al. , 2017) ) and seems to be more important when the convective envelope is thinner.Solar-like rotation profiles and anti-solar-like rotation profiles (with faster poles and a slower equator) are expected in G-and K-type stars depending on their rotation rates, based on three-dimensional 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 Wentzel-Kramers-Brillouin-Jeffreys (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 solid-body 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 solid-body 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 co-rotation resonances where the Doppler-shifted 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 K-type) stars.They also confirmed the existence of unstable inertial modes (i.e.modes with positive growth rates) at co-rotation 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 co-rotation resonances, particularly at low viscosities.Favier et al. (2014) also studied tidally forced inertial waves, but through non-linear 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 co-rotation 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 late-type 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 co-rotation 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 A144, page 2 of 25 A. Astoul et al.: The complex interplay between tidal inertial waves and zonal flows 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 stratified vertical 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 the wave 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 co-rotation 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 gravito-inertial 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 large-scale tidal perturbations and convective motions.In our model, we focus on the latitudinal differential rotation of the mean flow, 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 second-order 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 interpret energy flux exchanges between the waves and the mean flow.We use, in Sect.4, a three-layer numerical model to test our analytical predictions at critical levels.Frictional damping is included, and non-linear mean flow profiles are also used.Astrophysical applications with implications for low-mass 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.

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 gravito-inertial 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 two-dimensional 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 where Ω 0 is the rotation frequency of the star at the pole and f 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 f -plane 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 GM/R 3 , 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 small-scale 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.

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: where e ϕ is the azimuthal unit vector and r and θ are the radius and co-latitude, respectively.We introduce u 0 = r sin θΩ 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 co-latitude θ and at the pole).The shear contrast is positive for the Sun since the equator rotates faster than the pole, and negative for antisolar-like 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 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 (one-dimensional) 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. 2004Mathis et al. , 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 solid-body rotation used by Guenel et al. (2016a) was: where χ is the magnitude of the shear between the equator and the pole.Performing a second-order Taylor expansion around a fixed co-latitude θ 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 (5) 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): 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 solar-like 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 y-profile 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).

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 d dt = ∂ dt + u • ∇ is the total derivative operator.
All variables are then linearised at the first order: Zeroorder terms correspond to background equilibrium quantities, and first-order terms represent the leading perturbation.The local velocity, density, and pressure are therefore written as: where u = (u, v, w) in the local Cartesian basis.We have introduced the dimensionless parameter : 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: A144, page 4 of 25 A. Astoul et al.: The complex interplay between tidal inertial waves and zonal flows with n the unit vector parallel to the rotation axis.Projecting Eq. ( 12) into Cartesian coordinates, one can derive: 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 thermal-wind 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 Brunt-Väisälä frequency is where we introduce the dimensionless number F = gL/c 2 s , which is small when filtering acoustic waves.Consequently, the curl of Eq. ( 12) gives where we neglect the second-order 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 Brunt-Vä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 thermal-wind balance.

Equilibrium state of the background flow
It is worthwhile discussing the choice of keeping buoyancy forces in the zero-order momentum equation.Without gravitational forces, the momentum equation for mean dimensional variables is written as a geostrophic balance: This balance satisfies the Taylor-Proudman theorem (Rieutord 2015), namely the geostrophic flow is independent of the coordinate parallel to the rotation axis.When taking the x-axis (the only non-zero) projection of the curl of this equation, one obtains the following relationship: Without the vertical stratification embodied by the Brunt-Väisälä frequency, nor latitudinal stratification, the equilibrium of a y-dependent mean flow is thus not ensured for an incompressible fluid.An alternative to conserve the equilibrium without stratification would be to consider a z-dependence 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 f = 0 at the poles, the latitudinal stratification term will not appear in the perturbed fluid equations (as we can see from Eq. ( 18)).

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 non-dissipative fluid at various co-latitudes.For this purpose, we consider perturbations q in the normal mode 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.

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 second-order ODE for v: 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: where k ⊥ = k 2 x + k 2 z is the absolute wavenumber in the direction perpendicular to the y direction and σ = ω − k x U is the (dimensionless) Doppler-shifted 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 A144, page 5 of 25 A&A 647, A144 (2021) levels (see e.g.Bretherton 1966;Grimshaw 1975a).The critical level where the Doppler-shifted frequency equals zero (i.e.σ = 0) can be met when the mean flow matches the local phase velocity, and this is also known as 'co-rotation 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 co-latitudes other than the poles, the critical levels come in three flavours, the co-rotation σ = 0 and two other critical levels that are defined, in our model, by σ = ± f (we recall that f 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 Doppler-shifted frequency at critical levels other than the co-rotation resonance equals ±2Ω v , where Ω v is the vertical component of the rotation vector.

Dispersion relation, group, and phase velocities
The three-dimensional dispersion relation is a fourth-order equation in Doppler-shifted frequency when injecting wave-like 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 short-wavelength approximation in the meridional plane, as in Baruteau & Rieutord (2013) and Guenel et al. (2016a).This involves keeping only the second-order derivatives in the y and z directions, and it reduces the relation dispersion to a second-order equation when injecting plane wave-like solutions.In the local meridional plane, the differential equation reduces to a Poincaré-like equation: 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 ( f = 0 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 Doppler-shifted frequency σ: where ||k|| = k 2 y + k 2 z is the norm of the wave vector in the meridional plane (e.g. for fixed k x ), as in Baruteau & Rieutord (2013).Compared to solid-body rotation (see e.g.Rieutord 2015), an additional term (k 2 z f U ) 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 We can then specify the phase velocity in the meridional plane: In the same way, we can derive the expression for the group velocity in the meridional plane: We note that without differential rotation, the group velocity reduces to its well-known expression for solid-body rotation (e.g.see Rieutord 2015): Moreover, as in solid-body rotation, the group velocity (Eq.( 28)) and the phase velocity (Eq.( 27)) lie in perpendicular planes: When the box is located at the north pole (θ 0 = 0 in Fig. 1), by setting κ = 1 − U , we recover 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 ).

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 u φ → 0 while |u g • e y | → ∞ and |u 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 B (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 co-rotation.Now to get σ → ± f , we either need |k y | → ∞ at fixed k z , which implies u φ → 0 and u 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 u φ • e z → 0 and u g • e y → 0, while |u φ • e y | → f /k y and |u 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 co-rotation 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 = r sin θ 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 σ = ± f 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 co-rotation 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 σ = ± f .Hence, A144, page 6 of 25 the local critical levels at σ = ± f behave somewhat similarly to the co-rotation 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 co-rotation are similar but lead to different relationships for the phase and group velocities: (i) κ → 0 (i.e.U → 1), meaning that u φ → 0 and u g → 0: The wave is totally absorbed at corotation.(ii) |k y | → ∞ at fixed k z , which implies u φ → 0 and u 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 u φ → 0, u g • e y → 0 while |u g • e z | → κ/k y : the wave energy does not cross the co-rotation in the latitudinal direction (equivalent to the vertical paths of characteristic in global cylindrical geometry, as in Baruteau & Rieutord 2013).At the north pole1 , we actually have a perfect match with the conditions given by Baruteau & Rieutord (2013) when using a cylindrical rotation profile for the mean flow.

Energetical aspects
In this section, we examine the energetic balance associated with inertial waves in our inclined shear box model, without assuming the short-wavelength approximation.This energetic balance does not include potential energy because of the adiabaticity of the convective region, but two additional terms relative to the solid-body 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 first-order definition This first allows us to express the perturbed density from Eq. ( 18) as ρ = −ρ 0 ζ f U (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: 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 wave-like form in these directions.Further assuming that the box is δ thick in the y direction, the energetic balance yields: 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: 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 P shear can also be seen as the power transferred from the mean flow to the perturbation (or conversely) by the Reynolds stress: where we use partial integration and the periodicity of perturbations in the x and z directions.At the pole, f = 0, so we recover the definition of the Reynolds stress in Miles (1961), who studied the stability of a two-dimensional stratified y-sheared flow (i.e.τ = −ρ 0 uv).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 z-sheared flows.Moreover, we emphasise that the latitudinal flux of energy pv is not conserved, even in the inviscid free-wave 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 y-derivative of the latitudinal flux is: Using the same method as Broad (1995), we multiplied the x-projection of the inviscid force-free momentum equation by ζ: By multiplying by (U − c), the latitudinal flux of energy can thus be written as: By differentiating this relationship with respect to y, and by equalising with Eq. ( 36), one can obtain: that is (U − c) dτ dy = 0 with τ the Reynolds stress.Equation ( 39) is naturally satisfied at co-rotation, 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 three-dimensional stratified shear flows, Eq. ( 39) is not vectorial, because our base flow is unidirectional.

Polarisation relations
For the forthcoming analysis, it is useful to derive expressions of the perturbed projected velocities and the perturbed reduced pressure2 Π = p/ρ 0 , namely the polarisation relations (see Appendix B for more details).In the inviscid free-wave problem, these perturbed quantities can be written in terms of the latitudinal velocity, its derivative, and the shear: Without shear and at θ 0 = 0, we recover the polarisation relations in the solid-body rotation case (see e.g.Rieutord 2015).

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): which is the latitudinal flux averaged over vertical and longitudinal wavelengths divided by the Doppler-shifted 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 A = −τ/k x .By using the expression for the perturbed reduced pressure Π derived in the previous section, the wave action flux now reads: 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 A = 0 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 z-sheared mean flows without rotation (Booker & Bretherton 1967), with rotation under the traditional approximation (Jones 1967), and with rotation under the non-traditional 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 A is a measure of wave energy through a surface (in the (z, x) plane) since pv is the energy density transported by the group velocity3 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.

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 co-rotation σ = 0 and the critical levels σ = ± f when the box is tilted (for the co-rotation when the box is at the pole, see Sect.3.4).

Critical levels at σ = ± f
In this subsection, we treat both singularities, σ = ± f , 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 σ = ± f by approximating the ODE through its first-order 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 σ = ± f .For a linear mean flow profile U = Λy, with Λ a constant, y ± are given by: Without any assumption on the mean flow profile, the first-order Taylor expansion of the ODE (Eq.( 22)) near y ± is: with where the symbol ± refers to the regular singularities4 y + and y − and U ± is U evaluated at these singularities.The Frobenius method consists in injecting the power function (y − y ± ) λ into Eq.( 44), with λ a constant to be determined (see e.g.Morse & Feshbach 1953).The corresponding indicial equation is then with solutions: Therefore, the two independent solutions of Eq. ( 44) can be written as follows: 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,± : 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 non-growing wave towards infinity.The Taylor expansion of the base flow at the first order in y − y ± gives and, by definition, we have Consequently, the solution below the critical level is unambiguous in terms of the above solution coefficients and depends on In other words, when taking y − y ± to decrease from positive to negative values, its complex argument changes continuously from 0 to − sign k x U ± π.Thus, the appropriate path for determining the branch of (y − y ± ) λ ± passes under (above) y ± as long as k x U ± > 0 (k x U ± < 0) (the same reasoning can be found in Grimshaw 1975a).Therefore, the solution on both sides of the critical level y ± is: 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 wave-like solution with the varying latitudinal wavenumber ∓(k z f /k x U ± ) log |y − y ± |.Moreover, according to Eq. ( 42), the wave action flux on either side of y ± is: The group velocity gives the direction towards which the energy is transported (we recall that V g E = pv, with V g the group velocity and E the local energy density).By consequence, sign(V g ) = sign(Aσ) = ∓ sign( f k z ) for the solution featuring the coefficient b 0 .If k z f is positive, this wave transports energy downwards (upwards) across the critical level σ = f (σ = − f ).
If k z f is negative, the wave transports energy upwards (downwards) across the critical level σ = f (σ = − f ).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 T θ 0 , where with α k = k z /k x and R o = U ± , after passing through the critical level.Such a wave will always be attenuated since T θ 0 ≤ 1.The transmission factors T θ 0 = 10 • and T θ 0 = 80 • 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 fast-rotating 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 T θ 0 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 A144, page 9 of 25 A&A 647, A144 (2021) coefficient b 0 ).The solution of coefficient a 0 is not affected by the attenuation.This is the so-called 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 magneto-gravito-inertial waves in an inviscid and compressible z-sheared fluid.

Inertial wave crossing co-rotation
We performed the same analysis as in the previous section to treat the co-rotation 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: where U 0 and U 0 are the first and second derivatives of the mean flow profile U evaluated at the critical level y 0 .The singularity at the co-rotation 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) with v 1 (y) = +∞ n = 0 a n (y − y 0 ) n+2 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 includes all the solutions of Eq. ( 56), with c 0 = b 0 and a 1 = b 1 + a 0 determined by boundary conditions, and c n, n∈N * \{1} = b n + a n−1 determined by recurrence via the expansion of Eq. ( 22) around y 0 .As a result, the wave action flux given by Eq. ( 42) becomes just below and above the co-rotation, 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 co-rotation 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 −σ 2 k 2 ⊥ 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.

Inertial waves when the box is at the poles
When the box is located at the north or south pole, f = 0 and the ODE (Eq.( 22)) is greatly simplified.For θ 0 = 0, the dimensionless wave propagation equation becomes At the south pole (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, 2Ω = (0, f , f ) and there is a constant d f /dy = β, with f the Coriolis parameter (e.g.Miles 1961;Grimshaw 1975b;Gliatto & Held 2020).However, we cannot make a direct comparison at co-rotation, because the singularity in the equations for Rossby waves and inertial waves is not of the same order.We have a second-order pole around the co-rotation, while only first-order poles are found in the aforementioned studies.In fact, Eq. ( 60) is similar to the wave equation in stratified z-sheared flows (e.g.Jones 1968).
In our polar configuration, the y-coordinate is now the axial distance, and this means that the mean flow has a cylindrical profile.Such a 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 Proudman-Taylor theorem for fast-rotating 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.

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: 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 left-hand term in the bracket is o 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: where ỹ = 2k ⊥ (y − y 0 ), A144, page 10 of 25 A. Astoul et al.: The complex interplay between tidal inertial waves and zonal flows and A and B are complex constants given by boundary conditions.The Whittaker function M 0,−µ allows a quite straightforward analytic continuation: By consequence, the solution below the critical point y = y 0 is: Although the Whittaker functions do not feature precisely as wave-like 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 which will simply be denoted by R in the following.This can drastically change the behaviour of a wave passing through the co-rotation.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 Miles-Howard theorem for stratified z-sheared 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 Brunt-Vä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 R θ 0 = 0, π 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: Solutions are then fully evanescent for such shears.One can notice that Eq. ( 67) is the same far from co-rotation 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.

Frobenius method at the pole
Though analytic solutions are known, it is still useful to determine Frobenius solutions near co-rotation for two main reasons.First, these solutions can be derived for any mean flow profile near the co-rotation.Close to co-rotation, the mean flow is approximated by a Taylor expansion at the first-order U = U 0 (y − y 0 ).Secondly, Frobenius solutions may feature wave-like forms, which is helpful for physical interpretation.Therefore, Eq. ( 60) can be written near co-rotation as where R o = U 0 and R o = U 0 .The indicial equation gives: In the two next subsections, we examine the two cases where µ is imaginary or real.

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 co-rotation.The solutions of the indicial equation can be recast as The first-order solutions to Eq. ( 68) in the vicinity of y 0 are: for a 0 , b 0 ∈ C. 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 f = 0, the wave action flux Eq. ( 42) reduces to that is, injecting the solutions on both sides of the critical level: 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, τ = − k x A. Moreover, given Eq. ( 35), the Reynolds stress in our model reads τ = − ρ 0 uv.Using the polarisation relations for u, we recover the wave action flux in Eq. ( 72).
The pre-factor 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 Doppler-shifted frequency satisfies sign(σ) = − ς, as was the case for the co-rotation 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 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 ).

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 wave-like functions.The exponential form of solutions for R < 1/4 near the critical level reads and makes this region fully evanescent.Furthermore, the associated wave action flux is Without knowing the direction of the wave or the energy flux, since sign(V g ) = ς sign {a 0 b * 0 } , 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 solar-like stars and evolved stars.The goal of the method is to determine the reflection and transmission coefficients in a three-zone 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 over-reflected and/or over-transmitted 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 wave-like 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 over-reflection 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 over-reflection and over-transmission requires at minimum that A144, page 12 of 25 Table 2. 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.

Box
Critical level Attenuation Amplification Notes.The table indicates whether each critical level can cause attenuation or amplification of the (upward ↑ and downward ↓) travelling wave in the y direction, which depends, notably at the pole, on R = α k (1 − R o )/R 2 o and on the wave action flux above the critical level A + .The symbol indicates that no wave action flux is carried across the critical level.Moreover, the results are analogous in the other quadrants of the spherical body (with the direction of the attenuated wave through σ = ± f varying according to sign(k z f )).
the Rossby number is not the same in the whole domain by using, for instance, a non-linear 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.

A three-zone numerical model
In order to test the analytical predictions of the previous section, we built up a three-zone 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 combination of 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 wave-like solution crosses the shear region.

Description of the model
The mean flow profile that is used in the three-zone model is illustrated in Fig. 5.The zone with shear (zone II) is surrounded by two no-shear 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 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 shear flow, 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: 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 y-values: A144, page 13 of 25 A&A 647, A144 (2021) 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 details on 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 A T = A I−R , where A T is the wave action flux of the transmitted wave, and A I−R = A I − A R , where A I and A R 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 σ 2 = f 2 , we added a small friction, σ f = 10 −8 , to our set of units.Given the boundary conditions, the numerical solver deals with two first-order ODEs for v and Π, which take the form where 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)).

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: 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 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 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 : 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 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: ) .The amplitude of the incident and reflected waves can be written, in terms of A and B, as: where The transmission and reflection coefficients are then: 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 R o = −2 ± √ 8 ≈ {−4.8, 0.8}, which delineates the two regimes, as we can see in Fig. 3. Areas that are hatched do not possess the co-rotation point σ = 0.In addition, the wavenumber k 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 over-transmission are both possible (see Fig. 6).One should notice that for R o > 1, we always have by definition 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 over-reflection, nor overtransmission, in the hatched areas where there is no co-rotation point.This highlights the essential role of the critical level in inducing the over-reflection or over-transmission 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 over-transmission 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 wave-like functions.In particular, the amplification term (y − y 0 ) 1/2 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 hand-waving explanation does not formally demonstrate the origin of this amplification, it stresses the important role of the shear-region 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.

Wave action fluxes below and above the shear region
Since v and Π are continuous at the interfaces between the shear and no-shear 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 co-rotation 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 co-rotation (see Appendix D.2 for the detailed calculation): 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 |A T /A I−R | < 1 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, A T /A I−R < 0 means that |A I | 2 > |A R | 2 since the minus sign is taken in Eq. ( 90).Therefore, no over-reflection 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 and areas where the wave is over-transmitted for , the threshold between under-and over-reflection (around R o ≈ 0.9) is the same as for the reflection factor (in the right panel of Fig. 6).For the second one (|A T /A I−R | > 1), 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 over-reflection and over-transmission phenomena in the regime where R < 1/4.

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 The quantities A num , A T , A I−R , and A F 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 grey-shaded 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 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 first-order 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, a 0 (y − yc) 1/2+i|µ| corresponding to an upward wave is sufficient to correctly fit the numerical solution, meaning that the counter-propagating wave in zone I is reflected at y = 0.However, for the middle and right panels, it is not clear whether the first-order 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 over-transmitted but not over-reflected, 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 over-reflected by a factor of ∼1.5 and over-transmitted by a factor of ∼2, in concordance with the reflection and transmission coefficients plotted in Fig. 6.The wave action flux is negative, and |A I−R | < |A T | 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 z-sheared 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.

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 A144, page 17 of 25 A&A 647, A144 (2021) -2 0 10 8 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 dashed dotted and vertical dashed lines, respectively.
equation.Instead of going for an extensive numerical investigation of the parameter space, we instead on dynamics of inertial waves in our three-layer model as they cross the critical levels σ = ± f and σ = 0. Our results are presented in Fig. 10 for a box inclination of 10 • relative to the rotation axis, a shear fixed 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 σ = + f , σ = 0, and σ = {0, − f } (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 co-rotation y = y 0 , which is in line with the theoretical analysis in Sect.3.3.2for a constant shear.However, the wave is over-reflected and over-transmitted, 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 co-rotation 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 top-left and top-right 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 co-rotation 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 σ = ± f .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 over-transmission and over-reflection), which can be explained by the exponential form of the Frobenius series.

Numerical results with a non-constant shear
The choice of a linear mean flow profile allows the resolution of the ordinary differential equation at the pole and a simpler A144, page 18 of 25 Table 3. Reflection R and transmission T coefficients and the ratio of the wave action flux above and below a critical level A T /A I−R for a linear, square, and cubic mean flow (n = 1, 2, and 3, respectively).1.5 6.4 × 10 −1 7.0 × 10 −1 1.8 1.1 9.7 × 10 −1 2.9 −1.9 −1.8
implementation of the Frobenius method.Nevertheless, the correspondence between a global and a local mean flow as presented in Sect.2.2 involves higher-order terms than a simple linear dependence.Therefore, it is important to examine the effect of different mean flow profiles with non-zero 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 Moreover, we specify that the critical level in y < 0 for even values of n is not examined given our numerical set-up 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 A T /A I−R change quite significantly between n = 1, 2, and 3.At the pole, the ('geometric') over-transmission 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 over-reflection disappear when n = 3, whereas |A T /A I−R | > 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 counter-propagating wave discussed in the previous section.At the co-rotation 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 non-linear mean flow profiles, in agreement with the wave action flux derived analytically in Eq. ( 59).We speculate that the over-transmission and over-reflection are due to the polynomial form the solutions in the shear region.These wave amplifications may be related to shear instabilities in this particular three-layer 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 co-exist, 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 co-rotation 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 co-rotation y = y 0 .

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 ∆Ω: anti-solar-like rotation for ∆Ω < 0, cylindrical rotation for |∆Ω| 1, and solar-like rotation rotation for (not-too-low) positive ∆Ω (e.g.∆Ω/Ω 0 0.3 for the Sun).Several works based on three-dimensional 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 three-dimensional simulations that determines the rotation profile of a solar-like star.This criterion is based on the fluid Rossby number R of , defined as 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 anti-solar-like rotation, (ii) 0.3 < R of < 0.9 for solar-like rotation, and (iii) R of 0.3 for cylindrical rotation.Furthermore, they introduced the shear contrast ∆Ω S at the co-latitude 30 • as since the rotation frequency is often ill-defined at low colatitudes in three-dimensional numerical simulations in spherical geometry.From their three-dimensional simulations, Brun et al. (2017) obtained the following scaling: 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 We show in Fig. 11 the quantity ∆Ω/Ω 0 , expressed in Eq. ( 95), versus the age of solar-like stars for K to G spectral types.To compute this quantity, we used grids of the one-dimensional stellar evolution code STAREVOL (see Amard et al. 2019, for details of the code).In light of Fig. 11, the stars in the pre-main sequence (age 100 Myr) exhibit cylindrical rotation as they are fast rotating.During the main sequence, stars mostly feature solar-like rotation, while anti-solar-like 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 anti-solar 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 the following relationship: 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 σ = ± f , considering Fig. 2 (we recall that at the co-rotation the wave action flux is fully transmitted).A downward (upward) propagating wave through σ = f (σ = − f ) 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 anti-solar-like differential rotation since the transmission factor T θ 0 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 σ = ± f 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 radial direction, 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.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 Ψ = √ ρr 2 v r under the anelastic approximation, where ρ is the density, r the axial distance coordinate, and v r the axial velocity.For free inertial waves, their second-order differential equation is:

Mathis
where l Θ and m denote the equatorial and azimuthal wavenumbers, respectively, σ = ω + mΩ(r) is the (linear) Doppler-shifted frequency, and κ r is the 'axial' epicyclic frequency defined as κ 2 r = 4Ω 2 + 2Ωr dΩ dr .For cylindrical differential equation, the co-rotation resonance σ = 0 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: /(2Ω c ), the local Rossby number in cylindrical coordinates.By writing α 2 k = l 2 Θ /m 2 , 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 = 0, f = −1, and U = 0).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 A144, page 20 of 25 A. Astoul et al.: The complex interplay between tidal inertial waves and zonal flows   Mathis et al. (in prep.).Moreover, the convention of a plus sign in the Doppler-shifted 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 R = α 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 co-rotation: (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 over-reflected 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 to approximate, as best as possible, the behaviour of the Legendre polynomial P m l (cos Θ) 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 co-rotation 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 so-called 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 co-rotation 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 co-rotation.

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 co-rotation 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 co-rotation resonance (characterised by a zero Doppler-shifted wave frequency) and, more broadly, at critical levels (any singularities of the governing second-order 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 low-mass 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 co-latitude 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, low-mass 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 co-rotation resonance and two other critical levels that arise from the inclination between the gravity and the rotation vectors, which are defined by a Doppler-shifted 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 co-rotation.
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 non-dissipative 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 plane-parallel 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(Grimshaw , 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 A144, page 21 of 25 A&A 647, A144 (2021) (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 three-layer numerical model that comprises a shear zone, where the critical level is located, and two surrounding shear-free 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 power-law mean flow profile and varying the friction.This conclusion is also shared by the work of Alvan et al. (2013), who studied co-rotation 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 solar-like stars using the shear contrast (the rotation difference between the pole and the equator) and in giant gaseous planets through the local Rossby number (the ratio between the shear and the rotation frequency in cylindrical coordinates).We that throughout the lifetimes of K-and G-type 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 solar-like 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 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), andLindzen &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 co-rotation, analogous to the results of Acheson (1972) and Grimshaw (1975a) for hydromagnetic and gravito-inertial 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 Miles-Howard theorem for stratified shear flows (Miles & Howard 1964), which, for inertial waves such as those studied in this work, can be formulated as: > 1/4: wave attenuation; < 1/4: possible wave over-reflection, transmission. (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 Miles-Howard criterion, which does not involve wavenumbers.For this reason, the analogy between Eq. ( 99) and the Miles-Howard 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 sub-class of inertial waves, when neglecting the vertical or the radial perturbed velocity (Bretherton 1966;Lindzen & Tung 1978).The Miles-Howard 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 over-reflection and over-transmission can lead to potential shear instabilities for Ri < 1/4.Lindzen (1988) warns, however, that over-reflection and over-transmission 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 counter-propagating waves that can interact with one another to grow in amplitude with time (with conditions such as phase-locking).A parallel was drawn between over-reflection mechanisms and interacting counter-propagating 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, over-reflection was not retrieved for a non-linear 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 over-transmission for non-linear 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(Park et al. , 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 non-linearities in the perturbed hydrodynamical wave equations.Nonetheless, they should be considered in future studies since Barker & Ogilvie (2010) and Baruteau & Rieutord (2013) suggested important non-linear effects for inertial waves at co-rotation.Finally, magnetism may also play an important role in dissipating or redistributing angular momentum at critical levels through magnetic stresses (e.g.Wei 2016Wei , 2018;;Lin & Ogilvie 2018;Astoul et al. 2019).

Fig. 1 .
Fig. 1.Sketch of the local Cartesian box in the convective region of a low-mass star or giant planet.Global spherical coordinates, including the depth r 0 , the inclination of the box θ 0 , and the co-latitude θ of a point of interest M inside the box, are shown to facilitate the analogy between the spherical and the Cartesian geometries.

Fig. 2 .
Fig. 2. Transmission rate T θ 0 of a wave passing through any of the critical levels defined by σ = ± f as a function of the absolute value of the shear Rossby number |R o | = |U | (where U is scaled by 2Ω 0 ) and the ratio of the horizontal wavenumbers α k = k z /k x for a co-latitude of the box θ 0 = 10 • (left panel) and θ 0 = 80 • (right panel).

Fig. 3 .
Fig.3.Diagram showing how R θ 0 = 0, π ≡ R 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 R θ 0 = 0, π = 1/4 at the north and south poles, respectively.The purple domain shows where R θ 0 = 0, π < 1/4, and the white region where R θ 0 = 0, π > 1/4, regardless of whether the box is at the north or the south pole.In the dark red region, R θ 0 = 0 > 1/4 and R θ 0 = π < 1/4, and vice versa in the orange region.

FFig. 4 .
Fig. 4. Transmission rate T 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 ratio of 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 .

Fig. 5 .
Fig.5.Mean flow profiles used in the three-zone numerical model against y.In the no-shear regions I (y < 0) and III (y > 1), the mean flow is uniform and set respectively to U = 0 and U = Λ.The sheared region II (grey-shaded) can have a linear, square, or cubic mean flow profile.

Fig. 6 .Fig. 7 .
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. 9 for the analysis of the behaviour of the velocity in the three-layer model.

Fig. 8 .
Fig. 8. Ratio of the wave action flux above and below the critical 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 hatched zones represent areas without critical levels.Again, crosses indicate the set of parameters chosen to analyse the behaviour of the velocity in the three-layer model.

Fig. 9 .
Fig. 9. Numerical outputs of the three-layer 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, first-order Frobenius, and Whittaker velocities, respectively.Bottom: wave action flux against y.The quantities A num , A T , A I−R , and A F 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 grey-shaded 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).

Fig. 12 .
Fig. 12. Rossby numbers R o = r 2Ω dΩ dr for Jupiter and Saturn as a function of the axial distance r.

Table 1 .
Correspondence between local and global coordinate systems.
x , e y , e z ) (e ϕ , −e θ , e r ) Conical coordinate y/r 0 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.

Table 4 .
Regime of wave transmission at co-rotation deduced from Figs.3 and 4for three orbital states of a satellite around a giant gaseous planet.
is at the poles, y is the axial distance.That is why a polar configuration of the box best reproduces the 'equatorial' model of