© The Authors 2026

1. Introduction

From the Earth’s atmosphere and oceans to the radiative core of solar-type stars and the radiative envelope of early-type stars, the combined action of the buoyancy force, of the Coriolis acceleration, and of heat and viscous diffusions plays a key role in wave-driven momentum transport and chemical mixing (e.g. Lindzen 1981; Vadas & Fritts 2005; Lott & Guez 2013; MacKinnon et al. 2017; Sutherland et al. 2019; Ribstein et al. 2022; Achatz et al. 2024; Charbonnel & Talon 2005; Talon & Charbonnel 2005; Rogers et al. 2013; Rogers 2015; Rogers & McElwaine 2017; Pinçon et al. 2017; Varghese et al. 2024). In this framework, gravito-inertial waves (GIWs; Dintrans & Rieutord 2000), for which buoyancy and the Coriolis acceleration are restoring forces, deposit their angular momentum at the place where they are damped (e.g. Zahn et al. 1997; Vadas & Fritts 2005), where they encounter critical layers (e.g. Booker & Bretherton 1967; Alvan et al. 2013; Lott et al. 2015), and where they break (e.g. Sutherland 2001; Staquet & Sommeria 2002; Rogers et al. 2013; Mathis 2025). If many works have focused on simple gravity waves (e.g. Press 1981; Schatzman 1993; Zahn et al. 1997), the Coriolis acceleration must be taken into account in many geophysical and astrophysical configurations such as for weakly stratified oceanic layers (e.g. Gerkema & Shrira 2005; Gerkema et al. 2008; Shrira & Townsend 2010) and for rapidly rotating young low-mass stars and early-type stars (e.g. Pantillon et al. 2007; Mathis et al. 2008). A key theoretical challenge must then be addressed when the buoyancy frequency (N, also called the Brunt-Väisälä frequency) and the inertia frequency (2Ω; Ω being the rotation) are on the same order of magnitude, that is, N ∼ 𝒞(2Ω) with 1 < 𝒞 < 10. In this regime, the propagation of GIWs must be modelled in a bidimensional inseparable mathematical framework (e.g. Dintrans & Rieutord 2000; Gerkema & Shrira 2005). This situation is challenging both theoretically and numerically.

In this context, most of the work studying transport and mixing by GIWs in a geophysical or astrophysical context have adopted the so-called traditional approximation of rotation (TAR), where the local projection of the rotation vector along the latitudinal direction is neglected (e.g. Eckart 1961; Bildsten et al. 1996; Lee & Saio 1997; Mathis et al. 2008; Mathis 2009; Lee et al. 2014). It involves neglecting the component of the Coriolis acceleration along the direction of the entropy and chemical stratifications, which allows us to uncouple the horizontal and the vertical dynamics and a separation of variables in the mathematical modelling of GIWs. This simplified mathematical formalism can be applied when 2Ω​ < ​​ < ​N, but it leads to erroneous predictions when 2Ω ∼ N for GIWs propagation (e.g. Friedlander 1987; Gerkema & Shrira 2005; Fruman 2009; Prat et al. 2016, 2018), hydrodynamical instabilities (e.g. Zeitlin 2018; Park et al. 2021; Toghraei & Billant 2022; Chew et al. 2023; Park & Mathis 2025) and excitation of GIWs (e.g. Mathis et al. 2014; Augustson et al. 2020).

When geophysical and astrophysical fluid dynamicists have to evaluate and parameterise the effect of GIWs on the large-scale atmospheric and oceanic circulation and its long-term evolution and on the secular structural, chemical, and dynamical evolution of stars, it is therefore mandatory to evaluate when a non-traditional (NT) modelling must be adopted and when the TAR can be used. This is the objective of this work, where we design a local Cartesian prototype model that we present in Sect. 2 to allow us to explore the linear damping of GIWs (in Sect. 3) and their breaking (in Sect. 4) for any N/(2Ω) ratio. In Sect. 5, we provide the corresponding NT parameterisation for interactions of GIWs and mean flows, and we discuss the perspectives and conclusions of this work in Sect. 6.

2. Gravito-inertial waves and transport

Our goal is to provide a global understanding of the effect of the complete Coriolis acceleration on the interactions of the waves and mean flows in stellar, planetary, and geophysical flows. We focus on the two channels that drive these interactions: the linear damping of waves, and their non-linear breaking.

2.1. A non-traditional prototype model

To reach this objective, we chose to adopt a local approach that allowed us to unravel these physical mechanisms step by step. To do this, we chose to consider a Cartesian box, centred on a point M with spherical coordinates (r,Θ,φ) inside a stably stratified rotating stellar or planetary region, where r is the radius, Θ is the angle between the local effective gravity g_eff¹ and the rotation vector Ω (see Fig. 1), and φ is the azimuthal angle. Mx, My, and Mz are the axes along the local longitudinal (azimuthal), latitudinal, and vertical (along g_eff) directions, respectively, with the corresponding Cartesian coordinates {x,y,z} and unit-vector basis $\{{\widehat{\mathbf{e}}}_x,{\widehat{\mathbf{e}}}_y,{\widehat{\mathbf{e}}}_z\}$. We introduce a reduced horizontal coordinate, χ, that makes an angle α with respect to the x−axis: χ = xcos α + y sin α. This so-called f-plane (Pedlosky 1982) corotates with the rotation $\mathbf{\Omega }=\mathrm{\Omega }{\widehat{\mathbf{e}}}_{\mathrm{\Omega }}$, where Ω is the global mean planetary/stellar rotation with ${\widehat{\mathbf{e}}}_{\mathrm{\Omega }}=sin\mathrm{\Theta }{\widehat{\mathbf{e}}}_y+cos\mathrm{\Theta }{\widehat{\mathbf{e}}}_z$.

Fig. 1. Local Cartesian box in the f-plane reference frame.

To model the mean longitudinally averaged zonal flow (differential rotation), we followed Mathis (2009) and introduced in an inertial frame,

$$\begin{array}{c}\hfill {\mathit{V}}_0=(V_{\mathrm{\Omega }}+U(y,z)){\mathbf{e}}_x,\end{array}$$(1)

where V_Ω = r sin θ Ω is the zonal flow associated with the global mean rotation Ω at the point M with spherical coordinates {r,θ}, and U is the local longitudinally averaged zonal flow, which depends on the local latitudinal (y) and vertical (z) coordinates in the general case. We considered as a first step a local horizontally averaged zonal flow $\overline{U}\left(z\right)$ that only depends on z and is the local equivalent of the so-called shellular differential rotation introduced by Zahn (1992) in the global spherical case.

To treat the Coriolis acceleration in the equations of GIW dynamics written in this reference frame, we introduced the two components of 2Ω

$$\begin{array}{c}\hfill f=2\mathrm{\Omega }cos\mathrm{\Theta }\text{and}\stackrel{\sim }{f}=2\mathrm{\Omega }sin\mathrm{\Theta }\end{array}$$(2)

along the vertical and latitudinal directions, respectively. In the most general case, both components must be taken into account to obtain a correct treatment of GIW dynamics (Gerkema & Shrira 2005; Gerkema et al. 2008; Mathis et al. 2014). Then, we assumed that a constant f leads to the so-called NT f-plane, in contrast to the usual traditional approximation, in which the horizontal component $\stackrel{\sim }{f}$ is neglected (Eckart 1961). The validity of this approximation is discussed throughout this work.

This local approach is valid only for mechanisms that can be modelled locally in a box, such that

$$\begin{array}{c}\hfill \{L_x,L_y,L_z\}<<R,\end{array}$$(3)

where L_x, L_y, and L_z are the lengths of the box in the x, y, and z directions, respectively, and R is the radius of the body. In this framework, the linear damping of GIWs and their potential non-linear breaking are maximum near the so-called critical layer where their local Doppler-shifted frequency $\widehat{\omega }=\omega +k_x\overline{U}$ (ω is the frequency of the studied GIW, k_x is its azimuthal wavenumber with its phase (ω t+k_xx), and $\overline{U}$ is the local horizontally averaged (vertically sheared) zonal flow introduced above) vanishes and where angular momentum is deposited or extracted (Alvan et al. 2013). In our prototype NT model, we considered this critical layer to be within our box and that the interactions of the wave and the mean zonal flow occur on spatial scales such that Eq. (3) is verified. In our model, we also assumed the Cowling approximation (Cowling 1941), in which the fluctuations of the self-gravity induced by GIWs are neglected. This approximation is verified for GIWs that propagate towards their critical layer. They oscillate rapidly along the vertical direction. When their local vertical wavelength is smaller than the local density and pressure scale height ($H_{\rho }\left(z\right)=|\mathrm{d}z/\mathrm{d}ln\overline{\rho }\left(z\right)|$ and $H_{\mathrm{p}}\left(z\right)=|\mathrm{d}z/\mathrm{d}ln\overline{P}\left(z\right)|$, where $\overline{\rho }$ and $\overline{P}$ are the hydrostatic density and pressure, respectively), we locally assumed the Boussinesq approximation, in which the density fluctuations are only taken into account in the buoyancy force. This theoretical framework is coherent with the framework that was adopted for stars in the study of angular momentum transport by low-frequency progressive internal gravity waves. The framework is based on the seminal work by Zahn et al. (1997). We adopted the so-called quasi-adiabatic approximation (see also Vadas & Fritts 2005, for the Earth’s atmosphere). Because of the low values of the viscosity and the heat diffusion in stars and in other celestial bodies, the damping of the waves they induce along their propagation is treated as a first-order local perturbation of their adiabatic propagation. Their adiabatic propagation is generally computed assuming the anelastic approximation, which allows us to take high density contrasts in stars and planetary atmospheres into account while filtering higher-frequencies acoustic waves (we refer to Lindzen 1981; Auclair-Desrotour et al. 2018, for planetary atmospheres), but their linear damping is computed locally using the Boussinesq approximation, as we propose here. This simplification can be made because the terms of the heat and viscous diffusions, which scale with the Laplacian applied to the components of the velocity of the waves and their entropy and temperature fluctuations and therefore as the square of the total wave vector, dominate the terms relative to the variations in the background hydrostatic quantities.

In Sect. 2.2 we therefore derive the dissipative Poincaré equation, which describes the non-adiabatic propagation of GIWs in a local Boussinesq model. Next, in Sect. 2.3, we recall the main key results for their adiabatic propagation, which we use to study their linear damping in Sect. 3 and their convective and shear-induced breaking in Sect. 4.

2.2. The dissipative Poincaré equation

To derive the Poincaré equation that governs GIW dynamics, we write the linearised equations of motion of the studied stably stratified rotating stellar or planetary fluid in our local Cartesian setup assuming the Boussinesq and the Cowling approximations, as discussed in the previous section. We introduce the GIW velocity field $\mathit{u}=(u,v,w)={\mathit{u}}_H+{\mathit{u}}_V$, where u, v and w are the components in the local azimuthal, latitudinal, and vertical directions, respectively; ${\mathit{u}}_H=(u,v,0)$ and ${\mathit{u}}_V=(0,0,w)$ are the horizontal and vertical velocities. Next, we define the fluid buoyancy within the Boussinesq approximation following Gerkema et al. (2008),

$$\begin{array}{c}\hfill b=-{\overline{g}}_{\mathrm{eff}}\left(z\right)\frac{\rho \prime (\mathit{r},t)}{{\overline{\rho }}_{\mathrm{B}}},\end{array}$$(4)

where ρ_B is a given uniform and constant reference Boussinesq density, and the corresponding Brunt-Väisälä frequency

$$\begin{array}{cc}\hfill N^2\left(z\right)& =-\frac{{\overline{g}}_{\mathrm{eff}}}{{\overline{\rho }}_{\mathrm{B}}}(\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}z}+\frac{\overline{\rho }{\overline{g}}_{\mathrm{eff}}}{c_{\mathrm{s}}^2})\hfill \\ \hfill & =-{\overline{g}}_{\mathrm{eff}}(\frac{1}{{\overline{\rho }}_{\mathrm{B}}}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}z}-\frac{\overline{\rho }}{{\overline{\rho }}_{\mathrm{B}}}\frac{1}{{\mathrm{\Gamma }}_1\overline{P}}\frac{\mathrm{d}\overline{P}}{\mathrm{d}z}),\hfill \end{array}$$(5)

where ρ′ is the density fluctuation, $\overline{\rho }$ and $\overline{P}$ are the reference hydrostatic background density and pressure, respectively, ${\mathrm{\Gamma }}_1={(\partial ln\overline{P}/\partial ln\overline{\rho })}_{\overline{S}}$ is the first adiabatic exponent, where $\overline{S}$ is the macroscopic entropy, and $c_{\mathrm{s}}=\sqrt{({\mathrm{\Gamma }}_1\overline{P})/\overline{\rho }}$ is the sound speed. In the main text, we consider that buoyancy only arises from the entropy stratification. This is not the case in the oceans and in the interiors of stars and of giant planets, however, where the chemical stratification (the salinity in the case of the oceans) plays an important role (e.g. Thorpe 2005; Maeder 2009; Leconte & Chabrier 2012). This more general case is treated in Appendix C. The three linearised components of the momentum equation are given by

$$\begin{array}{c}\hfill \{\begin{array}{c}{\displaystyle {\partial }_{t}w-\stackrel{\sim }{f}u=-\frac{1}{{\overline{\rho }}_{\mathrm{B}}}{\partial }_{z}p\prime +b+\nu {\mathrm{\nabla }}^{2}w}\hfill \\ {\displaystyle {\partial }_{t}v+fu=-\frac{1}{{\overline{\rho }}_{\mathrm{B}}}{\partial }_{y}p\prime +\nu {\mathrm{\nabla }}^{2}v}\hfill \\ {\displaystyle {\partial }_{t}u-fv+\stackrel{\sim }{f}w=-\frac{1}{{\overline{\rho }}_{\mathrm{B}}}{\partial }_{x}p\prime +\nu {\mathrm{\nabla }}^{2}u}\hfill \end{array},\end{array}$$(6)

where t is the time, p′ the pressure fluctuation, and ν the viscosity. Next, we write the continuity equation in the Boussinesq approximation,

$$\begin{array}{c}\hfill {\partial }_{z}w+{\partial }_{y}v+{\partial }_{x}u=0,\end{array}$$(7)

and the equation for energy conservation,

$$\begin{array}{c}\hfill {\partial }_{t}b+N^2\left(z\right)w=\kappa {\mathrm{\nabla }}^{2}w,\end{array}$$(8)

where κ is the thermal diffusivity. Eliminating the horizontal components of the velocity, the pressure, and the buoyancy, we reduce the system to the dissipative Poincaré equation for the vertical velocity,

$$\begin{array}{c}\hfill D_{\kappa }{D}_{\nu }^2{\mathbf{\nabla }}^{2}w+D_{\kappa }{(2\mathbf{\Omega }·\mathbf{\nabla })}^{2}w+D_{\nu }\left[N^2{\mathrm{\nabla }}_{\perp }^{2}w\right]=0,\end{array}$$(9)

where $D_{\nu }=({\partial }_t-\nu {\mathbf{\nabla }}^2)$, $D_{\kappa }=({\partial }_t-\kappa {\mathbf{\nabla }}^2)$ and ${\mathrm{\nabla }}_{\perp }^2$ is the horizontal Laplacian.

2.3. Propagation of adiabatic gravito-inertial waves

To quantify the transport of momentum by GIWs in our prototype model for any value of N/2Ω, we follow in Sects. 3, 4, and 5 the methods by Vadas & Fritts (2005), Zahn et al. (1997), Lott & Guez (2013), and Mathis (2025). On the one hand, this means that for the linear damping of GIWs through the viscous and heat diffusion, we use the so-called quasi-adiabatic method. In this framework, dissipative processes are therefore treated as a first-order perturbation when compared to the adiabatic propagation of waves. On the other hand, we consider the adiabatic vertical wave number to compute the saturation velocity above which a linear GIW breaks (Lott & Guez 2013; Mathis 2025). As a first step, we thus considered the Poincaré equation (Eq. (9)) in the limit where {ν,κ} → 0. We expanded the vertical velocity of a monochromatic wave of frequency ω as

$$\begin{array}{c}\hfill w=W(\chi ,z)exp[\text{i}\omega t],\end{array}$$(10)

where we recall the definition of the reduced horizontal coordinate χ = xcos α + y sin α. We obtained the adiabatic Poincaré equation for GIWs,

$$\begin{array}{c}\hfill [N^2(z)-{\omega }^2+{\stackrel{\sim }{f}}_{\text{s}}^2]{\partial }_{\chi ,\chi }W+2f{\stackrel{\sim }{f}}_{\text{s}}{\partial }_{z,\chi }W+[f^2-{\omega }^2]{\partial }_{z,z}W=0,\end{array}$$(11)

where ${\stackrel{\sim }{f}}_{\text{s}}=\stackrel{\sim }{f}sin\alpha$. Following Gerkema & Shrira (2005), we introduced the expansion for ω ≠ f:

$$\begin{array}{c}\hfill W=\widehat{W}(z)exp\left[\text{i}{k}_{\perp }(\chi +\stackrel{\sim }{\delta }z)\right],\end{array}$$(12)

where

$$\begin{array}{c}\hfill \stackrel{\sim }{\delta }=\frac{f{\stackrel{\sim }{f}}_{\text{s}}}{{\omega }^2-f^2}·\end{array}$$(13)

Substitution of Eq. (12) into Eq. (11) leads to the equation of a simple harmonic oscillator along the vertical direction,

$$\begin{array}{c}\hfill \frac{{\text{d}}^2\widehat{W}}{\text{d}{z}^2}+k_V^2\left(z\right)\widehat{W}=0,\end{array}$$(14)

where

$$\begin{array}{c}\hfill k_V^2=k_{\perp }^2[\frac{N^2\left(z\right)-{\omega }^2}{{\omega }^2-f^2}+{\left(\frac{\omega {\stackrel{\sim }{f}}_{\text{s}}}{{\omega }^2-f^2}\right)}^2].\end{array}$$(15)

When the GIW vertical wavelength (λ_V = 2π/k_V) is very small compared to the density scale height (H_ρ), that is, k_V H_ρ​ > ​​ > ​1, the Jeffreys-Wentzel-Kramers-Brillouin (JWKB; we refer the reader to Fröman & Fröman 2005, for a detailed derivation) solutions can be used, where $\widehat{W}\left(z\right)\equiv A/k_V^{1/2}exp[i\epsilon {\displaystyle {\int }_{z_0}^z}{k}_V(z\prime )\mathrm{d}z\prime ]$. A is a prescribed amplitude depending on the wave excitation mechanism. For solar-type stars (where the excitation region is above the propagation zone), we have ε = −1; the phase propagates outwards while the energy propagates inwards in the super-inertial regime. For early-type stars (where the excitation region is below the propagation zone), we have ε = 1; the phase propagates inwards while the energy propagates outwards (Mathis 2025). For a uniform local value of the Brunt-Väisälä frequency (N), JWKB solutions simplify to plane-wave solutions $\widehat{W}(z)=A\prime exp\left(\text{i}\epsilon k_{V}z\right)$, where A′ is the corresponding amplitude. When k_V² > 0, we obtain propagative GIWs along the vertical direction, but evanescent GIWs when k_V² < 0. The vertical component of the velocity becomes in the propagative case

$$\begin{array}{c}\hfill w=A\prime exp[\text{i}({\stackrel{\sim }{k}}_{V}z+k_{\perp }\chi +\omega t)],\end{array}$$(16)

where we introduced the total vertical wave number

$$\begin{array}{c}\hfill {\stackrel{\sim }{k}}_V=\epsilon k_V+k_{\perp }\stackrel{\sim }{\delta }.\end{array}$$(17)

Its NT component $k_{\perp }\stackrel{\sim }{\delta }$ corresponds to the 2D behaviour of the propagation of GIWs when the complete Coriolis acceleration is taken into account (we refer the reader to Dintrans & Rieutord 2000; Mathis et al. 2014, for a detailed discussion). Their vertical propagation is coupled to their propagation along the horizontal direction in this case because ${\stackrel{\sim }{k}}_z$ depends on k_⊥. In the polar traditional case, this coupling vanishes because δ → 0, and the problem becomes separable.

From a local point of view, at a given position (z, θ), the propagation condition k_V² > 0, allows us to identify the possible frequency range for propagative GIWs (Gerkema & Shrira 2005; Mathis et al. 2014),

$$\begin{array}{c}\hfill {\omega }_{-}<\omega <{\omega }_{+},\end{array}$$(18)

where

$$\begin{array}{c}\hfill {\omega }_{\pm }=\frac{1}{\sqrt{2}}\sqrt{[N^2+4{\stackrel{\sim }{\mathrm{\Omega }}}^2]\pm \sqrt{{[N^2+4{\stackrel{\sim }{\mathrm{\Omega }}}^2]}^2-{(2fN)}^2}},\end{array}$$(19)

where we defined a modified rotation rate of the planet,

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathrm{\Omega }}\equiv \frac{1}{2}\sqrt{f^2+{\stackrel{\sim }{f}}_{\text{s}}^2}=\mathrm{\Omega }\sqrt{1-{sin}^2\mathrm{\Theta }{cos}^2\alpha }.\end{array}$$(20)

From a global point of view, super-inertial waves (with 2Ω < ω) propagate at all latitudes, while they are trapped within an equatorial belt with θ ∈ [θ_c,π−θ_c] with

$$\begin{array}{c}\hfill {\theta }_c=\mathrm{Arccos}\left[\frac{\omega }{2\mathrm{\Omega }}\sqrt{1+{\left(\frac{2\mathrm{\Omega }}{N}\right)}^2[1-{\left(\frac{\omega }{2\mathrm{\Omega }}\right)}^2]}\right]\end{array}$$(21)

in the sub-inertial regime (with ω < 2Ω; see, e.g. Prat et al. 2016). We note that in the traditional limit, in which 2Ω​ < ​​ < ​N, we recover the classical simplest result θ_c = Arccos[ω/(2Ω)] (e.g. Lee & Saio 1997).

2.4. Energy and momentum transport

The principal objective of our work is to examine the effect of taking the full Coriolis acceleration into account or assuming the TAR in modelling the interactions of GIWs and the mean zonal flows. Therefore, we introduced the flux of energy transported along the vertical direction by GIWs and the related flux of (angular) momentum. In this section, we thus focus on the energetics of GIW propagation.

In our setup, following André et al. (2017), the energy density flux in the vertical direction is ${\langle pw\rangle }_H$ in real variables, so that we obtain

$$\begin{array}{c}\hfill F_{\mathrm{E};V}\left(z\right)={\langle \mathrm{Re}(w)\mathrm{Re}(p)\rangle }_H=\frac{1}{2}\left(\widehat{W}{\widehat{P}}^{\ast }\right),\end{array}$$(22)

where ^* is the complex conjugate, Re the real part, and ${\langle ···\rangle }_H$ the horizontal average. Using the polarisation relation (A.9), in which $\widehat{P}$ is expressed as a function of $\widehat{W}$, and considering low-frequency GIWs for which we can assume the JWKB approximation, we obtain

$$\begin{array}{c}\hfill F_{\mathrm{E};V}\approx -\frac{1}{2}\epsilon \overline{\rho }\left(\frac{f^2-{\omega }^2}{\omega k_{\perp }^2}\right)k_V{|\widehat{W}|}^2.\end{array}$$(23)

In the simplest case of plane waves given in Eq. (16), this simplifies into

$$\begin{array}{c}\hfill F_{\mathrm{E};V}=-\frac{1}{2}\epsilon \overline{\rho }\left(\frac{f^2-{\omega }^2}{\omega k_{\perp }^2}\right)k_V{|A\prime |}^2,\end{array}$$(24)

where we recall that A′ is the amplitude of the vertical component of the velocity in this case. For a solar-type (early-type) star, we recover that in the super-inertial regime F_E; V < 0 (F_E; V > 0).

We recall that we introduced the mean background zonal flow in an inertial frame in Eq. (1),

$$\begin{array}{c}\hfill {\mathit{V}}_0=(V_{\mathrm{\Omega }}+\overline{U}\left(z\right)){\mathbf{e}}_x,\end{array}$$(25)

where V_Ω = r sin θ Ω is the zonal flow associated with the global mean rotation Ω at the point M with spherical coordinates {r,θ}, and $\overline{U}$ is the local horizontally averaged zonal flow. In the momentum equation projected on the local Cartesian basis (Eq. (6)), this vertically sheared zonal flow transforms the temporal variation ∂_t into ${\partial }_t+\overline{U}{\partial }_x$, where the second term is the advection term due to the mean flow, while the Coriolis acceleration along the longitudinal direction becomes $-fv+\stackrel{\sim }{f}w+w{\partial }_z\overline{U}$. When we assume as a first step that the ratio $R={\partial }_z\overline{U}/\left(2\mathrm{\Omega }\right)$ is low, the sole modification of the momentum equation is related to the advection term due to the mean flow that leads to the transformation of ω into the Doppler-shifted frequency $\widehat{\omega }=\omega +k_x\overline{U}$. The assumption that R is low is close to the weak differential rotation assumption adopted in the stellar case by Mathis et al. (2008) to unravel the effects of the global rotation of a stellar radiation zone onto the angular momentum transport by GIWs in a traditional framework.

In this approximation, the evolution equation for the mean zonal flow is given by

$$\begin{array}{c}\hfill \overline{\rho }\frac{\mathrm{d}\overline{U}}{\mathrm{d}t}={\partial }_z{F}_{\mathrm{AM};V}\left(z\right),\end{array}$$(26)

with the momentum flux transported by low-frequency GIWs along the vertical direction,

$$\begin{array}{cc}\hfill F_{\mathrm{AM};V}& =-\frac{k_x}{\widehat{\omega }}{F}_{\mathrm{E};V}=\frac{1}{2}\epsilon \overline{\rho }\left(\frac{f^2-{\widehat{\omega }}^2}{{\widehat{\omega }}^2}\right)cos\alpha \frac{k_V}{k_{\perp }}{|\widehat{W}|}^2\hfill \\ \hfill & =\frac{1}{2}\epsilon \overline{\rho }\left(\frac{f^2-{\widehat{\omega }}^2}{{\widehat{\omega }}^2}\right)cos\alpha \frac{k_V}{k_{\perp }}{|A\prime |}^2,\hfill \end{array}$$(27)

where we used the relation k_x = k_⊥cos α. The prograde (retrograde) waves are such that k_x < 0 (k_x > 0), and thus, α ∈ [π/2,3π/2] (α ∈ [0,π/2]U [3π/2, 2π]). For super-inertial GIWs that propagate in solar-type stars for which ε = −1, we recover that prograde (retrograde) waves deposit (extract) momentum in the radiative core (e.g. Mathis et al. 2008; Mathis 2009).

3. Linear damping

3.1. Linear spatial damping rate

The first channel that allows GIWs to exchange momentum with mean zonal flows is their linear damping by viscous and heat diffusions. We focus on low-frequency rapidly oscillating GIWs in space, which are most strongly damped (e.g. Zahn et al. 1997; Mathis 2009; Alvan et al. 2015), and which thus transport momentum most efficiently. For these waves, we applied the JWKB approximation.

In this framework, we formally write the vertical velocity as a quasi-plane wave,

$$\begin{array}{cc}\hfill w& =E\left(\mathit{r}\right)exp\left[i(\mathrm{\Phi }\left(\mathit{r}\right)+\omega t)\right]\hfill \\ \hfill & \equiv E\left(\mathit{r}\right)exp\left[i(\epsilon {\displaystyle {\int }_{z_0}^z}{\stackrel{\sim }{k}}_V\mathrm{d}z\prime +k_{\perp }\chi +\omega t)\right],\hfill \end{array}$$(28)

where E is an envelope function that varies slower than the phase Φ, and ${\stackrel{\sim }{k}}_V$ and k_⊥ are the total vertical and horizontal wave numbers, respectively.

We followed the path of Zahn et al. (1997), who computed the linear damping considering the local dispersion relation derived from the dissipative wave equation within the Boussinesq approximation. Substituting Eq. (28) in Eq. (9), we obtained the dispersion relation:

$$\begin{array}{cc}\hfill (\omega -i\kappa k^2){(\omega -i\nu k^2)}^2{k}^2& =(\omega -i\kappa k^2)({\stackrel{\sim }{f}}_s^2{k}_{\perp }^2+2f{\stackrel{\sim }{f}}_s{k}_{\perp }{\stackrel{\sim }{k}}_V+f^2{\stackrel{\sim }{k}}_V^2)\hfill \\ \hfill & +(\omega -i\nu k^2)N^2{k}_{\perp }^2.\hfill \end{array}$$(29)

We assumed the quasi-adiabatic approximation, in which we considered the viscous friction and the heat diffusion as first-order perturbations of the adiabatic propagation (see Press 1981; Zahn et al. 1997; Vadas & Fritts 2005, for stellar interiors and the Earth’s atmosphere, respectively). Equation (29) reads

$$\begin{array}{cc}& -i{k}^2{\omega }^{-1}[(\kappa +2\nu ){\omega }^2{k}^2-\kappa (f^2{\stackrel{\sim }{k}}_V^2+2f{\stackrel{\sim }{f}}_s{\stackrel{\sim }{k}}_V{k}_{\perp }+{\stackrel{\sim }{f}}_s^2{k}_{\perp }^2)-\nu N^2{k}_{\perp }^2]\hfill \\ \hfill & =A{k}_{\perp }^2+2B{k}_{\perp }{\stackrel{\sim }{k}}_V+C{\stackrel{\sim }{k}}_V^2,\hfill \end{array}$$(30)

where $A=N^2-{\omega }^2+{\stackrel{\sim }{f}}_s^2$, $B=f{\stackrel{\sim }{f}}_s$, and C = f² − ω².

To derive the spatial damping, we assumed that the frequency ω is real and expanded the vertical wave number as

$$\begin{array}{c}\hfill {\stackrel{\sim }{k}}_V\equiv {\stackrel{\sim }{k}}_{V_0}\left(z\right)+{\stackrel{\sim }{k}}_{V_1}\left(z\right)={\stackrel{\sim }{k}}_{V_0}\left(z\right)+i{\widehat{\tau }}_j\left(z\right),\end{array}$$(31)

with ${\stackrel{\sim }{k}}_{V_0}$ the adiabatic vertical wave number and ${\widehat{\tau }}_j$, where j ≡ {NT, T,Ω = 0} indicates the adopted approximation, the first-order local linear damping rate. Equation (30) is then expanded as

$$\begin{array}{cc}& -i({\stackrel{\sim }{k}}_{V_0}^2+k_{\perp }^2){\omega }^{-1}[(\kappa +2\nu ){\omega }^2({\stackrel{\sim }{k}}_{V_0}^2+k_{\perp }^2)\hfill \\ \hfill & -\kappa (f^2{\stackrel{\sim }{k}}_{V_0}^2+2f{\stackrel{\sim }{f}}_s{\stackrel{\sim }{k}}_{V_0}{k}_{\perp }+{\stackrel{\sim }{f}}_s^2{k}_{\perp }^2)-\nu N^2{k}_{\perp }^2]\hfill \\ \hfill & =A{k}_{\perp }^2+2B{k}_{\perp }{\stackrel{\sim }{k}}_{V_0}+i2B{k}_{\perp }{\widehat{\tau }}_{\mathrm{NT}}+C{\stackrel{\sim }{k}}_{V_0}^2+i2C{\stackrel{\sim }{k}}_{V_0}{\widehat{\tau }}_{\mathrm{NT}}.\hfill \end{array}$$(32)

Keeping the zeroth-order terms, we recovered the adiabatic dispersion relation

$$\begin{array}{c}\hfill A{k}_{\perp }^2+2B{k}_{\perp }{\stackrel{\sim }{k}}_{V_0}+C{\stackrel{\sim }{k}}_{V_0}^2=0,\end{array}$$(33)

while first-order terms allowed us to derive the vertical damping rate,

$$\begin{array}{cc}\hfill {\widehat{\tau }}_{\mathrm{NT}}& =-\frac{1}{2\omega }\frac{{\stackrel{\sim }{k}}_{V_0}^2+k_{\perp }^2}{B{k}_{\perp }+C{\stackrel{\sim }{k}}_{V_0}}\hfill \\ \hfill & \times (\kappa [{\omega }^2({\stackrel{\sim }{k}}_{V_0}^2+k_{\perp }^2)-(f^2{\stackrel{\sim }{k}}_{V_0}^2+2f{\stackrel{\sim }{f}}_s{\stackrel{\sim }{k}}_{V_0}{k}_{\perp }+{\stackrel{\sim }{f}}_s^2{k}_{\perp }^2)]\hfill \\ \hfill & +\nu [2{\omega }^2({\stackrel{\sim }{k}}_{V_0}^2+k_H^2)-N^2{k}_{\perp }^2]).\hfill \end{array}$$(34)

When considering a vertically sheared flow, we substitute $\omega \equiv \widehat{\omega }$ as soon as the ratio $R={\partial }_z\overline{U}/\left(2\mathrm{\Omega }\right)$ remains low.

The global damping rate as defined by Press (1981), Schatzman (1993), and Zahn et al. (1997) is then given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{G}}\left(z\right)=\epsilon {\displaystyle {\int }_{z_0}^z}|{\widehat{\tau }}_j|\mathrm{d}z\prime .\end{array}$$(35)

3.2. An example: The stellar case

To demonstrate and illustrate the mandatory use of NT formalisms in the general case and constrain the domain of validity of the TAR as in Auclair-Desrotour et al. (2018), we considered the case of stellar radiation zones where the Prandtl number P_r = ν/K​ < ​​ < ​1 (e.g. P_r ≈ 10⁻⁶ for the radiative core of the Sun). We therefore focused on the radiative damping as the sole cause of the linear damping of GIWs (e.g. Press 1981; Schatzman 1993; Zahn et al. 1997; Mathis et al. 2008; Mathis 2009). In this limit, the linear NT radiative damping is given by

$$\begin{array}{c}\hfill {\widehat{\tau }}_{\mathrm{NT}}=-\frac{\kappa }{2}\frac{1}{\omega }{k}_{\perp }^3\frac{[{(\delta +\sqrt{\mathcal{D}})}^2+1][{\omega }^2[{(\delta +\sqrt{\mathcal{D}})}^2+1]-\mathcal{F}]}{B+C(\delta +\sqrt{\mathcal{D}})},\end{array}$$(36)

where

$$\begin{array}{c}\hfill \mathcal{D}=\frac{B^2-AC}{C^2}\end{array}$$(37)

and

$$\begin{array}{c}\hfill \mathcal{F}=f^2{(\delta +\sqrt{\mathcal{D}})}^2+2f{\stackrel{\sim }{f}}_s(\delta +\sqrt{\mathcal{D}})+{\stackrel{\sim }{f}}_s^2.\end{array}$$(38)

In the traditional case for which $\stackrel{\sim }{f}={\stackrel{\sim }{f}}_s=B=\delta =0$, this reduces to

$$\begin{array}{c}\hfill {\widehat{\tau }}_{\mathrm{T}}=\frac{\kappa }{2}{k}_{\perp }^3\frac{N^2(N^2-f^2)}{\omega {({\omega }^2-f^2)}^{3/2}{(N^2-{\omega }^2)}^{1/2}}.\end{array}$$(39)

In the low-frequency (ω​ < ​​ < ​N) non-rotating case (f = 0) case, we recover the classical result ${\widehat{\tau }}_{\mathrm{\Omega }=0}=1/2\times \kappa k_{\perp }^3\times N^3/{\omega }^4$ (e.g. Press 1981; Schatzman 1993; Zahn et al. 1997).

In Fig. 2 we study the variation in the ratios of the vertical linear damping rate and of the adiabatic vertical wave number computed assuming the TAR by those computed in the NT case as a function of the wave frequency normalised by the inertial frequency (ω/2Ω) in a weakly stratified case for which S = N/2Ω = 5 (left) and in a strongly stratified case for which S = N/2Ω = 50 (right) for different colatitudes θ = {0,π/8,π/4,3π/8,π/2}. We chose the specific case α = π/2 (i.e. the axisymmetric GIWs for which k_x = 0) to maximise the NT effects since ${\stackrel{\sim }{f}}_s=\stackrel{\sim }{f}sin\alpha$. Except at the pole in the traditional case) and at the equator, these ratios decrease when the normalised frequency ω/2Ω decreases. This decrease strengthens when the sub-inertial regime is entered (i.e. ω < 2Ω) and tends to ∞ with $\{{k_V^2}_{\mathrm{T}}/{k_V^2}_{\mathrm{NT}},{\widehat{\tau }}_{\mathrm{T}}/{\widehat{\tau }}_{\mathrm{NT}}\}\to 0$ when ω → ω₋. This means that the TAR can strongly underestimate the vertical wave number, and as a consequence, the linear damping, in the sub-inertial regime. This demonstrates the necessity to adopt the NT framework when studying the transport of momentum and matter by GIWs. This is also true in the strongly stratified regime, in which, however, the NT prediction converges very rapidly towards the traditional one (for which N/(2Ω)​ > ​​ > ​1 while ω​ < ​​ < ​N) in the super-inertial regime.

Fig. 2. Ratios ${\widehat{\tau }}_{\mathrm{T}}/{\widehat{\tau }}_{\mathrm{NT}}$ and ${k_V^2}_{\mathrm{T}}/{k_V^2}_{\mathrm{NT}}$ as a function of the normalised frequency ω/2Ω for different colatitudes θ ≡ {0,π/8,π/4,3π/8,π/2} in the weakly (N/2Ω = 5; left panel) and strongly stratified (N/2Ω = 50; right panel) cases.

In Fig. 3 we illustrate these predictions with a complementary point of view. We plot the ratios ${k_V^2}_{\mathrm{T}}/{k_V^2}_{\mathrm{NT}}$ (left raw) and ${\widehat{\tau }}_{\mathrm{T}}/{\widehat{\tau }}_{\mathrm{NT}}$ (right raw) as a function of the ratio N/(2Ω) and of the colatitude θ for a sub-inertial wave for which ω = 1.5Ω (first line) and a super-inertial wave for which ω = 2.5Ω (second line). On the one hand, we again verified that the NT predictions converge to the traditional ones, as expected in the strongly stratified regime in which N/(2Ω)​ > ​​ > ​1 while ω​ < ​​ < ​N (e.g. Friedlander 1987; Mathis et al. 2008; Mathis 2009; Mathis et al. 2014). On the other hand, we again found a strong underestimation of the vertical wave number, and as a consequence, of the linear damping in the sub-inertial regime, as we showed in Fig. 2. In the sub-inertial regime, the strongest underestimation occurs for weak N/(2Ω) ratios close to the critical colatitude θ_c. This might be related to the so-called wave focusing that appears in the NT case when GIWs become strongly focused and sheared close to θ_c. This can lead to strong mixing in weakly stratified oceanic layers (Shrira & Townsend 2010). In the super-inertial regime, the strongest underestimation (with moderate values ∼0.6–0.7) occurs for weak N/(2Ω) ratios at mid-latitudes.

Fig. 3. Ratios ${k_V^2}_{\mathrm{T}}/{k_V^2}_{\mathrm{NT}}$ (left column) and ${\widehat{\tau }}_{\mathrm{T}}/{\widehat{\tau }}_{\mathrm{NT}}$ (right column) as a function of the ratio N/2Ω and of the colatitude θ in the sub-inertial (ω = 1.5Ω; first line) and in the super-inertial (ω = 2.5Ω; second line) regimes. In the sub-inertial regime, GIWs are trapped in an equatorial belt while they are propagating at all colatitudes in the super-inertial regime.

4. Wave breaking

4.1. Convective wave breaking

4.1.1. Formalism

We considered the second mechanism through which GIWs exchange momentum with mean flows: convective wave breaking (CWB). To derive the transported flux of momentum, we followed the method proposed by Lott & Guez (2013) for non-orographic gravity waves and generalised to GIWs by Ribstein et al. (2022) in the traditional case (see also Mathis 2025, for a global deep spherical shell). It is important to underline that the predictions obtained within these theoretical frameworks have been successfully compared to in situ measurements of wave-triggered Reynolds stresses in the stratosphere (Lott et al. 2023). This simultaneously demonstrates the realism and the robustness of the proposed parameterisation and the relevance of the method allowing its derivation. The conditions in the Earth’s atmosphere, where the Prandtl number P_r is close to unity (P_r ≈ 0.71), strongly differ from those of stellar radiation zones in which P_r​ < ​​ < ​1. Because the adiabatic heat equation (Eq. (40)) used to derive the convective breaking condition does not depend on P_r, however, we hope that this approach can also be applied in stellar interiors, but also in oceanic layers in which P_r is above unity (P_r ≈ 7).

As a first step, we considered the heat transport equation expressed as a function of the fluctuation of the temperature (T′),

$$\begin{array}{c}\hfill D_{t}T\prime +\mathrm{\Gamma }w=0\text{with}\mathrm{\Gamma }=\frac{\overline{T}{N}^2}{\overline{g}},\end{array}$$(40)

where $D_t={\partial }_t+\overline{U}{\partial }_x$. Following Lindzen (1981), Lott & Guez (2013), Ribstein et al. (2022), and Mathis (2025), the condition for the wave convective breaking is

$$\begin{array}{c}\hfill |{\partial }_{z}T\prime |>\mathrm{\Gamma }.\end{array}$$(41)

Using the Fourier expansion in time in the linearised heat transport equation (40) and considering low-frequency GIWs for which the JWKB approximation can be used, we derived the saturation amplitude of the vertical velocity for which |∂_zT′| = Γ,

$$\begin{array}{c}\hfill |\widehat{W}{|}_{\mathrm{sat}}^{\mathrm{CWB}}=\frac{\widehat{\omega }}{k_V},\end{array}$$(42)

above which GIWs are subject to convective breaking. We obtained a result similar in its mathematical form to both Lott & Guez (2013) and Mathis (2025). Moreover, as pointed out by Mathis (2025), this result is equivalent to the breaking criteria proposed by Press (1981), Goodman & Dickson (1998), and Barker & Ogilvie (2010). This can be explained by the fact that the Coriolis acceleration does not intervene directly in the heat transport equation; its action is through the modification of the vertical wave number k_V. This allowed us, using Eq. (23), to compute the carried flux of energy along the vertical direction,

$$\begin{array}{c}\hfill F_{\mathrm{E}}^{\mathrm{CWB}}=-\frac{1}{2}\epsilon \overline{\rho }\frac{(f^2-{\widehat{\omega }}^2)}{\widehat{\omega }{k}_{\perp }^2}{k}_V{\left(\frac{\widehat{\omega }}{k_V}\right)}^2.\end{array}$$(43)

Finally, we derived the momentum flux along this direction using Eq. (24),

$$\begin{array}{cc}\hfill F_{\mathrm{AM}}^{\mathrm{CWB}}& =-\frac{k_x}{\widehat{\omega }}{F}_{\mathrm{E}}\hfill \\ \hfill & =\frac{1}{2}\epsilon \overline{\rho }\left(\frac{f^2-{\widehat{\omega }}^2}{{\widehat{\omega }}^2}\right)cos\alpha \frac{k_{\mathrm{V}}}{k_{\perp }}{\left(\frac{\widehat{\omega }}{k_V}\right)}^2,\hfill \end{array}$$(44)

where k_x = k_⊥cos α. Introducing

$$\begin{array}{cc}\hfill F_r& =\frac{\widehat{\omega }}{N},R_o=\frac{\widehat{\omega }}{2\mathrm{\Omega }},\hfill \\ \hfill R_{o;V}& =\frac{\widehat{\omega }}{f}=\frac{1}{cos\theta }{R}_o,\text{and}{R}_{o;H}=\frac{\widehat{\omega }}{\stackrel{\sim }{f}}=\frac{1}{sin\theta }{R}_o,\hfill \end{array}$$(45)

we write

$$\begin{array}{cc}\hfill F_{\mathrm{AM}}^{\mathrm{CWB}}& =\frac{1}{2}\epsilon \overline{\rho }(R_{o;V}^{-2}-1)cos\alpha \hfill \\ \hfill & \times {[\frac{F_r^{-2}-1}{1-R_{o;V}^{-2}}+\frac{R_{o;H}^{-2}{sin}^2\alpha }{{(1-R_{o;V}^{-2})}^2}]}^{-1/2}{\left(\frac{\widehat{\omega }}{k_{\perp }}\right)}^2;\hfill \end{array}$$(46)

sub-inertial (super-inertial) waves correspond to R_o < 1 (R_o > 1).

In the non-rotating case, this reduces to the result obtained by Lott & Guez (2013),

$$\begin{array}{c}\hfill |F_{\mathrm{AM}}^{\mathrm{CWB}}(\mathrm{\Omega }=0)|=\frac{1}{2}\overline{\rho }cos\alpha F_r{\left(\frac{\widehat{\omega }}{k_{\perp }}\right)}^2.\end{array}$$(47)

Using the expression for the flux of momentum transported along the vertical direction given in Eq. (27), we also explicitly computed the ratio of its value in the NT case by the one in the traditional case,

$$\begin{array}{c}\hfill \frac{F_{\mathrm{AM};\mathrm{NT}}}{F_{\mathrm{AM};\mathrm{T}}}=\frac{{k_V}_{\mathrm{NT}}}{{k_V}_{\mathrm{T}}}{\left(\frac{|{\widehat{W}}_{\mathrm{NT}}|}{|{\widehat{W}}_{\mathrm{T}}|}\right)}^2.\end{array}$$(48)

Using the saturation amplitude obtained in Eq. (42), we obtained

$$\begin{array}{c}\hfill \frac{F_{\mathrm{AM};\mathrm{NT}}^{\mathrm{CWB}}}{F_{\mathrm{AM};\mathrm{T}}^{\mathrm{CWB}}}=\left(\frac{{k_V}_{\mathrm{T}}}{{k_V}_{\mathrm{NT}}}\right)<1.\end{array}$$(49)

4.1.2. Exploration of the parameter space

In Fig. 4 we examine the ratio of the momentum flux carried at the convective breaking by GIWs to the flux carried by gravity waves (corresponding to the non-rotating case Ω ≡ 0) as a function of the colatitude θ for fixed wave Froude ($F_{\mathrm{r}}=\widehat{\omega }/N$) and Rossby ($R_{\mathrm{o}}=\widehat{\omega }/\left(2\mathrm{\Omega }\right)$) numbers in the NT (left panel) and in the traditional (right panel) cases. For decreasing wave Rossby numbers (i.e. for increasing rotation at fixed frequency) these ratios decrease for both cases. This can be interpreted keeping in mind that rotation inhibits the convective instability and the triggered transport of energy and as a consequence of momentum (Stevenson 1979; Augustson & Mathis 2019). In the sub-inertial regime (i.e. for R_o < 1), the flux of angular momentum vanishes outside an equatorial belt for θ < θ_c and θ > π − θ_c. This corresponds to the equatorial trapping of GIWs by the Coriolis acceleration in the sub-inertial regime (e.g. Dintrans & Rieutord 2000; Gerkema & Shrira 2005; Mathis et al. 2014). The key difference between the traditional and the NT regime is the artificial convergence of the ratio to unity at the equator predicted by the TAR. This convergence is artificial because the TAR cannot be applied at the equator: the projection of the rotation vector along the vertical direction vanishes there. The correct prediction is therefore obtained in the NT framework when the full Coriolis acceleration is taken into account. In this regime, the carried momentum flux diminishes at all the colatitudes when compared to the non-rotating case, in particular, at the equator. This difference between the NT and the traditional cases is evaluated in Fig. 5, where we plot the ratio of the flux predicted in the NT regime and of the one predicted in the traditional regime as predicted in Eq. (49). In the sub-inertial regime, it is well below unity since k_VT < k_VNT. The reason is that when the full Coriolis acceleration is taken into account, the convective instability of the waves and the related transport is inhibited more effectively. In the super-inertial regime, the ratio converges to unity. This agrees with the results obtained in the previous section where we studied linear radiative damping, where the NT prediction converges to the traditional one when R₀ → ∞. As a general conclusion, it is therefore mandatory to use the NT framework to study wave-mean flow interactions, in particular, in the equatorial regions. On the one hand, these regions are of particular importance in planetary atmospheres where these interactions play a key role in the general circulation, for instance, with the quasi-biennial oscillation (e.g. Baldwin et al. 2001, and references therein). On the other hand, for rapidly rotating stars, sub-inertial GIWs trapped in the equatorial belt efficiently deposit angular momentum below the stellar surface, which can lead to mass loss (Ando 1985; Rogers 2015; Neiner et al. 2020).

Fig. 4. Ratio of the vertical flux of momentum computed with the full Coriolis acceleration (left panel) and assuming the TAR (right panel) with its value computed in the non-rotating case (with Ω ≡ 0) as a function of the colatitude (θ) for a fixed value of the Froude wave number (Fr = ω/N = 0.25) and different values of the wave Rossby number (Ro = {0.1,0.5,1,2,100}). In the NT case, we fixed α = π/4.

Fig. 5. Ratio of the vertical flux of momentum computed with the full Coriolis acceleration with its value computed in the traditional case as a function of the colatitude (θ) for a fixed value of the wave Froude number (Fr = ω/N = 0.25) and different values of the wave Rossby number (Ro = {0.1,0.5,1,2,100}). In the NT case, we fixed α = π/4.

4.2. The shear-induced GIW breaking

As shown by Garcia Lopez & Spruit (1991) and Mathis (2025), internal gravity waves and GIWs are not only subject to convective breaking. They can become unstable to their own vertical shear, which can lead to their breaking. We call this shear-driven wave breaking (SWB). Following Mathis (2025), we adapted the method built by Lott et al. (2012) and Lott & Guez (2013) for the study of the CWB to the case of the SWB. To do this, we considered the instability criteria for the vertical shear of GIWs from which we derived an estimate of the saturation amplitude for the vertical component of the velocity and the corresponding flux of momentum transported along the vertical direction.

Before we derived these quantities, it is important to report the progresses achieved recently in the study of the instability of a vertical shear in a stably stratified rotating fluid with the full Coriolis acceleration. In a recent theoretical article, Park & Mathis (2025) have demonstrated that this vertical shear is submitted to two different types of instabilities: the inflectional instability, which corresponds to the classical vertical shear instability studied in stellar physics (e.g. Zahn 1983, 1992; Maeder & Meynet 1996; Prat & Lignières 2013, 2014; Garaud & Kulenthirarajah 2016), and the inertial instabilities, such as the instability of Rayleigh-unstable differential rotation profiles and the Goldreich-Schubert-Fricke instability in stellar physics (e.g. Goldreich & Schubert 1967; Fricke 1968; Barker et al. 2020; Dymott et al. 2023). Park & Mathis (2025) demonstrated that the inflectional instability mostly develops around the equator for non-axisymmetric perturbations, while baroclinic instabilities develop at mid-latitudes for stably stratified rotating fluids with a weak heat diffusion and around the poles for stably stratified rotating fluids with a strong heat diffusion, as in the case of stellar radiation zones. We focused on the interactions of GIWs with zonal flows, in particular, in the sub-inertial regime in which GIWs are trapped in a belt around the equator. We thus chose to focus on the vertical-shear inflectional instability, while vertical-shear inertial instabilities are important for the potential breaking of Rossby waves that can propagate at high latitudes. One of the key results by Park & Mathis (2025) is that the theoretical inflectional instability criteria around the equator derived within a linear stability analysis are close to the criterion for the non-rotating case. This was also observed in direct numerical simulations around the equator of vertical shear instabilities by Chang & Garaud (2021). We therefore considered the following instability criteria:

$$\begin{array}{c}\hfill \frac{N^2}{{\left|{\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}\left[\sqrt{{\langle \mathrm{Re}\left({\mathit{u}}_H\right)·\mathrm{Re}\left({\mathit{u}}_H\right)\rangle }_H}\right]}\right|}^2}<R{i}_{\mathrm{c}}\text{and}\end{array}$$(50)

$$\begin{array}{c}\hfill \frac{N^2}{{\left|{\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}\left[\sqrt{{\langle \mathrm{Re}\left({\mathit{u}}_H\right)·\mathrm{Re}\left({\mathit{u}}_H\right)\rangle }_H}\right]}\right|}^2}\left(\frac{vl}{K}\right)<R{i}_{\mathrm{c}}\end{array}$$(51)

for stably stratified rotating fluids that are the seat of weak and strong heat diffusion, respectively. We recall that ${\langle ···\rangle }_H$ is the horizontal average with ${\langle \mathrm{Re}\left(\mathit{X}\right)·\mathrm{Re}\left(\mathit{X}\right)\rangle }_H=1/2(\mathit{X}·{\mathit{X}}^{\ast })$, and we introduce $u_H=\sqrt{{\mathit{u}}_H·{\mathit{u}}_H^{\ast }}$, the amplitude of the horizontal component of GIWs, and v and l the characteristic vertical velocity and length scale that must be considered in the regime of strong heat diffusion. As in Mathis (2025), we identify that $v\equiv \sqrt{{\langle \mathrm{Re}\left({\mathit{u}}_V\right)·\mathrm{Re}\left({\mathit{u}}_V\right)\rangle }_H}$, the amplitude of the vertical component of the studied GIW, and $l\equiv 2\pi k_V^{-1}$ is its vertical wave number. In the sections below, each of these regimes is studied.

4.2.1. The regime of weak heat diffusion

This first regime corresponds to the geophysical case where the heat diffusion does not dominate the viscous diffusion (we recall that the Prandtl number, P_r, is of the order of unity in the Earth’s atmosphere, but above unity in the oceans, with P_r ∼ 7). For low-frequency GIWs for which the JWKB approximation can be used, the Richardson criteria given in Eq. (50) becomes:

$$\begin{array}{c}\hfill \frac{N^2}{1/2k_V^2{\left|u_H\right|}^2}<R{i}_{\mathrm{c}}.\end{array}$$(52)

Following the method presented by Mathis (2025), we express the horizontal average of the horizontal component of the velocity of GIWs using polarisation relations, Eqs. (A.7) and (A.8), and the JWKB approximation. We obtain

$$\begin{array}{c}\hfill {\left|u_H\right|}^2={\left|\widehat{U}\right|}^2+{\left|\widehat{V}\right|}^2={\left(\frac{k_V}{k_{\perp }}\right)}^2\frac{f^2+{\widehat{\omega }}^2}{{\widehat{\omega }}^2}{\left|\widehat{W}\right|}^2.\end{array}$$(53)

The threshold instability condition allowed us to derive the expression for the saturation amplitude of the vertical velocity of GIWs when their vertical shear triggers their breaking,

$$\begin{array}{c}\hfill {\left({\left|\widehat{W}\right|}^2\right)}_{\mathrm{sat}}^{\mathrm{SWB}}=\frac{2}{R{i}_{\mathrm{c}}}\frac{N^2}{f^2+{\widehat{\omega }}^2}{\left(\frac{k_{\perp }}{k_V}\right)}^2{\left(\frac{\widehat{\omega }}{k_V}\right)}^2.\end{array}$$(54)

In the non-rotating case, we verified that we recovered the result obtained by Mathis (2025) in the simplest case of internal gravity waves,

$$\begin{array}{c}\hfill {\left({\left|\widehat{W}\right|}^2\right)}_{\mathrm{sat}}^{\mathrm{SWB}}(\mathrm{\Omega }=0)=\frac{2}{R{i}_{\mathrm{c}}}{\left(\frac{\widehat{\omega }}{k_V}\right)}^2.\end{array}$$(55)

Using again the expression for the flux of momentum transported along the vertical direction given in Eq. (27), we computed the ratio of its value in the NT case by the one in the traditional case,

$$\begin{array}{c}\hfill \frac{F_{\mathrm{AM};\mathrm{NT}}}{F_{\mathrm{AM};\mathrm{T}}}=\frac{{k_V}_{\mathrm{NT}}}{{k_V}_{\mathrm{T}}}{\left(\frac{|{\widehat{W}}_{\mathrm{NT}}|}{|{\widehat{W}}_{\mathrm{T}}|}\right)}^2.\end{array}$$(56)

Using the saturation amplitude obtained in Eq. (54), we obtained

$$\begin{array}{c}\hfill \frac{F_{\mathrm{AM};\mathrm{NT}}^{\mathrm{SWB}}}{F_{\mathrm{AM};\mathrm{T}}^{\mathrm{SWB}}}={\left(\frac{{k_V}_{\mathrm{T}}}{{k_V}_{\mathrm{NT}}}\right)}^3<1.\end{array}$$(57)

As for the CWB, we therefore obtained a predicted flux when the full Coriolis acceleration was taken into account that is weaker than the flux predicted within the TAR. Since k_VT < k_VNT, the shear instability can be favoured in the NT case, but with a saturation amplitude for the vertical velocity that is weaker than in the traditional prediction, which leads to a weaker momentum flux. The change in the power of k_VT/k_VNT of the ratio $F_{\mathrm{AM};\mathrm{NT}}^{\mathrm{SWB}}/F_{\mathrm{AM};\mathrm{T}}^{\mathrm{SWB}}$ from unity obtained for the CWB to 3 is due to the modulation factor $N^2/(f^2+{\widehat{\omega }}^2){(k_{\perp }/k_V)}^2$ in Eq. (54).

4.2.2. The regime of strong heat diffusion

This regime corresponds to the case of the bulk of stellar radiative regions, where the Prandtl number is small. We again considered low-frequency GIWs for which the JWKB approximation can be applied. Following the same method as in the weak heat diffusion regime, the vertical shear instability criteria (51) becomes

$$\begin{array}{c}\hfill \frac{N^2}{1/2k_V^2{\left|u_H\right|}^2}\frac{1/2{\left|\widehat{W}\right|}^{2}2\pi k_V^{-1}}{K}<R{i}_{\mathrm{c}}.\end{array}$$(58)

With the polarisation relation Eq. (53), this leads to the saturation amplitude for the vertical component of the velocity,

$$\begin{array}{c}\hfill {\left|\widehat{W}\right|}_{\mathrm{sat}}^{\mathrm{SWB}}=\frac{2\pi }{R{i}_{\mathrm{c}}}\frac{N^2}{f^2+{\widehat{\omega }}^2}{\left(\frac{k_{\perp }}{k_V}\right)}^2\frac{\widehat{\omega }}{K{k}_V^2}\frac{\widehat{\omega }}{k_V}.\end{array}$$(59)

Using Eq. (48), we finally obtained the ratio of the value of the momentum flux transported along the vertical direction predicted in the NT case by the one in the traditional case,

$$\begin{array}{c}\hfill \frac{F_{\mathrm{AM};\mathrm{NT}}^{\mathrm{SWB}}}{F_{\mathrm{AM};\mathrm{T}}^{\mathrm{SWB}}}={\left(\frac{{k_V}_{\mathrm{T}}}{{k_V}_{\mathrm{NT}}}\right)}^9<1.\end{array}$$(60)

It is thus more strongly overestimated in the traditional regime for shear-driven breaking of GIWs in the strong heat diffusion stellar regime than in their shear-driven breaking in the weak heat diffusion atmospheric or oceanic regimes and in the case of their convective breaking.

5. Wave transport parameterisation

We have demonstrated that the full Coriolis acceleration must be taken into account in the general case when the conditions for the validity of the TAR are not satisfied. When this approximation is applied outside of its range of applicability, the linear damping of GIWs in the sub-inertial regime is underestimated and predicts a momentum deposition at an altitude at which it has been already damped in reality, while the transported momentum flux by GIW convective and shear-driven breakings is overestimated. As when rotation is ignored (Lott & Guez 2013) and in the traditional case (Ribstein et al. 2022), we thus have to provide a parametrisation of momentum transport by GIWs. Following Lott & Guez (2013), we thus derived the transported momentum flux,

$$\begin{array}{cc}\hfill F_{\mathrm{AM}}\left(z\right)& =-\frac{k_x}{\widehat{\omega }}{F}_{\mathrm{E}}\left(z\right)\hfill \\ \hfill & ={\displaystyle \sum _{\omega ,k_{\perp }}}\min \{F_{\mathrm{AM}}(z=z_{\mathrm{o}})exp[-{\tau }_{\mathrm{G}}\left(z\right)],\hfill \\ \hfill & F_{\mathrm{AM};\mathrm{NT}}^{\mathrm{CWB}}\left(z\right),F_{\mathrm{AM};\mathrm{NT}}^{\mathrm{SWB}}\left(z\right)\}.\hfill \end{array}$$(61)

It is given by the minimum value between the classical quasi-adiabatic flux that applies during the GIW propagation before the altitude of the breaking and the saturated values of the flux at this place and after.

This formalism works as follows. We considered a given altitude z. On the one hand, if the quasi-adiabatic velocity of a GIW is lower than the breaking saturation velocities (i.e. before the altitude at which the wave breaks), the corresponding quasi-adiabatic momentum flux is lower than the flux carried by the breaking of the wave, and the min function therefore selects the quasi-adiabatic value. On the other hand, if the quasi-adiabatic velocity of a GIW has formally become higher than one of the breaking saturation velocities (i.e. after the wave has broken), the quasi-adiabatic flux becomes formally higher than the flux carried by the breaking, which is naturally selected by the min function. This allowed us to compute the variation in the mean zonal flow $\overline{U}$ that varies with altitude, following Eq. (26).

In our formalism, we assumed that GIWs can be convectively unstable and that their own vertical shear can be unstable (Thorpe 1999; Sutherland 2001). The formalism is based on local linear stability criteria and selects the instability with the lowest threshold velocity, leading to one of the instabilities and related breaking. It was shown in the literature, however, that (i) internal gravity waves and GIWs can be submitted to other types of instabilities, such as parametric sub-harmonic instabilities (e.g. Ghaemsaidi & Mathur 2019), (ii) the trigger of the different instabilities depends on the properties of the initial perturbations (e.g. Lombard & Riley 1996; Liu et al. 2010), and (iii) several wave instabilities (e.g. the convective and the shear-induced instabilities) can occur simultaneously in a flow (e.g. Lombard & Riley 1996; Howland et al. 2021). All these mechanisms are related to subtle wave-wave non-linear interactions. This shows that devoted non-linear numerical simulations must be undertaken to better characterise these complex non-linear mechanisms and improve their parametrisation for a long timescale evolution (e.g. Fruman et al. 2014; Gervais et al. 2021).

6. Discussion and conclusions

We have built a prototype local Cartesian model to study the interactions between gravito-inertial (inertia-gravity) waves and stellar (or planetary and geophysical) mean flows. This allowed us to treat any stratification, mean rotation rate, heat diffusion, and viscosity. We considered two of the key mechanisms that allow waves to exchange momentum with mean flows: their linear damping by viscosity and heat diffusion (e.g. Press 1981; Zahn et al. 1997; Vadas & Fritts 2005), and their non-linear breaking because of their convective (Lindzen 1981; Lott et al. 2012; Lott & Guez 2013) and shear (e.g. Garcia Lopez & Spruit 1991; Mathis 2025) instabilities. Our local model can be considered as a first predictive tool of what occurs in regions in stellar, planetary, atmospheric, and oceanic layers close to critical layers where the linear damping of the waves drastically increases until they break non-linearly. For the linear damping, we demonstrated that it increases when we compared the rotating case to the non-rotating case, as predicted in the literature (e.g. Pantillon et al. 2007; Mathis et al. 2008). This is due to the increase in the wave number. The key new interesting result, obtained here when the heat diffusion dominated (i.e. for stellar radiation zones), is that the linear damping predicted using the traditional approximation is completely underestimated in the sub-inertial regime when compared to the NT prediction, in which the full Coriolis acceleration is taken into account. This means that when the buoyancy associated with the stable stratification does not dominate the component of the Coriolis acceleration in the direction of the entropy or chemical stratification, we must abandon the traditional approximation, which would predict a deposit of momentum far away from the excitation region of waves than it is in reality. This result is of great importance when studying transport and mixing in weakly stably stratified rotating stellar radiative regions in rapidly rotating young late-type stars and early-type stars. The same conclusion applies for weakly stratified atmospheric and oceanic layers where the damping is driven by the combination of heat diffusion and viscosity or dominated by viscosity, respectively (note that a similar conclusion has been obtained recently by Slepyshev 2025). When considering the convective and shear-induced non-linear breaking of the waves, we found that the saturation velocity at which the transition to breaking occurs decreases because the wave number increases. This leads to a decrease in the efficiency of the momentum transport. For the CWB, this can be explained by the loss of efficiency of convection in rapidly rotating systems (e.g. Stevenson 1979; Barker et al. 2014; Augustson & Mathis 2019; Currie et al. 2020). For the shear-induced breaking, the inflectional instability of a vertically sheared flow is localised around the equator with an instability criterion similar to the non-rotating case (Park & Mathis 2025). The efficiency loss of the induced momentum transport in this case then arises because the NT vertical wave vector is larger than the traditional one outside the equator and the poles. Once again, we note a strong difference between the predictions obtained using the traditional approximation or going beyond it. The traditional approximation overestimates the momentum transport because it underestimates the vertical wave number, which increases because of the NT terms. This result is of great importance again for rapidly rotating stars where GIWs propagate, while these waves are key players that drive and regulate the atmospheric and oceanic general circulation on our Earth (e.g. MacKinnon et al. 2017; Sutherland et al. 2019; Ribstein et al. 2022; Achatz et al. 2024). Key prospects of this work are now to combine this complete NT modelling to a corresponding wave-excitation model (Mathis et al. 2014; Augustson et al. 2020) to be able to generalise it to global geometry, and to take magnetic fields into account (e.g. Rogers & MacGregor 2011; Mathis & de Brye 2012). These would be key developments for models of the stellar structure and secular evolution and for atmospheric/oceanic general circulation models in which transport of momentum and macroscopic mixing of chemicals have to be parameterised in a robust and tractable way.

Acknowledgments

S.M. warmly thanks the referee, Pr. A. Barker, for his constructive and detailed report, which allows him to improve the initial manuscript. S.M. acknowledges support from the European Research Council (ERC) under the Horizon Europe programme (Synergy Grant agreement 101071505: 4D-STAR), from the CNES SOHO-GOLF and PLATO grants at CEA-DAp, and from PNPS (CNRS/INSU). While partially funded by the European Union, views and opinions expressed are however those of the author only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

References

Appendix A: Polarisation relations

The main goal of our so-called prototype local model is to quantify the GIWs - zonal mean flow interactions and the related flux of energy and momentum these waves are transporting. To reach this objective, we need to compute the product of the vertical velocity of GIWs by their fluctuation of pressure. To do so in an efficient way, the best method is to derive the expression of the horizontal components of the velocity, the fluctuation of pressure and the buoyancy as a function of the vertical component of the velocity through the so-called polarisation relations.

To derive these latters, we follow Auclair-Desrotour et al. (2018) and we express each component of the velocity and each perturbation of a thermodynamical quantity as this has been done in Eqs. (10) and (12) for the vertical velocity:

$$\begin{array}{c}\hfill \stackrel{\sim }{x}(\mathit{r},t)=\mathrm{Re}\left\{X(\chi ,z){\text{e}}^{\text{i}\omega t}\right\}\text{where}X=\widehat{X}(z)exp\left[\text{i}{k}_{\perp }(\chi +\stackrel{\sim }{\delta }z)\right]\end{array}$$(A.1)

and Re is the real part of a complex number. Since in the NT case, the Poincaré equation is not separable, $\widehat{X}$ carries only part of the vertical dependence of the field $\stackrel{\sim }{x}$, the other part being included in $exp\left(\text{i}{k}_{\perp }\stackrel{\sim }{\delta }z\right)$. A separable solution will be $X=\widehat{X}(z)exp\left(\text{i}{k}_{\perp }\chi \right)$ without any dependence of the vertical function on k_⊥. For a global model going beyond the local model presented here, the term $exp\left(\text{i}{k}_{\perp }\stackrel{\sim }{\delta }z\right)$ corresponds to the mixed derivative ∂_z, χ in Eq. (9), which is responsible for the non-separability of the Poincaré equation in the general case.

The linearised momentum equation, Eq. (6), becomes in its inviscid limit:

$$\begin{array}{cc}& \text{i}\omega U-fV+\stackrel{\sim }{f}W=-\frac{\text{i}{k}_{\perp }cos\alpha }{\overline{\rho }}P,\hfill \end{array}$$(A.2)

$$\begin{array}{cc}& \text{i}\omega V+fU=-\frac{\text{i}{k}_{\perp }sin\alpha }{\overline{\rho }}P,\hfill \end{array}$$(A.3)

$$\begin{array}{cc}& \text{i}\omega W-\stackrel{\sim }{f}U=-\frac{1}{\overline{\rho }}\frac{\partial P}{\partial z}+B,\hfill \end{array}$$(A.4)

while we obtain for the equations of conservation of mass (Eq. 7) and energy in its adiabatic limit (Eq. 7), respectively:

$$\begin{array}{cc}& \text{i}{k}_{\perp }cos\alpha U+\text{i}{k}_{\perp }sin\alpha V+\frac{\partial W}{\partial z}=0,\hfill \end{array}$$(A.5)

$$\begin{array}{cc}& \text{i}\omega B+N^{2}W=0.\hfill \end{array}$$(A.6)

Using Eq. (A.1), we have

$$\begin{array}{c}\hfill \frac{\partial W}{\partial z}=({\widehat{W}}^{\prime }+\text{i}{k}_{\perp }\stackrel{\sim }{\delta }\widehat{W})exp\left[\text{i}{k}_{\perp }(\chi +\stackrel{\sim }{\delta }z)\right],\end{array}$$

where the prime denotes differentiation with respect to z, and all the fields can be expressed in term of $\widehat{W}$ as

$$\begin{array}{cc}\hfill \widehat{U}& =\frac{{\stackrel{\sim }{f}}_{\text{s}}(fcos\alpha -\text{i}\omega sin\alpha )}{f^2-{\omega }^2}\widehat{W}+\frac{fsin\alpha +\text{i}\omega cos\alpha }{\omega k_{\perp }}{\widehat{W}}^{\prime },\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill \widehat{V}& =\frac{{\stackrel{\sim }{f}}_{\text{s}}(fsin\alpha +\text{i}\omega cos\alpha )}{f^2-{\omega }^2}\widehat{W}-\frac{fcos\alpha -\text{i}\omega sin\alpha }{\omega k_{\perp }}{\widehat{W}}^{\prime },\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill \widehat{P}& =\text{i}\frac{{\rho }_0{\stackrel{\sim }{f}}_{\text{c}}}{k_{\perp }}\widehat{W}+\text{i}\frac{{\rho }_0(f^2-{\omega }^2)}{\omega k_{\perp }^2}{\widehat{W}}^{\prime },\hfill \end{array}$$(A.9)

$$\begin{array}{cc}\hfill \widehat{B}& =\text{i}\frac{N^2}{\omega }\widehat{W},\hfill \end{array}$$(A.10)

where ${\stackrel{\sim }{f}}_{\text{c}}=\stackrel{\sim }{f}cos\alpha$. Then, each field has to be multiplied by $exp\left[\text{i}(k_{\perp }(\chi +\stackrel{\sim }{\delta }z)+\omega t)\right]$ to get the complete solution. The physical solution is obtained by taking the real part of the complete complex solution. For low-frequency GIWs, for which the JWKB assumption can be assumed, we thus obtain

$$\begin{array}{c}\hfill \widehat{P}\approx -\epsilon \overline{\rho }\frac{(f^2-{\omega }^2)}{\omega k_{\perp }^2}{k}_V\widehat{W}.\end{array}$$(A.11)

Appendix B: Comparison with the non-rotating regime

In the main text, we present and discuss the ratios k_V²_T/k_V²_NT and ${\widehat{\tau }}_{\mathrm{T}}/{\widehat{\tau }}_{\mathrm{NT}}$ to identify when it is necessary to take into account the full Coriolis acceleration. It is also interesting to show as a complement in Figs. A.1 and A.2 the ratio k_V²_NT/k_V²(Ω = 0), k_V²_T/k_V²(Ω = 0), ${\widehat{\tau }}_{\mathrm{NT}}/\widehat{\tau }(\mathrm{\Omega }=0)$ and ${\widehat{\tau }}_{\mathrm{T}}/\widehat{\tau }(\mathrm{\Omega }=0)$ respectively in the weakly (left panels) and strongly (right panels) stratified cases in which S = N/2Ω = 5 and 50, respectively. This allows us to identify how both k_V² and $\widehat{\tau }$ are increasing at low-frequencies with rotation. This leads to linear dampings closer to the excitation region and to transported vertical flux of momentum by convective and shear-induced wave breakings weaker in the rotating case than in the non-rotating case (we refer the reader to Sects. 3 and 4, respectively). These trends are predicted both when assuming the TAR and in a fully NT framework. Finally, we also recover that ${k_V^2}_{\mathrm{T}}/{k_V^2}_{\mathrm{NT}}<1$ and ${\widehat{\tau }}_{\mathrm{T}}/{\widehat{\tau }}_{\mathrm{NT}}<1$.