© A. Astoul et al. 2021

1 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 1966, 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, 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, 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 low-mass (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 co-rotation 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 case of a viscous medium. This distinction is analogous to the distinction between shear layers and attractors of characteristics that are kinds of singularities for the governing equation of inertial waves in a spherical shell. The aforementioned singularities can act very differently, with either severe absorption at the critical level (as in Booker & Bretherton 1967, for stratified vertical shear flows) or no attenuation if the wave propagates in a peculiar direction (Jones 1967; Acheson 1972; Grimshaw 1975a, for stratifiedvertical shear flows with rotation and magnetism). In other cases, a critical level may even give rise to wave amplification under certain conditions related to the first and second derivatives of the mean flow velocity (Lindzen & Tung 1978; Lindzen & Barker 1985, for barotropic and stratified shear flows, respectively). These studies all used an invariant quantity (the Reynolds stress or thewave action for rotating or magnetic flows) as a diagnostic tool to interpret the role of the critical level in terms of energy transmission and to quantify exchanges between the wave and the mean flow (Eliassen & Palm 1961; Bretherton 1966).

In light of these various studies, it is necessary to carefully consider 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 meanflow, varying the box orientation to model either cylindrical or conical rotation. The behaviour of free inertial waves in this framework is then examined near critical levels using both analytical and numerical approaches.

This paper is organised as follows. In Sect. 2, we describe the local shear model with its main assumptions and the system of governing equations. In Sect. 3, we establish a 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 interpretenergy 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.

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 r0, 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.

2 Local Cartesian model including differential rotation

2.1 Presentation of the model

The local model takes the form of an inclined sheared box, centred at a point C of a convective shell, as illustrated in Fig. 1. The inclined box model has already been used by Auclair Desrotour et al. (2015) to analytically characterise the properties of tidal 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 semi-convective 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 θ₀ relative to the rotation axis, the rotation vector in the local coordinate system is $$2{\Omega }_0=(0,2{\Omega }_0\sin {\theta }_0,2{\Omega }_0\cos {\theta }_0)=2{\Omega }_0(0,\tilde{f},f),$$(1)

where Ω₀ is the rotation frequency of the star at the pole and $\tilde{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 $\sqrt{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.

2.2 Mean flow profile

In global spherical geometry, the mean flow based on a conical rotation profile Ω(θ) is written (e.g.in Guenel et al. 2016a) as: $$u=\Omega \times r=r\sin \theta \Omega (\theta )e_{\phi },$$(2)

where e_φ is the azimuthal unit vector and r and θ are the radius and co-latitude, respectively. We introduce u₀ = rsinθΩ₀e_φ, the mean flow at a point M inside the box (see Fig. 1), without differential rotation. We also use the shear contrast δΩ = Ω(θ) − Ω₀ (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 anti-solar-like rotating stars. Using the notations of Fig. 1, the centre C of the box is located at a distance r₀ sinθ₀ from the rotation axis. Accordingly, the latitudinal coordinate of the point M in the local frame is $$y=r\sin ({\theta }_0-\theta ).$$(3)

It should be noted that the radial coordinate r of the point M in spherical geometry can be written as r = r₀ + 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. 2004, 2018; Decressin et al. 2009; Alvan et al. 2013). Since y∕r and thus θ₀ − θ 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: $$\delta \Omega (\theta )={\Omega }_0\chi {\sin }^2\theta ,$$(4)

where χ is the magnitude of the shear between the equator and the pole. Performing a second-order Taylor expansion around a fixed co-latitude θ₀, such that θ = θ₀− y∕r₀, and at a specified depth r₀ inside the convective region, the local mean flow U can be recast as $$\begin{array}{ll}u-u_0\hfill & =r_0{\Omega }_0\chi {\sin }^3\theta \hfill \\ \hfill & \simeq {\Omega }_0\chi [r_0{\sin }^3{\theta }_0-3\cos {\theta }_0{\sin }^2{\theta }_{0}y\frac{3}{2}\hfill \\ \hfill & +\frac{y^2}{r_0}\left(3\cos {\theta }_0^2\sin {\theta }_0-\frac{3}{2}{\sin }^3{\theta }_0\right)+O(y^3)].\hfill \end{array}$$

We point out that the Taylor expansion must be pushed further at the pole θ₀ = 0 (and at the pole θ₀ = π with an opposite sign): $$u-u_0\simeq -\frac{{\Omega }_0\chi }{r_0^2}{y}^3+O(y^4).$$(6)

Accordingly, we can approximate a conical shear as a linear mean flow at the first order when the box is tilted. We recall that conical shear has been observed in the solar convective zone and is expected in slowly and moderately rotating 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).

2.3 System of equations

To derive the system of governing equations for tidal waves in the local reference frame, we made several hypotheses. Stratification terms, which usually drive the propagation of internal gravity waves, were kept for the sake of clarity and will be methodically kept or removed after applying the Boussinesq approximation and setting the equations for inertial waves. Moreover, we assume that the action of turbulence can be modelled as a Rayleigh friction term in the momentum equation with an effective frictional damping rate σ_f. This simplifies the analytical solution of the fluid equations compared to the usual modelling of turbulence as an effective viscous force (see in particular Ogilvie 2009). The momentum, continuity, and thermodynamic equations for tidal waves in a differentially rotating Cartesian framework are thus: $$\begin{array}{l}\frac{\text{d}u}{\text{d}t}+2{\Omega }_0\times u=-\frac{\nabla p}{\rho }+g-{\sigma }_{f}u+f,\\ \frac{\partial \rho }{\partial t}+\nabla \cdot (\rho u)=0,\\ \frac{\text{d}\rho }{\text{d}t}-\frac{1}{c_{\text{s}}^2}\frac{\text{d}p}{\text{d}t}=0,\end{array}$$

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 $\frac{\text{d}}{\text{d}t}=\frac{\partial }{\text{d}t}+u\cdot \nabla$ is the total derivative operator.

All variables are then linearised at the first order: Zero-order terms correspond to background equilibrium quantities, and first-order terms represent the leading perturbation. The local velocity, density, and pressure are therefore written as: $$\{\begin{array}{cc}u& =U(y)e_x+u^{\prime }\\ \rho & ={\rho }_0+ϵ{\rho }^{\prime }\\ p& =p_0+ϵp^{\prime }\end{array},$$(10)

where u = (u, v, w) in the local Cartesian basis. We have introduced the dimensionless parameter ϵ: $$ϵ=\frac{{(2{\Omega }_0)}^{2}L}{g},$$(11)

where we use 1∕(2Ω₀), 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Ω₀ L to scale velocity and ρ_TgL to scale pressure, with ρ_T the reference density. The dimensionless momentum equation of the mean flow is: $$ϵn\times U=-\frac{\nabla p_0}{{\rho }_0}-e_z,$$(12)

with n the unit vector parallel to the rotation axis. Projecting Eq. (12) into Cartesian coordinates, one can derive: $$\{\begin{array}{cc}{\partial }_x{p}_0& =0,\\ {\partial }_y{p}_0& =-ϵ{\rho }_{0}fU,\\ {\partial }_z{p}_0& =-{\rho }_0+ϵ{\rho }_0\tilde{f}U.\end{array}$$(13)

At the leading order in ϵ, one can recognise the hydrostatic balance and, at the first order, the geostrophic balance (the set is akin to the 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 $$N^2={ϵ}^{-1}\left(F\frac{{\partial }_z{p}_0}{{\rho }_0}-{\partial }_z\ln {\rho }_0\right),$$(14)

where we introduce the dimensionless number $F=gL/c_{\text{s}}^2$, which is small when filtering acoustic waves. Consequently, the curl of Eq. (12) gives $$\{\begin{array}{cc}& {\partial }_x{\rho }_0=0,\\ & {\partial }_y\ln {\rho }_0=ϵ\left(\tilde{f}{\partial }_{y}U-fUF\right),\\ & {\partial }_z\ln {\rho }_0=ϵ\left(F\tilde{f}U-N^2\right)-F,\end{array}$$(15)

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: $$\begin{array}{l}\frac{\text{d}u}{\text{d}t}+v{\partial }_{y}U+n\times u^{\prime }=-\frac{\nabla p^{\prime }}{{\rho }_0}-\frac{{\rho }^{\prime }}{{\rho }_0}{e}_z-{\sigma }_f{u}^{\prime }+f,\\ \nabla \cdot u^{\prime }=0,\\ \frac{\text{d}{\rho }^{\prime }}{\text{d}t}=-{\rho }_{0}v\tilde{f}{\partial }_{y}U,\end{array}$$

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.

2.4 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: $$2{\rho }_0{\Omega }_0\times U=-\nabla p_0.$$(19)

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: $$(2{\Omega }_0\cdot \nabla )({\rho }_{0}U)=0.$$(20)

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 $\tilde{f}=0$ at the poles, the latitudinal stratification term will not appear in the perturbed fluid equations (as we can see from Eq. (18)).

3 Dynamics of inertial waves at critical levels: analytical predictions

In this section, we analytically investigate the behaviour of inertial waves at critical levels in a non-dissipative fluid at various co-latitudes. For this purpose, we consider perturbations q in the normal mode $$q(x,y,z,t)=q(y)\exp i(k_{x}x+k_{z}z-\omega t)+\text{c}.\text{c}.,$$(21)

where ω is the complex inertial frequency, k_x and k_z are the real streamwise and vertical wavenumbers, respectively, and c.c. is the complex conjugate.

3.1 Wave propagation equation in the latitudinal direction

Using the modal form Eq. (21) for ρ, p, and u, we solved the set of hydrodynamic equations, Eqs. (16)–(18), for the latitudinal velocity v. Considering free inertial waves (i.e. without forcing terms), the set of perturbation equations can be recast into a single second-order ODE for v: $$A{v}^{″}+B{v}^{\prime }+Cv=0,$$(22)

where the prime now denotes the derivative according to y, and A, B, and C are the coefficients that can be simplified without friction as follows: $$\begin{array}{cc}A=& {\sigma }^2-{\tilde{f}}^2,\\ B=& -2k_x{\tilde{f}}^2\frac{U^{\prime }}{\sigma }-2i{k}_z\tilde{f}f,\\ C=& -k_{\perp }^2{\sigma }^2+k_z^{2}f(f-U^{\prime })-2i{k}_x{k}_z\tilde{f}f\frac{U^{\prime }}{\sigma }\\ & -2k_x^2{\tilde{f}}^2\frac{{U^{\prime }}^2}{{\sigma }^2}+\frac{k_x{U}^{″}}{\sigma }\left({\sigma }^2-{\tilde{f}}^2\right),\end{array}$$(23)

where $k_{\perp }=\sqrt{k_x^2+k_z^2}$ is the absolute wavenumber in the direction perpendicular to the y direction and σ = ω − k_xU 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 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 co-rotation 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 $\sigma =\pm \tilde{f}$ (we recall that $\tilde{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.

3.2 Propagation properties

3.2.1 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: $$\left({\sigma }^2-{\tilde{f}}^2\right){\partial }_{y,y}v-2\tilde{f}f{\partial }_{y,z}v+\left[{\sigma }^2+f\left(U^{\prime }-f\right)\right]{\partial }_{z,z}v=0,$$(24)

where we recover the Poincaré equation (for the propagation of inertial waves in the inviscid limit; Cartan 1922) in the meridional plane when there is no shear (U′ = 0) and at the poles ($\tilde{f}=0$ and f = 1). Moreover, we set v ∝ exp−i(k_zz − k_xx) so as to write the wave dispersion relation for the Doppler-shifted frequency σ: $${\sigma }^2=\frac{1}{\left|\right|k|{|}^2}\left[{\left(n\cdot k\right)}^2-k_z^{2}f{U}^{\prime }\right],$$(25)

where $\left|\right|k\left|\right|=k_y^2+k_z^2$ 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_z^{2}f{U}^{\prime }$) is present, which accounts for the latitudinal shear. Assuming that σ² takes positive values (as in Baruteau & Rieutord 2013; Guenel et al. 2016a), we therefore introduce $$\gamma =\sqrt{{(n\cdot k)}^2-k_z^{2}f{U}^{\prime }}.$$(26)

We can then specify the phase velocity in the meridional plane: $$v_{\varphi }=\frac{\sigma }{\left|\right|k|{|}^2}k=\pm \frac{\gamma }{\left|\right|k|{|}^3}k.$$(27)

In the same way, we can derive the expression for the group velocity in the meridional plane: $$v_{\text{g}}={\nabla }_k\sigma =\pm \frac{\gamma }{\left|\right|k|{|}^3}\left\{-k+{\left(\frac{\left|\right|k\left|\right|}{\gamma }\right)}^2\left[n(n\cdot k)-k_{z}f{U}^{\prime }{e}_z\right]\right\}.$$(28)

We note that without differential rotation, the group velocity reduces to its well-known expression for solid-body rotation (e.g. see Rieutord 2015): $$v_{\text{g}}=\pm \frac{k\wedge (n\wedge k)}{\left|\right|k|{|}^3}.$$(29)

Moreover, as in solid-body rotation, the group velocity (Eq. (28)) and the phase velocity (Eq. (27)) lie in perpendicular planes: v_g ⋅v_ϕ = 0.

When the box is located at the north pole (θ₀ = 0 in Fig. 1), by setting κ = 1 − U′, we recover $$\{\begin{array}{cc}{\sigma }^2& =\frac{k_z^2}{\left|\right|k|{|}^2}{\kappa }^2\\ v_{\varphi }& =\pm \kappa \frac{k_{z}k}{\left|\right|k|{|}^3}\\ v_{\text{g}}& =\pm \kappa \frac{k_y}{\left|\right|k|{|}^3}(-k_z{e}_y+k_y{e}_z),\end{array}$$(30)

as in Latter & Balbus (2009) and Baruteau & Rieutord (2013), where κ can be identified with the epicyclic frequency of Baruteau & Rieutord (2013) and k_y corresponds to the cylindrical component of the wavenumber (k_s).

3.2.2 Phase and group velocity at singularities

In this section, we derive the conditions required to meet singularities in terms of wavenumbers and shear, and we examine the implications for the phase and group velocities. When the box is inclined, for σ → 0 we must have γ → 0, meaning that v_ϕ →0 while  |v_g ⋅e_y|→∞ and |v_g ⋅e_z|→∞ (possibility 1).

Guenel et al. (2016a) found similar results by studying the propagation of free inertial waves in a global frame with conical shear, namely, when their parameter $\mathcal{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 $\sigma \to \pm \tilde{f}$, we either need |k_y|→∞ at fixed k_z, which implies v_ϕ →0 and v_g →0 and means that inertial waves cannot get through the critical level (possibility 2), or |k_z|→ 0 at fixed k_y, which gives v_ϕ ⋅e_z → 0 and v_g ⋅e_y → 0, while $|v_{\varphi }\cdot e_y|\to \tilde{f}/k_y$ and |v_g ⋅e_z|→ f∕k_y; the wave may then cross the critical level with some preferential direction (possibility 3).

Again, these conditions share some similarities with those observed for 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 = rsinθ and y ~ r₀(θ₀ − θ). 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 $\sigma =\pm \tilde{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 co-rotation σ = 0 or with the other critical levels at $\sigma =\pm \tilde{f}$. Hence, the local critical levels at $\sigma =\pm \tilde{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 co-rotation 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 v_ϕ →0 and v_g →0: The wave is totally absorbed at co-rotation. (ii) |k_y|→∞ at fixed k_z, which implies v_ϕ →0 and v_g →0: We can draw the same conclusion as in the previous case and this case is also analogous to the second possibility when the box is tilted. (iii) |k_z | → 0 at fixed k_y, which gives v_ϕ →0, v_g ⋅ e_y → 0 while |v_g ⋅e_z|→ κ∕k_y: the wave energy does not cross the 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 pole¹, we actually have a perfect match with the conditions given by Baruteau & Rieutord (2013) when using a cylindrical rotation profile for the mean flow.

3.2.3 Energetical aspects

In this section, we examine the energetic balance associated with inertial waves in our inclined shear box model, without assuming the 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_xc, where c is the longitudinal phase velocity (e.g. as in Booker & Bretherton 1967), we can use the first-order definition $$\begin{array}{ll}v\hfill & =\frac{\text{d}\zeta }{\text{d}t}=(U-c){\zeta }_x,\hfill \\ w\hfill & =\frac{\text{d}\eta }{\text{d}t}=(U-c){\eta }_x.\hfill \end{array}$$(31)

This first allows us to express the perturbed density from Eq. (18) as $\rho =-{\rho }_0\zeta \tilde{f}{U}^{\prime }$ (as a reminder, the symbol ′ has been dropped for perturbed quantities). Then, by multiplying the momentum equation, Eq. (16), by ρ₀ u, we get the energy balance equation: $$(U-c){\partial }_x{e}_{\text{k}}+\nabla \cdot (pu)={\rho }_0\left[-v{U}^{\prime }u+\zeta \tilde{f}{U}^{\prime }w-{\sigma }_f{u}^2+f\cdot u\right],$$(32)

where e_k = ρ₀u²∕2 is the kinetic energy density and pu is the so-called 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: $${\mathcal{P}}_{\text{ext}}={\mathcal{P}}_{\text{shear}}+D_{\text{visc}}+{\mathcal{P}}_{\text{f}},$$(33)

where we introduce, from left to right, the power of the external pressures at the boundaries − δ∕2 and δ∕2 on the perturbed latitudinal flow, the work of the shear, the viscous dissipation, and the forcing power, which read, respectively: $$\begin{array}{cc}& {\mathcal{P}}_{\text{ext}}={\left[\tilde{f}\overline{pv}\tilde{f}\right]}_{-\delta /2}^{\delta /2},\\ & {\mathcal{P}}_{\text{shear}}=-{\displaystyle \underset{-\delta /2}{\overset{\delta /2}{\int }}{\rho }_0}{U}^{\prime }\left(\overline{uv}-\tilde{f}\overline{w\zeta }\right)\text{d}y,\\ & D_{\text{visc}}=-{\displaystyle \underset{-\delta /2}{\overset{\delta /2}{\int }}{\rho }_0}{\sigma }_f\overline{u^2}\text{d}y,\\ & {\mathcal{P}}_{\text{f}}={\displaystyle \underset{-\delta /2}{\overset{\delta /2}{\int }}{\rho }_0}\overline{f\cdot u}\text{d}y,\end{array}$$(34)

where the bar represents the average in the (z, x) plane over one period. We note that the energy density and the acoustic flux in the x and z directions drop out in Eq. (33) when integrating because of the wave periodicity in those directions. The quantity ${\mathcal{P}}_{\text{shear}}$ can also be seen as the power transferred from the mean flow to the perturbation (or conversely) by the Reynolds stress: $$\tau =-{\rho }_0(\overline{uv}-\tilde{f}\overline{w\zeta })=-{\rho }_0\overline{v\left(u+\tilde{f}\eta \right)},$$(35)

where we use partial integration and the periodicity of perturbations in the x and z directions. At the pole, $\tilde{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. $\tau =-{\rho }_0\overline{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 $\overline{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: $$\frac{\text{d}}{\text{d}y}\overline{pv}=-{\rho }_0{U}^{\prime }\overline{v\left(u+\tilde{f}\eta \right)}.$$(36)

Using the same method as Broad (1995), we multiplied the x-projection of the inviscid force-free momentum equation by ζ: $$-p_x\zeta ={\rho }_0(U-c)\zeta u_x+{\rho }_0(U^{\prime }-f)v\zeta +{\rho }_0\tilde{f}w\zeta .$$(37)

By multiplying by (U − c), the latitudinal flux of energy can thus be written as: $$\overline{pv}=-{\rho }_0(U-c)\overline{v\left(u+\tilde{f}\eta \right)}.$$(38)

By differentiating this relationship with respect to y, and by equalising with Eq. (36), one can obtain: $$-{\rho }_0(U-c)\frac{\text{d}}{\text{d}y}\left\{\overline{v\left(u+\tilde{f}\eta \right)}\right\}=0,$$(39)

that is $(U-c)\frac{\text{d}\tau }{\text{d}y}=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.

3.2.4 Polarisation relations

For the forthcoming analysis, it is useful to derive expressions of the perturbed projected velocities and the perturbed reduced pressure² Π = p∕ρ₀, 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: $$\begin{array}{cc}\Pi & =\frac{1}{\sigma k_{\perp }^2}\left\{iA{v}^{\prime }+\left[i{k}_x\frac{U^{\prime }}{\sigma }A+f(k_z\tilde{f}-i{k}_x\sigma )\right]v\right\},\\ w& =\frac{1}{\sigma k_{\perp }^2}\left\{\left(i\sigma k_z-k_x\tilde{f}\right)v^{\prime }+k_x\left[i{k}_z(U^{\prime }-f)-k_x\tilde{f}\frac{U^{\prime }}{\sigma }\right]v\right\},\\ u& =\frac{1}{\sigma k_{\perp }^2}\left\{\left(i\sigma k_x+k_z\tilde{f}\right)v^{\prime }+k_z\left[i{k}_z(f-U^{\prime })+k_x\tilde{f}\frac{U^{\prime }}{\sigma }\right]v\right\}.\end{array}$$(40)

Without shear and at θ₀ = 0, we recover the polarisation relations in the solid-body rotation case (see e.g. Rieutord 2015).

3.2.5 Conservation of the wave action flux

While the latitudinal flux of energy is not conserved in the whole domain, there is a conserved quantity, called the wave action flux, as introduced in Grimshaw (1975a): $$\mathcal{A}=\frac{\overline{pv}}{\sigma }=\frac{{\rho }_0}{2}\Re \left\{\frac{\Pi v^{*}}{\sigma }\right\},$$(41)

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 $\mathcal{A}=-\tau /k_x$. By using the expression for the perturbed reduced pressure Π derived in the previous section, the wave action flux now reads: $$\mathcal{A}=\frac{1}{2}{\rho }_0\Re \left\{\frac{iA{v}_y{v}^{*}+f\tilde{f}{k}_z|v{|}^2}{{\sigma }^2{k}_{\perp }^2}\right\}.$$(42)

Unlike the latitudinal flux of energy, but similarly to the Reynolds stress, this wave action flux is conserved along the latitudinal direction. One can demonstrate that ${\mathcal{A}}^{\prime }=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 $\mathcal{A}$ is a measure of wave energy through a surface (in the (z, x) plane) since $\overline{pv}$ is the energy density transported by the group velocity³ V_g in the latitudinal direction (e.g. Bretherton & Garrett 1968; Mathis & de Brye 2012). It should be underlined that the wave action flux has been defined in the inviscid limit and is not conserved when the friction is taken into account.

3.3 Inertial waves at critical levels when the box is tilted

In this section, we analytically investigate waves passing through the various critical levels in the tilted box. We examine the behaviour of the waves around the co-rotation σ = 0 and the critical levels $\sigma =\pm \tilde{f}$ when the box is tilted (for the co-rotation when the box is at the pole, see Sect. 3.4).

3.3.1 Critical levels at $\sigma =\pm \tilde{f}$

In this subsection, we treat both singularities, $\sigma =\pm \tilde{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 $\sigma =\pm \tilde{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 $\sigma =\pm \tilde{f}$. For a linear mean flow profile U = Λy, with Λ a constant, y_± are given by: $$y_{\pm }=\frac{\omega \mp \tilde{f}}{k_x\Lambda }.$$(43)

Without any assumption on the mean flow profile, the first-order Taylor expansion of the ODE (Eq. (22)) near y_± is: $$v^{″}+\frac{\mathcal{B}}{y-y_{\pm }}{v}^{\prime }=0,$$(44)

with $$\mathcal{B}=1\pm i\frac{k_{z}f}{k_x{U^{\prime }}_{\pm }},$$(45)

where the symbol ± refers to the regular singularities⁴ y₊ and y₋ and ${U^{\prime }}_{\pm }$ is U′ evaluated at these singularities. The Frobenius method consists in injecting the power function ${(y-y_{\pm })}^{\lambda }$ into Eq. (44), with λ a constant to be determined (see e.g. Morse & Feshbach 1953). The corresponding indicial equation is then $$\lambda (\lambda -1)+\mathcal{B}\lambda =0,$$(46)

with solutions: $${\lambda }_{\pm }=\left\{0,\mp i\frac{k_{z}f}{k_x{U}_{{\pm }^{\prime }}}\right\}.$$(47)

Therefore, the two independent solutions of Eq. (44) can be written as follows: $$v_{1,\pm }(y)={\sum }_{n=0}^{+\infty }{a}_n{(y-y_{\pm })}^n\text{ and }{v}_{2,\pm }(y)={\sum }_{n=0}^{+\infty }{b}_n{(y-y_{\pm })}^{n+{\lambda }_{\pm }},$$(48)

where a_n and b_n are complex constants. Both solutions are valid in the vicinity of the critical level around which they are built in the complex plane, up to the next singularity if it exists. The coefficients a₀ and b₀ 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,±: $$v_{\pm }(y)\simeq a_0+b_0{(y-y_{\pm })}^{{\lambda }_{\pm }}.$$(49)

Owing to the existence of a branch point at y_± (since λ_± is complex),reconnecting solutions on either side of the critical levels is not straightforward. This requires both physical and mathematical arguments (see in particular Miles 1961; Booker & Bretherton 1967; Ringot 1998). In order to remove the degeneracy of the path from positive to negative y − y_± (i.e. choosing either e^+iπ or e^−iπ), we made use of a complex inertial frequency ω = ω_R + iω_I, assuming the radiation condition ω_I > 0. This condition ensures a non-growing wave towards infinity. The Taylor expansion of the base flow at the first order in y − y_± gives $$y-y_{\pm }=\frac{U(y)-U_{\pm }}{{U^{\prime }}_{\pm }},$$(50)

and, by definition, we have $$U_{\pm }=\frac{\omega \mp \tilde{f}}{k_x}.$$(51)

Consequently, the solution below the critical level is unambiguous in terms of the above solution coefficients and depends on $$\text{sign}[\Im \{y-y_{\pm }\}]=-\text{sign}\left[k_x{U^{\prime }}_{\pm }\right].$$(52)

In other words, when taking y − y_± to decrease from positive to negative values, its complex argument changes continuously from 0 to $-\text{sign}\left[k_x{U^{\prime }}_{\pm }\right]\pi$. Thus, the appropriate path for determining the branch of ${(y-y_{\pm })}^{{\lambda }_{\pm }}$ passes under(above) y_± as long as $k_x{U^{\prime }}_{\pm }>0$ ($k_x{U^{\prime }}_{\pm }<0$) (the same reasoning can be found in Grimshaw 1975a). Therefore, the solution on both sides of the critical level y_± is: $$\begin{array}{lll}\hfill & \hfill & v_{\pm }(y)\simeq \hfill \\ \hfill & \hfill & \{\begin{array}{ccc}& a_0+b_0|y-y_{\pm }{|}^{\mp i{k}_{z}f/k_x{U^{\prime }}_{\pm }}& \text{ for }y>y_{\pm }\\ & a_0+b_0|y-y_{\pm }{|}^{\mp i{k}_{z}f/k_x{U^{\prime }}_{\pm }}\text{exp}(\mp \pi k_{z}f/\left|k_x{U^{\prime }}_{\pm }\right|)& \text{ for }y<y_{\pm }\end{array}.\hfill \end{array}$$(53)

The remaining issue is now to know in which direction the wave is propagating. The second part of the solution can be assimilated to a wave-like solution with the varying latitudinal wavenumber $\mp (k_{z}f/k_x{U^{\prime }}_{\pm })\log |y-y_{\pm }|$. Moreover, according to Eq. (42), the wave action flux on either side of y_± is: $$\mathcal{A}=\frac{{\rho }_0{k}_{z}f}{2\tilde{f}{k}_{\perp }^2}\{\begin{array}{ccc}& |a_0{|}^2-|b_0{|}^2& \text{ for }y>y_{\pm }\\ & |a_0{|}^2-|b_0{|}^2\text{exp}(\mp 2\pi k_{z}f/|k_x{U^{\prime }}_{\pm }|)& \text{ for }y<y_{\pm }\end{array}\text{}\text{}\text{}.$$(54)

The group velocity gives the direction towards which the energy is transported (we recall that $V_{\text{g}}\mathcal{E}=\overline{pv}$, with V_g the group velocity and $\mathcal{E}$ the local energy density). By consequence, $\text{sign}(V_{\text{g}})=\text{sign}(\mathcal{A}\sigma )=\mp \text{sign}(f{k}_z)$ for the solution featuring the coefficient b₀. If k_zf is positive, this wave transports energy downwards (upwards) across the critical level $\sigma =\tilde{f}$ ($\sigma =-\tilde{f}$). If k_z f is negative, the wave transports energy upwards (downwards) across the critical level $\sigma =\tilde{f}$ ($\sigma =-\tilde{f}$). In all cases, the action flux of the wave with the amplitude |b₀| will be transmitted (in the direction given by the sign of k_zf and the critical level y₊ or y₋) by a factor $T_{{\theta }_0}$, where $$T_{{\theta }_0}=\exp \left(-2\pi \frac{\left|k_{z}f\right|}{\left|k_x{U^{\prime }}_{\pm }\right|}\right)=\exp \left(-2\pi \frac{\left|{\alpha }_k\cos {\theta }_0\right|}{\left|R_{\text{o}}\right|}\right),$$(55)

with α_k = k_z∕k_x and $R_{\text{o}}={U^{\prime }}_{\pm }$, after passing through the critical level. Such a wave will always be attenuated since $T_{{\theta }_0}\le 1$. The transmission factors $T_{{\theta }_0={10}^{°}}$ and $T_{{\theta }_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 θ₀ + 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 θ₀ = kπ with k ∈ {0, 1} are not well described by the attenuated factor $T_{{\theta }_0}$ and require a specific treatment, as discussed in Sect. 3.4.

It is important to note that, with fixed parameters {k_z, θ₀, R_o}, the attenuation of the wave action flux is specific to a single direction of wave propagation (i.e. the solution featuring the coefficient b₀). The solution of coefficienta₀ 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.

Fig. 2 Transmission rate $T_{{\theta }_0}$ of a wave passing through any of the critical levels defined by $\sigma =\pm \tilde{f}$ as a functionof the absolute value of the shear Rossby number |Ro| = |U′| (where U′ is scaled by 2Ω0) and the ratioof the horizontal wavenumbers αk = kz∕kx for a co-latitude of the box θ0 = 10° (left panel) and θ0 = 80° (right panel).

3.3.2 Inertial wave crossing co-rotation

We performed the same analysis as in the previous section to treat the co-rotation point y₀ where σ = 0 (i.e. U(y₀) = U₀ = ω∕k_x). The linearised ODE (Eq. (22)) near the co-rotation using the Taylor expansion of σ and U at the lowest orders is: $$\begin{array}{cc}& v^{″}+\left(-\frac{2}{y-y_0}+2i{k}_z\frac{f}{\tilde{f}}-\frac{U_{0^{″}}}{{U^{\prime }}_0}\right)v^{\prime }\\ & +\left[\frac{2}{{(y-y_0)}^2}-\frac{2i{k}_{z}f}{\tilde{f}(y-y_0)}+\frac{U_{0^{″}}}{U_{0^{\prime }}(y-y_0)}\right]v=0,\end{array}$$(56)

where $U_{0^{\prime }}$ and $U_{0^{\prime \prime }}$ are the first and second derivatives of the mean flow profile U evaluated at the critical level y₀. 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₂ of Eq. (56) that includes a logarithmic part, such as (e.g. Schmid et al. 2002) $$v_2(y)={\sum }_{n=0}^{+\infty }{b}_n{(y-y_0)}^{n+1}+L\ln (y-y_0)v_1(y),$$(57)

with $v_1(y)={\sum }_{n=0}^{+\infty }{a}_n{(y-y_0)}^{n+2}$ the first solution and a_n, b_n, and L complex coefficients. However, when injecting v₁ + v₂ into Eq. (56), one finds L = 0, meaning that a sole polynomial solution in the form $$v(y)={\sum }_{n=0}^{+\infty }{c}_n{(y-y_0)}^{n+1}$$(58)

includes all the solutions of Eq. (56), with c₀ = b₀ and a₁ = b₁ + a₀ determined by boundary conditions, and $c_{n, n\in {\mathbb{N}}^{*}\backslash \left\{1\right\}}=b_n+a_{n-1}$ determined by recurrence via the expansion of Eq. (22) around y₀. As a result, thewave action flux given by Eq. (42) becomes $$\mathcal{A}=\frac{{\rho }_0}{2k_x^2{U^{\prime }}_0^2{k}_{\perp }^2}\left[{\tilde{f}}^2\Im \left\{c_0^{*}{c}_1\right\}+f\tilde{f}{k}_z|c_0{|}^2\right]$$(59)

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 $-{\sigma }^2{k}_{\perp }^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.

3.4 Inertial waves when the box is at the poles

When the box is located at the north or south pole, $\tilde{f}=0$ and the ODE (Eq. (22)) is greatly simplified. For θ₀ = 0, the dimensionless wave propagation equation becomes $$v^{″}+\left[\frac{k_z^2(1-U^{\prime })}{{\sigma }^2}+\frac{k_x{U}^{″}}{\sigma }-k_{\perp }^2\right]v=0.$$(60)

At the southpole (i.e. θ₀ = π), 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\Omega =(0,\tilde{f},f)$ and there is a constant df∕dy = β, with f the Coriolis parameter (e.g. Miles 1961; Grimshaw 1975b; Gliatto & Held 2020). However, we cannot make a direct comparison at co-rotation, because the singularity in the equations for Rossby waves and inertial waves is not of thesame 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. Sucha rotation profile is expected in giant planets such as Jupiter and Saturn (Kaspi et al. 2017; Galanti et al. 2019, respectively) as a natural outcome of the 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.

3.4.1 Analytical solutions with constant shear

Analytical solutions of the ODE Eq. (60) are difficult to find for general profiles of the mean flow (e.g. a quadratic mean flow profile). A linear mean flow profile, on the other hand, has analytic solutions, which is why we use such a profile in this section (i.e. U = R_oy, with R_o the shear Rossby number that is taken as constant here). Equation (60) then becomes: $$v^{″}+\left[\frac{{\alpha }_k^2(1-R_{\text{o}})/R_{\text{o}}{}^2}{{(y-y_0)}^2}-k_{\perp }^2\right]v=0,$$(61)

where y₀ = ω∕(k_xR_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 ${\alpha }_k^2(1+R_{\text{o}})/R_{\text{o}}{}^2$ 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: $$v(y)=A{M}_{0,\mu }(\tilde{y})+B{M}_{0,-\mu }(\tilde{y}),$$(62)

where ỹ = 2k_⊥(y − y₀), $$\mu =\sqrt{\frac{1}{4}-\frac{{\alpha }_k^2(1-R_{\text{o}})}{R_{\text{o}}{}^2}},$$(63)

and A and B are complex constants given by boundary conditions. The Whittaker function M_0,−μ allows a quite straightforward analytic continuation: $$M_{0,\mu }({\text{e}}^{-i\pi }\tilde{y})=-i{\text{e}}^{-i\pi \mu }{M}_{0,\mu }(\tilde{y}).$$(64)

By consequence, the solution below the critical point y = y₀ is: $$v_{\text{W}}(y)=-i\varsigma A{\text{e}}^{-i\varsigma \pi \mu }{M}_{0,\mu }(-\tilde{y})-i\varsigma B{\text{e}}^{i\varsigma \pi \mu }{M}_{0,-\mu }(-\tilde{y}).$$(65)

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 $$R_{{\theta }_0=0,\pi }={\alpha }_k^2(1\mp R_{\text{o}})/R_{\text{o}}{}^2,$$(66)

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_{{\theta }_0=0,\pi }$ 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: $$v^{″}-k_{\perp }^{2}v=0,\text{ with }v\propto {\text{e}}^{\pm k_{\perp }y}.$$(67)

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.

3.4.2 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^{\prime }}_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 $$v^{″}+\frac{{\alpha }_k^2(1-R_{\text{o}})/R_{\text{o}}{}^2}{{(y-y_0)}^2}v=0,$$(68)

where $R_{\text{o}}={U^{\prime }}_0$ and $R_{\text{o}}{}^{\prime }={{U^{\prime }}^{\prime }}_0$. The indicial equation gives: $$\lambda =\frac{1}{2}\pm \mu .$$(69)

In the two next subsections, we examine the two cases where μ is imaginary or real.

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

3.4.3 Theoretical stable regime (R > 1/4)

We address here the case where R > 1∕4. The same analysis from Sect. 3.3 can be carried out to determine how a wave behaves upon crossing the co-rotation. The solutions of the indicial equation can be recast as $$\lambda =\frac{1}{2}\pm i\left|\mu \right|.$$(70)

The first-order solutions to Eq. (68) in the vicinity of y₀ are: $$\{\begin{array}{ccc}& a_0{(y-y_0)}^{\frac{1}{2}+i\left|\mu \right|}+b_0{(y-y_0)}^{\frac{1}{2}-i\left|\mu \right|},& y>y_0\\ & -i\varsigma \left[a_0{(y_0-y)}^{\frac{1}{2}+i\left|\mu \right|}{e}^{\varsigma \pi \left|\mu \right|}+b_0{(y_0-y)}^{\frac{1}{2}-i\left|\mu \right|}{e}^{-\varsigma \pi \left|\mu \right|}\right],& y<y_0\end{array}$$(71)

for $a_0,b_0\in ℂ$. 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 $\tilde{f}=0$, the wave action flux Eq. (42) reduces to $$\mathcal{A}=-\frac{{\rho }_0}{2k_{\perp }^2}\Im \left\{v_y{v}^{*}\right\},$$(72)

that is, injecting the solutions on both sides of the critical level: $$\mathcal{A}=-\frac{{\rho }_0}{2k_{\perp }^2}\left|\mu \right|\{\begin{array}{ccc}& |a_0{|}^2-|b_0{|}^2,& y>y_0\\ & -|a_0{|}^2{\text{e}}^{2\varsigma \pi \left|\mu \right|}+|b_0{|}^2{\text{e}}^{-2\varsigma \pi \left|\mu \right|},& y<y_0\end{array}.$$(73)

This formulation is quite similar to the expression of the Reynolds stress (τ) in vertically stratified mean flows, which can be found in Booker & Bretherton (1967) in Cartesian geometry. We recall that, indeed, $\tau =-k_x\mathcal{A}$. Moreover, given Eq. (35), the Reynolds stress in our model reads $\tau =-{\rho }_0\overline{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 co-rotation 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₀ | 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 $$T_{\text{F}}=\exp \{-2\pi |\mu \left|\right\}.$$(74)

This attenuation factor is shown versus R_o and α_k in Fig. 4. We observe that the wave is largely absorbed at the critical level and thus deposits most, if not all, its energy for most couples (α_k, R_o).

Fig. 4 Transmission rate TF of the wave action flux across the co-rotation when the box is located at the pole for R > 1∕4. It is displayed against the Rossby number Ro and the ratioof the horizontal wavenumbers αk = kz∕kx. The forbidden region where R < 1∕4 is shown in white, and the dark red cone corresponds to values of TF that are lower than 10−9.

3.4.4 Possible unstable regime (R < 1/4)

We now deal with the case where R < 1∕4 (i.e. μ is real). Contrary to the situation where R > 1∕4, we can no longer assimilate solutions to wave-like functions. The exponential form of solutions for R < 1∕4 near the critical level reads $$v(y)=\{\begin{array}{ccc}& a_0{(y-y_0)}^{\frac{1}{2}+\left|\mu \right|}+b_0{(y-y_0)}^{\frac{1}{2}-\left|\mu \right|},& y>y_0\\ & -i\varsigma [{\text{e}}^{-i\varsigma \pi \left|\mu \right|}{a}_0{(y-y_0)}^{\frac{1}{2}+\left|\mu \right|}& \\ & +{\text{e}}^{i\varsigma \pi \left|\mu \right|}{b}_0{(y-y_0)}^{\frac{1}{2}-\left|\mu \right|}],& y<y_0\end{array}$$(75)

and makes this region fully evanescent. Furthermore, the associated wave action flux is $$\mathcal{A}=-\frac{{\rho }_0}{k_{\perp }^2}\mu \{\begin{array}{ccc}& \Im \left\{a_0{b}_0^{*}\right\},& y>y_0\\ & -\Im \left\{a_0{b}_0^{*}{\text{e}}^{-2i\varsigma \pi \mu }\right\},& y<y_0\end{array}.$$(76)

Without knowing the direction of the wave or the energy flux, since $\text{sign}(V_{\text{g}})=\varsigma \text{sign}\left(\Im \left\{a_0{b}_0^{*}\right\}\right)$, 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 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.

4 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 first-order ODEs satisfied by v and Π, the combinationof which led to the wave propagation equation, Eq. (22). By imposing boundary conditions such that waves satisfy the dispersion relations (see also Appendix D.1), we examine the dynamics of inertial waves propagating in the shear region. Also, whenever possible, we analytically calculate the wave transmission and reflection coefficients as the wave-like solution crosses the shear region.

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 shearedregion II (grey-shaded) can have a linear, square, or cubic mean flow profile.

4.1 Description of the model

The mean flow profile that is used in the 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 $$U(y)=\{\text{}\begin{array}{cccccc}& 0 & \text{for}& y<0:& \text{zoneI},& \\ & \Lambda y^n & \text{for}& 0\le y\le 1:& \text{zoneII},& \\ & \Lambda & \text{for}& y>1:& \text{zoneIII},& \end{array}$$(77)

where U(y) is continuous at each interface and n is an integer: n = 1 for a linear shear flow, n = 2 for a square shearflow, or n = 3 for a cubic shear flow (see also Fig. 5). In zone I, we assume that there is an incident wave that enters the shear zone as well as a wave that is reflected at the interface between zones I and II or in zone II, that is: $$v(y)=A_I\exp (i{k}_{I}y)+A_R\exp (i{k}_{R}y),$$(78)

where A_I and A_R are the amplitudes and k_I and k_R are the wavenumbers of the incident and reflected waves, respectively. We further imposed as a boundary condition in zone III a transmitted wave that propagates towards positive y-values: $$v(y)=A_T\exp (i{k}_{T}y),$$(79)

where A_T and k_T are the amplitude and wavenumber of the transmitted wave, respectively. We imposed A_T = 1 without loss of generality and computed the remaining amplitudes A_I and A_R. More detailson the solutions and the dispersion relations of the waves in zones I and III can be found in Appendix D.1. We ensured that the transmitted wave carries energy upwards by deriving the wave action flux in zones I and III (see Appendix D.2).

We imposed the continuity of the latitudinal velocity and reduced pressure at the interfaces (at y = 0 and y = 1). By doing so, the wave action flux is continuous at both interfaces. Thus, in the absence of critical points, the wave action flux is conserved in the whole domain, namely ${\mathcal{A}}_{\text{T}}={\mathcal{A}}_{\text{I-R}}$, where ${\mathcal{A}}_{\text{T}}$ is the wave action flux of the transmitted wave, and ${\mathcal{A}}_{\text{I-R}}={\mathcal{A}}_{\text{I}}-{\mathcal{A}}_{\text{R}}$, where ${\mathcal{A}}_{\text{I}}$ and ${\mathcal{A}}_{\text{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 ${\sigma }^2={\tilde{f}}^2$, we added a small friction, σ_f = 10⁻⁸, 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 $$\{\text{}\text{}\text{}\text{}\text{}\begin{array}{cc}& F_{1,v}v+F_{1,p}\Pi =A{v}^{\prime },\\ & F_{2,v}v+F_{2,p}\Pi =A{\Pi }^{\prime },\end{array}$$(80)

where $$\begin{array}{cc}{F}_{1,v}& =\tilde{f}\frac{U^{\prime }}{\sigma }\left(k_x\tilde{f}+i{k}_{z}s\right)-(U^{\prime }-f)\left(k_{x}s+i{k}_z\tilde{f}\right),\\ F_{1,p}& =-i{k}_{\perp }^{2}s,\\ F_{2,v}& =is\left[s^2+f(U^{\prime }-f)\right]-i{\tilde{f}}^2\left(s+f\frac{U^{\prime }}{\sigma }\right),\\ F_{2,p}& =i{k}_{z}f\tilde{f}-k_{x}fs,\end{array}$$(81)