The effects of non-linearities on tidal flows in the convective envelopes of rotating stars and planets in exoplanetary systems

ABSTRACT

In close exoplanetary systems, tidal interactions drive orbital and spin evolution of planets and stars over long time-scales. Tidally forced inertial waves (restored by the Coriolis acceleration) in the convective envelopes of low-mass stars and giant gaseous planets contribute greatly to the tidal dissipation when they are excited and subsequently damped (e.g. through viscous friction), especially early in the life of a system. These waves are known to be subject to non-linear effects, including triggering differential rotation in the form of zonal flows. In this study, we use a realistic tidal body forcing to excite inertial waves through the residual action of the equilibrium tide in the momentum equation for the waves. By performing 3D non-linear hydrodynamical simulations in adiabatic and incompressible convective shells, we investigate how the addition of non-linear terms affects the tidal flow properties, and the energy and angular momentum redistribution. In particular, we identify and justify the removal of terms responsible for unphysical angular momentum evolution observed in a previous numerical study. Within our new set-up, we observe the establishment of strong cylindrically sheared zonal flows, which modify the tidal dissipation rates from prior linear theoretical predictions. We demonstrate that the effects of this differential rotation on the waves neatly explains the discrepancies between linear and non-linear dissipation rates in many of our simulations. We also highlight the major role of both corotation resonances and parametric instabilities of inertial waves, which are observed for sufficiently high tidal forcing amplitudes or low viscosities, in affecting the tidal flow response.

planet-star interactions, hydrodynamics, stars: low-mass, planets and satellites: gaseous planets, waves, instabilities

1 INTRODUCTION

Over the last 25 yr or so, nearly 5000 exoplanets have been discovered around mostly late-type stars from M to F spectral types, thanks to spatial and ground-based observational campaigns (see e.g. Perryman 2018, for a review). Several hundred of these exoplanets,¹ known as Hot Jupiters, are giant gaseous planets orbiting in less than 10 d around their host stars. In such compact systems, tides inside the star and the planet are very likely to have played a major role in modifying their rotational periods, their obliquities (also known as spin–orbit angles) and shaping the orbital architecture of the planet (typically the eccentricity and the semimajor axis; see e.g. Ogilvie 2014, 2020; Mathis 2019). In addition to these secular dynamical modifications, star–planet tidal interactions can also bring about internal structural changes by tidal heating, and are for example suspected to affect the size of the planet by inducing radius inflation (e.g. Millholland 2019, for close sub-Neptune planets).

One component of the tidal response in a body due to the tidal gravitational potential of its perturber (the star or the planet) is the quasi-hydrostatic tidal bulge, and its associated large-scale flow, which is the so-called equilibrium tide (Love 1911; Zahn 1966a, 1989; Eggleton, Kiseleva & Hut 1998; Hansen 2010, 2012). Signatures of ellipsoidal deformations and equilibrium tides have, for example, been detected by photometry at the surfaces of HAT-P 7, WASP-18, and WASP-12, around the latter of which orbits a Hot Jupiter undergoing orbital decay (Welsh et al. 2010; Shporer et al. 2019; Maciejewski et al. 2016, 2020; Yee et al. 2020; Turner, Ridden-Harper & Jayawardhana 2021). The inward migration of WASP-12 b is thought to be due to the second type of tides, the dynamical tides, which are the oscillatory responses to the tidal potential that arise when considering the non-hydrostatic (i.e. inertial) terms in the equations of motion for tidal flows (Zahn 1966c, 1970, 1975; Ogilvie & Lin 2004). In the case of WASP-12, tidally excited gravity waves (which are restored by the buoyancy force) are thought to propagate within the radiative zone, and if the star has a radiative core they can geometrically focus towards the stellar centre, and break non-linearly, thus dissipating their kinetic energy (e.g. Goodman & Dickson 1998; Ogilvie & Lin 2007; Barker & Ogilvie 2010, 2011; Weinberg et al. 2017; Barker 2020). The exchange of angular momentum with the planet due to the tidal dissipation inside the star is what is causing the planetary infall in this picture. However, this dissipative mechanism is likely to be most efficient in systems featuring old (towards the end of the main sequence) and slowly rotating stars hosting a massive planet (Barker 2020; Lazovik 2021).

Another powerful mechanism to transfer angular momentum is the dissipation of inertial waves (restored by the Coriolis acceleration), which can propagate in the convective envelopes of low-mass stars and giant gaseous planets (e.g. Ogilvie & Lin 2004, 2007; Wu 2005a, b). The dissipation of tidal inertial waves in stars is particularly efficient in the early stages of the life of a star (typically during the pre-main sequence of late-type stars) as was shown for example by Mathis (2015), Barker (2020), and Lazovik (2021), using the frequency-averaged formalism developed by Ogilvie (2013). It can then have major consequences, for example, on the planetary semimajor axis, hastening the infall of the planet if it orbits inside the corotation radius (where the planetary orbital period equals the rotational period of the star), or contribute significantly to pushing it away in the opposite case (in comparison with dissipation of the equilibrium tide alone, see also Bolmont & Mathis 2016; Gallet et al. 2017). It can also explain the circularization periods of solar-type binary stars, which was a long-standing theoretical problem (Barker 2022). Inertial waves excited inside planets are also believed to be important in explaining the preferentially circular orbits of the shortest period hot Jupiters through tidally driven orbital circularization (Ogilvie & Lin 2004; Barker 2016).

However, giving robust predictions for the dynamical evolution of a close star–planet system is particularly challenging. The assessment of the tidal dissipation has been shown to significantly depend on many parameters like the stellar angular velocity and metallicity, and the masses of both the star and the planet (e.g. Mathis 2015; Gallet et al. 2017; Bolmont et al. 2017; Barker 2020; Lazovik 2021). Moreover, the current formalism used for tidal inertial waves is based on a frequency average of the tidal dissipation, even though the dissipation is known to be strongly dependent on the tidal forcing frequency in the linear regime when these waves are excited (Ogilvie & Lin 2004; Ogilvie 2005, 2009; Rieutord & Valdettaro 2010). Although this formalism is an efficient and rather straightforward tool to provide information on the long-term tidal evolution of exoplanetary (and close binary) systems, the frequency-dependent tidal dissipation gives us an instant picture of tidal effects at a given evolutionary stage of the system featuring specific tidal forcing frequencies. The resonance locking mechanism, for instance, relies on the continuous match over time between the forcing frequencies of the perturber and some resonant and highly dissipative modes inside the perturbed body (Witte & Savonije 1999; Fuller, Luan & Quataert 2016). This mechanism has been successfully applied to the Saturnian system to provide an explanation for the high tidal dissipation rate in Saturn inferred from the rapid orbital expansion measured for Titan (using inertial modes, with important consequences for the formation scenario of Saturn’s satellites, Lainey et al. 2020) and also applied to extra-solar giant gaseous planets orbiting F-type stars with convective cores (using gravity modes, Ma & Fuller 2021). Furthermore, the nature, propagation, and dissipation properties of tidal flows are also modified by the magnetism and differential rotation both present in the convective envelopes of low-mass stars and giant gaseous planets (Baruteau & Rieutord 2013; Guenel et al. 2016b, a; Wei 2016, 2018; Lin & Ogilvie 2018). Different powerful and frequency-dependent tidal dissipation mechanisms can then emerge like the Ohmic dissipation of magnetic energy overtaking the viscous dissipation of kinetic energy in various low-mass stars throughout their lives (see Astoul et al. 2019), or enhanced dissipation at corotation resonances (or critical layers, where the phase speed of the wave matches the local velocity of the background flow) in differentially rotating bodies (see Astoul et al. 2021, for an analytical characterization of this phenomenon).

Inertial wave propagation also depends strongly on the geometry of the container. While in the special case of a uniformly rotating homogeneous full sphere (or spheroid), the governing Poincaré equation for inertial oscillations is separable in the coordinates and allows for a complete set of regular solutions (i.e. normal modes, Bryan 1889; Ogilvie & Lin 2004; Wu 2005a, b; Lee 2020), the problem is in general ill-posed in other geometries due to the conflicting boundary conditions (Rieutord, Georgeot & Valdettaro 2000, 2001). In uniformly rotating spherical shells (relevant for models of the convective envelopes of stars and planets), singular modes show up in the inviscid limit, taking the form of limit cycles, which are periodic and have straight paths of characteristics (also called attractors), reflecting from the shell boundaries. In the presence of viscosity, these singular solutions are regularized and both their energy and viscous dissipation are focused along attractors, developing oscillating shear layers, which are centred around them with a certain width depending on the viscosity (Kerswell 1995; Rieutord & Valdettaro 1997, 2018). The activation of these shear layers is determined by the size of the shell compared to the core and by the forcing frequency. They are often associated with intense peaks of tidal dissipation when integrating over the shell (e.g. Ogilvie 2009), and these resonant peaks could be associated with the presence of hidden global modes, by analogy with those observed in full sphere geometry (Lin & Ogilvie 2021).

On top of the effects listed above, non-linear fluid effects should start to play a role in (at least) the most compact exoplanetary systems. For instance, the non-linear interactions between the large-scale equilibrium tide and smaller-scale inertial waves gives rise to the elliptical instability. This is thought to be a highly dissipative mechanism able to synchronize the planet and damp its obliquity (but not excite it) and eccentricity, but probably only for the very-shortest Hot Jupiters, and which has been extensively studied in experiments and numerical simulations (e.g. Kerswell 2002; Le Bars et al. 2010; Barker & Lithwick 2013, 2014; Cébron et al. 2013; Barker 2016). Interestingly, non-linear effects are likely to affect small-scale tidal waves for even smaller tidal amplitudes than for the large-scale equilibrium tidal flows, as discussed for gravity modes in e.g. Barker & Ogilvie (2010), Barker (2011), and see also the astrophysical discussion of Section 4 later in this paper. In his pioneering work, Tilgner (2007) demonstrated in numerical computations that the non-linear self-interaction of inertial waves in shear layers in spherical shells can generate strong axisymmetric zonal (or geostrophic) flows in spherical shells. Along these lines, Morize et al. (2010) and Favier et al. (2014, hereafter Paper I) verified through experiments and numerical simulations, respectively, that tidal forcing is able to trigger this non-linear phenomenon, through the non-uniform deposition of angular momentum inside a convective tidally deformed sphere or shell, akin to the convective regions of stars and giant planets. In such a previously poorly studied and complicated differentially rotating environment, the dissipative properties of tidal inertial waves can be deeply affected (as was also shown in the linear study of Baruteau & Rieutord 2013, with cylindrical differential rotation, which depends only on the distance from the axis of rotation). Moreover, approaching astrophysically relevant regimes for compact exoplanetary systems, which have low viscosities and high tidal forcing amplitudes, favours such non-linear effects and is also conducive to the occurrence of various kind of fluid instabilities, like parametric or shear instabilities (Jouve & Ogilvie 2014, Paper I).

In this paper, we revisit the effects of non-linearities on tidally forced inertial waves in direct numerical simulations of a convective shell, modelling the convective envelopes of a low-mass star or a giant gaseous planet subjected to an imposed tidal potential. In this study, we continue along the lines of Paper I by simulating an incompressible shell, relegating a study of realistic density profiles to future work. In contrast to previous numerical studies, we use a more realistic tidal forcing to excite inertial waves through the residual action of the equilibrium tide acting as an effective body force (e.g. Ogilvie 2013). In this way, it is possible to analytically and numerically unveil the role of the non-linear terms in redistributing energy between the different scales and components of the tidal flows, uncovering the origin of some unexpected evolution of angular momentum observed for some tidal frequencies in Paper I. In addition, more robust predictions (in the astrophysical context) can be made for the impact of non-linearities on tidal dissipation rates and on scaling laws for the emerging zonal flows.

In Section 2, we describe the analytical model governing the excitation, propagation, and dissipation of tidal inertial waves in an adiabatic and incompressible convective shell. We also derive the energy and angular momentum balances of tidal flows in this framework. Then, in Section 3, we perform and describe new results obtained from direct numerical simulations of tidal waves in spherical shells, comparing our results with the prior study of Paper I. Therein, we also identify the formation and effects of differential rotation, including corotation resonances, and the occurrence of parametric instabilities of inertial waves. Finally, we evaluate the relevance of our work for exoplanetary systems in Section 4, and we discuss the limitations of our model and conclude in Section 5.

2 MODEL FOR TIDALLY FORCED INERTIAL WAVES

2.1 Tidal flow decomposition and governing equations

We consider the convective shell of an undistorted spherical body, an idealized representation of the convective envelope of a low-mass star or a giant gaseous planet, which is subjected to the tidal potential Ψ due to an orbiting companion. The thickness of the convective shell is determined by the radial aspect ratio (the ratio of the inner to outer radii) which is set to α = 0.5 in the rest of this study. This choice is directly relevant for the envelopes of certain low-mass stars and giant planets with extended dilute cores (e.g. as observed for Jupiter), and is made so that the peculiar inertial wave behaviour in shells is captured without adopting a very thin shell. The body is assumed to have an initial uniform rotation |$\boldsymbol \Omega =\Omega \boldsymbol e_z$| along the vertical axis with unit vector |$\boldsymbol e_z$|⁠, with constant Ω that is small compared to the critical angular velocity of the body in order to neglect centrifugal effects. Inertial waves, which are restored by the Coriolis pseudo-force, can thus be excited between the cut-off frequencies ±2Ω in the fluid frame. For their treatment, the convection inside the shell is treated only in so far as it is assumed to have led to an adiabatic stratification (and it possibly provides a turbulent viscosity), which is, for instance, a reasonable assumption for most of the solar convective envelope (e.g. Christensen-Dalsgaard 2021), and probably also for the outermost convective envelope of Jupiter and Saturn (Debras & Chabrier 2019; Militzer, Wahl & Hubbard 2019). Indeed, the squared Brunt-Väisälä frequency² is in reality probably negative in the bulk of the convective envelope, thus driving convective motions, but is mostly negligible in absolute value compared to Ω², except possibly in a small region near the surface where the density is very low and the convection is no longer efficient. The convective shell is also considered incompressible as a first approach, with a uniform density ρ.

In our model, tidal flows are treated as Eulerian perturbations of the hydrostatic equilibrium of the body using the non-wavelike/wavelike decomposition (equivalent to the equilibrium/dynamical tide³) introduced by Ogilvie (2013). Accordingly, we use in the rest of the paper the subscripts _w and _nw to designate the tidal response related to the wavelike (dynamical) and non-wavelike (equilibrium) tides, respectively. In the frame corotating with angular frequency |$\boldsymbol{\Omega }$| and spherical polar coordinates (r, θ, φ) centred on the body at time t, the tidal potential can be written as

$$\begin{eqnarray} \Psi (r,\theta ,\varphi ,t)=A\, \mathrm{Re}\left[\left(\frac{r}{R}\right)^lY^m_l(\theta ,\varphi)\mathrm{e}^{-i\omega t}\right], \end{eqnarray}$$

(1)

where we introduce the tidal amplitude A = GM₂R²/a³, which depends on the mass M₂ and semimajor axis a of the perturber (with G the gravitational constant), R the radius of the tidally perturbed body, |$Y^m_l$| the orthonormalized spherical harmonic of degree l and order m (the azimuthal wavenumber) with m ≤ l, and ω the tidal forcing frequency in the rotating frame. Note that we only consider the dominant quadrupolar tidal component of Ψ for l = m = 2. This approximation is valid for quasi-coplanar and quasi-circular systems (i.e. for asynchronous tides, e.g. Ogilvie 2014), but is also the dominant component for an eccentric orbit. The tidal potential satisfies Laplace’s equation (ΔΨ = 0) inside the perturbed body whose tidally modified gravitational potential is Φ_nw.

As its name suggests, the non-wavelike tide is the non-oscillatory, instantaneous and hydrostatic, fluid response to the tidal potential, and its associated flow. By defining its displacement |$\boldsymbol \xi _\mathrm{nw}$|⁠, flow velocity |$\boldsymbol u_\mathrm{nw}=\partial _t\boldsymbol \xi _\mathrm{nw}$|⁠, pressure perturbation p_nw, and density perturbation ρ_nw, the non-wavelike tide thus satisfies the closed set of linearized equations:

$$\begin{eqnarray} \begin{aligned} &\partial _t\boldsymbol u_\mathrm{nw}=-\nabla \left(\frac{p_\mathrm{nw}}{\rho }+\Phi _\mathrm{nw}+\Psi \right),\ \ & p_\mathrm{nw}=-\rho \left(\Phi _\mathrm{nw}+\Psi \right),\\ &\rho _\mathrm{nw}=-\rho \boldsymbol \nabla \cdot \boldsymbol \xi _\mathrm{nw}, &\Delta \Phi _\mathrm{nw}=4\pi G\rho _\mathrm{nw}, \end{aligned} \end{eqnarray}$$

(2)

where the equations correspond respectively to the momentum, hydrostatic equilibrium, continuity, and Poisson equations. These are completed by the inner and outer boundary conditions:

$$\begin{eqnarray} \left\lbrace \begin{aligned} &\boldsymbol \xi _\mathrm{nw}\cdot \boldsymbol n=0 &\text{at}~~ r=\alpha R,\\ &\boldsymbol \xi _\mathrm{nw}\cdot \boldsymbol n=-\frac{\Phi _\mathrm{nw}+\Psi }{g} &\text{at}~~ r=R, \end{aligned}\right. \end{eqnarray}$$

(3)

with |$\boldsymbol n$| the unit vector normal to the boundaries and g the constant surface gravitational acceleration of the body. Note that here we adopt an impermeable inner boundary, which more appropriately describes tides in the envelopes of giant planets containing a solid core than in a low-mass star, where we must match the non-wavelike tidal displacement to the equilibrium tide in the radiative core (e.g. Barker 2020).

From equations (2) and (3), the non-wavelike displacement can be derived from a potential such that |$\boldsymbol \xi _\mathrm{nw}=-\boldsymbol \nabla X$| (Ogilvie 2013) where

$$\begin{eqnarray} X(r,\theta ,\varphi ,t)=\frac{C_\mathrm{t}R^2}{2(1-\alpha ^5)}\left[\left(\frac{r}{R}\right)^2+\frac{2}{3}\alpha ^5\left(\frac{R}{r}\right)^3\right]\mathrm{e}^{-i\omega t}Y_2^2(\theta ,\varphi), \\ \end{eqnarray}$$

(4)

with C_t a constant embodying the amplitude of the tidal forcing. The dimensionless parameter C_t can be expressed as

$$\begin{eqnarray} C_\mathrm{t}=(1+k_2)\epsilon , \end{eqnarray}$$

(5)

introducing two useful astrophysical quantities: the usual dimensionless tidal amplitude parameter ϵ = (R/a)³M₂/M₁ with M₁ the mass of the perturbed body (Ogilvie 2014; Barker 2016), along with the real part of the quadrupolar Love number |$k_2=\frac{3\rho }{5\overline{\rho }-3\rho }$| with |$\overline{\rho }$| the mean density in the whole body (core plus envelope). The quantity ϵ quantifies the ratio of the amplitude of the tidal potential over the self-gravity of the perturbed body at its surface, while k₂ accounts for its non-dissipative tidal deformation and is determined by calculating the real part (in-phase component) of the ratio of the amplitude of the perturbed gravitational potential to the tidal potential, Φ_nw/Ψ, for l = m = 2 (Love 1911; Ogilvie 2013). Since k₂ is constrained such that 0 < k₂ ≤ 3/2, with its maximum reached for a full homogeneous sphere (when |$\rho =\overline{\rho }$|⁠), the value of C_t using α = 0.5 is also a measure of the magnitude of the tidal forcing with a similar value to ϵ (their astrophysical values will be discussed further in Section 4).

With the non-wavelike perturbation now being defined as above, the remaining vortical action of the rotation on this tide |$\boldsymbol f_\mathrm{t}= -2\boldsymbol \Omega \wedge \boldsymbol u_\mathrm{nw}$| is used to excite the non-zero frequency (i.e. wavelike) tidal response (Ogilvie & Lin 2004; Ogilvie 2005, 2013). The non-linear equation of motion for the wavelike (dynamical) tides in the fluid frame is thus:

$$\begin{eqnarray} \partial _t\boldsymbol u_\mathrm{w}+2\boldsymbol \Omega \wedge \boldsymbol u_\mathrm{w}+(\boldsymbol u\cdot \boldsymbol \nabla)\boldsymbol u =-\frac{\boldsymbol \nabla p_\mathrm{w}}{\rho }+\nu \Delta \boldsymbol u_\mathrm{w} +\boldsymbol f_\mathrm{t}, \end{eqnarray}$$

(6)

where |$\boldsymbol u=\boldsymbol u_\mathrm{w}+\boldsymbol u_\mathrm{nw}$| is the total perturbed flow and ν is an effective viscosity accounting for the sum of the molecular kinematic viscosity and an eddy viscosity modelling the action of turbulent convective motions on tidal flows. We point out that the molecular viscosity in convective regions may be negligible compared to the turbulent one. While realistic values of the kinematic Ekman number E = ν/(ΩR²) usually range between⁴ E = 10⁻¹⁰ and 10⁻²⁰ in convective envelopes, the turbulent one is often taken to be 10⁻⁴ − 10⁻⁶ when calculated from mixing-length theory in stellar evolution models (e.g. as in Astoul et al. 2019; Fournier, Gizon & Hyest 2022). It should also be stressed that while the action of turbulent convective motions on the equilibrium tides has long been modelled as a turbulent eddy viscosity,⁵ whether or not this is fully justified, the modelling of the action of turbulent motions on the wavelike tides is still (even more) highly speculative.

In equation (6), we neglect the dissipative term associated with the non-wavelike tidal flow (⁠|$\nu \Delta \boldsymbol u_\mathrm{nw}$|⁠) since this term is estimated to be small relative to those retained when inertial waves are excited (see also the last paragraph of Section 2.2). In addition, we assume that the non-wavelike tide is perfectly maintained on the time-scale of our simulations, which are much shorter than tidal evolutionary time-scales and thus probe the ‘instantaneous’ tidal response for a given system (similar reasoning has been used in Barker & Astoul 2021). However, it is worth noting that we have so far kept the non-wavelike tidal flow in the non-linear advection term |$(\boldsymbol u\cdot \boldsymbol \nabla)\boldsymbol u$|⁠, which thus consists of four terms, and we will discuss this choice in Section 3.2. The equation of motion equation (6) is also combined with the continuity equation for the wavelike perturbations |$\boldsymbol \nabla \cdot \boldsymbol u_\mathrm{w}=0$|⁠, assuming incompressibility. To complete this model, we adopt stress-free and impenetrable boundary conditions for the wavelike flow at both the inner r = αR and outer r = R boundaries. Namely, no tangential stress (⁠|$\boldsymbol n\wedge ([\sigma ]\boldsymbol n)=\boldsymbol 0$| with [σ] the viscous stress tensor) and no radial velocity |$\boldsymbol u_\mathrm{w}\cdot \boldsymbol n=0$| at the boundaries (e.g. Rieutord 2015). Though a no-slip boundary condition might be preferred at the inner boundary in giant gaseous planets if they possess a solid core, this choice of boundary condition does not seem to significantly affect the propagation of inertial modes in a spherical shell (Rieutord & Valdettaro 2010, Paper I), and our choice is also more appropriate for modelling stars.

We non-dimensionalize our simulations using the length-scale R, the time-scale Ω⁻¹, and the density value ρ.

2.2 Kinetic and angular momentum balances

We derive in this section the energy and angular momentum balances in our model, and analyse the energetic redistribution and exchanges between the wavelike and non-wavelike tidal flows, which are of prime importance for analysing our non-linear simulations (Section 3), and motivate our choices for neglecting certain non-linear terms driving unrealistic fluxes through the boundary (Section 3.2).

From the scalar product between |$\rho \boldsymbol u_\mathrm{w}$| and the momentum equation (6) and after spatial integration 〈 · 〉 over the shell’s volume V, one can get the kinetic energy |$K=\langle \rho \boldsymbol u_\mathrm{w}^2/2\rangle$| balance for tidal inertial waves:

$$\begin{eqnarray} \partial _t K= I_\mathrm{nw-w}+I_\mathrm{w-nw}-D_\nu +P_\mathrm{t}, \end{eqnarray}$$

(7)

where we have introduced the non-linear energy transfer terms |$I_\mathrm{i-j}=-\langle \rho \boldsymbol u_\mathrm{w}\cdot (\boldsymbol u_i\cdot \boldsymbol \nabla)\boldsymbol u_j\rangle$| with i, j ∈ {w, nw} between the wavelike and non-wavelike flows, the tidal power

$$\begin{eqnarray} P_\mathrm{t}=\langle \rho \boldsymbol u_\mathrm{w}\cdot \boldsymbol f_\mathrm{t}\rangle , \end{eqnarray}$$

(8)

reflecting the energy injected into the wavelike flows by the tidal work, which is damped by the tidal viscous dissipation

$$\begin{eqnarray} D_\nu =-\langle \rho \nu \boldsymbol u_\mathrm{w}\cdot \Delta \boldsymbol u_\mathrm{w}\rangle , \end{eqnarray}$$

(9)

of the wavelike flow.

In the derivation of equation (7), several terms have vanished due to the assumptions of our model. The energy transfer term involving the wavelike/wavelike non-linearity does not appear in the energy balance because:

$$\begin{eqnarray} I_\mathrm{w-w}=-\langle \rho \boldsymbol u_\mathrm{w}\cdot (\boldsymbol u_\mathrm{w}\cdot \boldsymbol \nabla)\boldsymbol u_\mathrm{w}\rangle = -\frac{\rho }{2}\int _{\delta V}\boldsymbol u_\mathrm{w}^2~\boldsymbol u_\mathrm{w}\cdot \boldsymbol n\, \mathrm{d}S=0, \end{eqnarray}$$

(10)

when applying the divergence theorem with a surface element dS and a unit vector |$\boldsymbol n$| normal to the shell boundaries, together with incompressibility and impenetrability at the boundaries δV for the wavelike flow (⁠|$\boldsymbol \nabla \cdot \boldsymbol u_\mathrm{w}=0$| and |$\boldsymbol u_\mathrm{w}\cdot \boldsymbol n=0$|⁠, respectively). For the same reasons, there is no term associated with the wavelike pressure. The transfer term associated with the non-wavelike/non-wavelike non-linear advection also cancels out:

$$\begin{eqnarray} I_\mathrm{nw-nw}&=&-\langle \rho \boldsymbol u_\mathrm{w}\cdot (\boldsymbol u_\mathrm{nw}\cdot \boldsymbol \nabla)\boldsymbol u_\mathrm{nw}\rangle \\ &=&-\left\langle \rho \boldsymbol u_\mathrm{w}\cdot \left[\boldsymbol \nabla (\boldsymbol u_\mathrm{nw}^2/2)-\boldsymbol u_\mathrm{nw}\wedge (\boldsymbol \nabla \wedge \boldsymbol u_\mathrm{nw})\right]\right\rangle \\ &=&-\frac{\rho }{2}\int _{\delta V} \boldsymbol u_\mathrm{nw}^2~\boldsymbol u_\mathrm{w}\cdot \boldsymbol n\, \mathrm{d}S=0, \end{eqnarray}$$

(11)

where we use in the second line standard vector calculus identities as in Barker & Astoul [2021, equation (7)], along with the fact that the non-wavelike flow is curl-free and the wavelike flow is impenetrable at the boundaries.

By using the divergence-free assumption for the non-wavelike flow |$\boldsymbol \nabla \cdot \boldsymbol u_\mathrm{nw}=0$|⁠, the first ‘mixed’ energy transfer term I_nw-w in equation (7) can be recast as

$$\begin{eqnarray} I_\mathrm{nw-w}=-\langle \rho \boldsymbol u_\mathrm{w}\cdot (\boldsymbol u_\mathrm{nw}\cdot \boldsymbol \nabla)\boldsymbol u_\mathrm{w}\rangle =-\frac{\rho }{2}\int _{\delta V}~\boldsymbol u_\mathrm{w}^2\boldsymbol u_\mathrm{nw}\cdot \boldsymbol n\, \mathrm{d}S. \end{eqnarray}$$

(12)

In a tri-axial ellipsoidal body, which more realistically corresponds to a centrifugally and tidally deformed star or planet, the streamlines of the non-wavelike flow are supposed to be tangent to the surface boundary in the bulge frame (rotating with the orbit, for asynchronous tides). Thus, we expect that the surface integral in equation (12) also vanishes on average in such a deformed geometry in the bulge frame, since |$\boldsymbol u_\mathrm{nw}\cdot \boldsymbol n=0$| there (Barker, Braviner & Ogilvie 2016; Barker & Astoul 2021). In the fluid frame, we do not expect this term to vanish identically but we will explain later that its magnitude is expected to be much smaller than I_w-w (in a point-wise sense). This integral does not vanish in our undeformed spherical shell model due to advection of kinetic energy through the spherical boundaries by the non-wavelike flow. This has important consequences for the evolution of the angular momentum, as will be seen in the following. Finally, the last non-zero energy transfer term gives

$$\begin{eqnarray} I_\mathrm{w-nw}=-\langle \rho \boldsymbol u_\mathrm{w}\cdot (\boldsymbol u_\mathrm{w}\cdot \boldsymbol \nabla)\boldsymbol u_\mathrm{nw}\rangle =\langle \rho \boldsymbol u_\mathrm{nw}\cdot (\boldsymbol u_\mathrm{w}\cdot \boldsymbol \nabla)\boldsymbol u_\mathrm{w}\rangle , \end{eqnarray}$$

(13)

using the same arguments as for equation (10). The quantity I_w-nw hence features the Reynolds stresses involving correlations between wavelike flows components and gradients of the non-wavelike flow.

Note that, if we define |$\mathcal {K}_{i}=\rho |\boldsymbol {u}_i^2|/2$| where i ∈ {w, nw},

$$\begin{eqnarray} -\langle \rho \boldsymbol u_\mathrm{w}\cdot (\boldsymbol u\cdot \boldsymbol \nabla)\boldsymbol u\rangle &=&\sum _{i,j\in \lbrace \mathrm{w,nw}\rbrace }I_\mathrm{i,j} \\ &=&-\int _{\delta V}\mathcal {K}_\mathrm{w}\ \boldsymbol u_\mathrm{nw}\cdot \boldsymbol n\, \mathrm{d}S~+I_\mathrm{w-nw}, \end{eqnarray}$$

(14)

and similarly

$$\begin{eqnarray} -\langle \rho \boldsymbol u_\mathrm{nw}\cdot (\boldsymbol u\cdot \boldsymbol \nabla)\boldsymbol u\rangle &= -\int _{\delta V}(\mathcal {K}_\mathrm{nw}+\boldsymbol u_\mathrm{nw}\cdot \boldsymbol u_\mathrm{w})\ \boldsymbol u_\mathrm{nw}\cdot \boldsymbol n\, \mathrm{d}S-I_\mathrm{w-nw}. \\ \end{eqnarray}$$

(15)

From these two equations, I_w-nw is clearly identified as an energy transfer term between the wavelike and non-wavelike flows. It takes the opposite sign in the equation for the energy of the non-wavelike flow, if we were to allow this component to evolve. Equations (14) and (15) together imply

$$\begin{eqnarray} \begin{aligned} -\langle \rho \boldsymbol u\cdot (\boldsymbol u\cdot \boldsymbol \nabla)\boldsymbol u\rangle =&-\int _{\delta V}\mathcal {K}_\mathrm{tot}\ \boldsymbol u_\mathrm{nw}\cdot \boldsymbol n\, \mathrm{d}S, \end{aligned} \end{eqnarray}$$

(16)

where |$\mathcal {K}_\mathrm{tot}=\rho |\boldsymbol u_\mathrm{w}+\boldsymbol u_\mathrm{nw}|^2/2$|⁠. Note also that

$$\begin{eqnarray} \langle \rho \boldsymbol {u}_\mathrm{w} \cdot \boldsymbol {u}_{\mathrm{nw}}\rangle =-\langle \rho \boldsymbol {u}_\mathrm{w} \cdot \nabla \dot{X}\rangle =-\int _{\delta V}\rho \dot{X}\boldsymbol {u}_\mathrm{w} \cdot \boldsymbol n\, \mathrm{d}S=0, \end{eqnarray}$$

(17)

upon applying the divergence theorem and the boundary conditions, implying that

$$\begin{eqnarray} K_{\mathrm{tot}}=\langle \mathcal {K}_\mathrm{tot}\rangle =\langle \rho |\boldsymbol u_\mathrm{w}+\boldsymbol u_\mathrm{nw}|^2\rangle /2=\langle \mathcal {K}_\mathrm{w}+\mathcal {K}_{\mathrm{nw}}\rangle =K + K_{\mathrm{nw}}, \end{eqnarray}$$

(18)

where |$K=\langle \mathcal {K}_\mathrm{w}\rangle$| (as above) and |$K_\mathrm{nw}=\langle \mathcal {K}_{\mathrm{nw}}\rangle$|⁠.

If we were to impose a radial velocity u_r at the outer boundary, like in Ogilvie (2009) or Paper I, instead of a tidal body force, the (non-linear) energy balance would be:

$$\begin{eqnarray} \partial _t K_\mathrm{tot}=-\int _{\delta V(r=R)}(\mathcal {K}_\mathrm{tot}+p_\mathrm{w})\ u_r\, \mathrm{d}S~~~-D_\nu , \end{eqnarray}$$

(19)

where the first term inside the integral is omitted in linear calculations (since it comes from non-linear terms), and the second term represents tidal forcing in their model.

The total angular momentum of the fluid in the shell is

$$\begin{eqnarray} \boldsymbol L=\int _V \rho \boldsymbol r\wedge \boldsymbol u\, \mathrm{d}V, \end{eqnarray}$$

(20)

where |$\boldsymbol r$| is the position vector. By taking the vector product between |$\boldsymbol r$| and the momentum equation (6) and integrating over the volume, we get the angular momentum balance

$$\begin{eqnarray} \partial _t\boldsymbol L=-2\int _V\rho \boldsymbol r\wedge (\boldsymbol \Omega \wedge \boldsymbol u)\, \mathrm{d}V-\int _V \rho \boldsymbol r\wedge (\boldsymbol u\cdot \boldsymbol \nabla \boldsymbol u)\, \mathrm{d}V, \end{eqnarray}$$

(21)

where the pressure torque vanishes, as well as the viscous torque by the adoption of stress-free boundary conditions (Paper I; with the difference that |$\boldsymbol u=\boldsymbol u_\mathrm{w}+\boldsymbol u_\mathrm{nw}$| here Jones et al. 2011). By injecting the new variable |$\tilde{\boldsymbol u}=\boldsymbol u+\boldsymbol \Omega \wedge \boldsymbol r$| similarly as in Paper I, equation (21) yields

$$\begin{eqnarray} \partial _t\boldsymbol L=-\int _V\boldsymbol \rho r\wedge \left[(\tilde{\boldsymbol u}\cdot \boldsymbol \nabla)\tilde{\boldsymbol u}\right]\, \mathrm{d}V=\int _{\partial V}\rho (\boldsymbol r\wedge \tilde{\boldsymbol u})\boldsymbol u_\mathrm{nw}\cdot \boldsymbol n\, \mathrm{d}S\ne \boldsymbol {0}, \\ \end{eqnarray}$$

(22)

where the incompressibility of the total flow |$\boldsymbol u$| has been used in both equalities along with the impenetrability condition for |$\boldsymbol u_\mathrm{w}$|⁠. Therefore, the total angular momentum is not conserved over time in our spherical model, and its evolution is entirely controlled by the surface flux sustained by the ‘unrealistic’ radial non-wavelike flow through the boundaries in our spherical model (since we look at short intervals of time relative to tidal evolutionary time-scales, and thus assume |$\boldsymbol u_\mathrm{nw}$| to be perfectly maintained). As we will show partly below and in the next section, this non-vanishing flux term, due to either imposed radial velocity boundary conditions, as in Paper I [see equation (19)], or through the inclusion of the non-linearity terms involving the non-wavelike tide in our spherical model, is important and causes unexpected tidal evolution, such as the sustained unrealistic tidal de-synchronization observed in the aforementioned study. Note that the tidal effective forcing itself does not contribute to changing the total angular momentum (in the context of our model that focuses on an ‘instantaneous snapshot’ of the system’s evolution) since |$\langle \boldsymbol {r}\wedge \boldsymbol {f}_\mathrm{t} \rangle =\boldsymbol {0}$| (due to the integral over azimuthal angles).

We can crudely estimate the importance of I_nw-w compared to the other terms in the energy balance in equation (7) using physical scaling arguments. The non-wavelike tide has a typical velocity magnitude u_nw ∼ RωC_t and length-scale L ∼ R, because it represents large-scale flows in the convective envelope, its time-scale is the inverse of the tidal frequency ω, and its displacement scales as ξ_nw ∼ RC_t (equation (4)). In a spherical shell, the scales of the wavelike tide are controlled by the structure of the inertial wave shear layers, whose thicknesses depend on the Ekman number (in linear theory). The wavelike tide has corresponding velocity u_w ∼ E^−βRωC_t and a length-scale ℓ ∼ E^αR, where α, β > 0 are coefficients depending on the particular (singular) mode in the shell: for shear layers α = 1/3 (if ℓ is the wavelength of the oscillation inside the shear layer) or α = 1/4 (if ℓ is the width of the shear layer) and β = 1/6; for Ekman boundary layers near the critical latitude⁶ α = 2/5 and β = 1/5 (these scaling have been analytically and numerically verified for instance in Kerswell 1995; Rieutord & Valdettaro 2018, and Paper I). Given the low values of the turbulent or microscopic Ekman number E ≪ 1, in any of these cases we have ℓ ≪ L and u_w ≫ u_nw in spherical shells. Furthermore, we also have the following point-wise approximations

$$\begin{eqnarray} \mathcal {I}_\mathrm{nw-w}&\sim \rho u_\mathrm{w}^2u_\mathrm{nw}/\ell , & \mathcal {I}_\mathrm{w-nw}\sim \rho u_\mathrm{w}^2u_\mathrm{nw}/L, \\ \mathcal {I}_\mathrm{w-w}&\sim \rho u_\mathrm{w}^3/\ell , & \mathcal {I}_\mathrm{nw-nw}\sim \rho u_\mathrm{w}u_\mathrm{nw}^2/L, \\ \mathcal {D}_\nu &\sim \rho \nu u_\mathrm{w}^2/\ell ^2, & \mathcal {P}_\mathrm{t}\sim \Omega u_\mathrm{w} u_\mathrm{nw}, \end{eqnarray}$$

where |$\mathcal {I}_\mathrm{i-j}$| and |$\mathcal {P}_\mathrm{t}$| are the energy transfer terms (and |$\mathcal {D}_\nu$| the viscous term) before computing the volume integral. Hence, we typically expect the local ordering

$$\begin{eqnarray} \begin{aligned} ~~\mathcal {I}_\mathrm{w-w} &\gg \mathcal {I}_\mathrm{nw-w} \gg \mathcal {I}_\mathrm{w-nw}\gg \mathcal {I}_\mathrm{nw-nw}, \end{aligned} \end{eqnarray}$$

(23)

in magnitude for shear layers and critical latitudes, for the typical Ekman numbers and tidal amplitudes in stars and planets. In addition, while these orderings are also expected to hold in our simulations for the typical values of the Ekman number E ∼ 10⁻⁴ − 10⁻⁶ and tidal forcing parameters C_t ∼ 10⁻² − 10⁻¹ we use in our simulations, we may not expect |$\mathcal {I}_\mathrm{nw-w}\ll \mathcal {I}_\mathrm{w-w}$|⁠. By consequence, if no cancellations in the surface integral of I_nw-w occur, it is expected that this term will play a prominent role in the energy and angular momentum balances equations (7) and (22), and could drive unrealistic evolution. This provides a first reason to motivate our neglect of the mixed (non-wavelike/wavelike and wavelike/non-wavelike) non-linearities over the wavelike/wavelike one in simulations. If we neglect this term, we emphasize that the total angular momentum is conserved over time |$\partial _t\boldsymbol L=0$|⁠, and the energy balance now reads:

$$\begin{eqnarray} \partial _t K=-D_\nu +P_\mathrm{t}. \end{eqnarray}$$

(24)

Morover, given the integral and local relationships equations (11) and (23), the non-wavelike/non-wavelike non-linearity contributes negligibly to the energy exchanges in the body. It only modifies the pressure, which is constrained by incompressibility, and therefore this term has no effect in our model with our adopted boundary conditions. Finally, assuming the same magnitude of the viscosity holds for both the wavelike and non-wavelike flows, our neglect of the dissipative term for the latter is further justified since νΔu_nw/νΔu_w ∼ E^{2α + β} ≪ 1.

3 NON-LINEAR SIMULATIONS

We now describe and analyse non-linear simulations of tidally forced inertial waves governed by equation (6) for an incompressible flow satisfying stress-free boundary conditions in a convective shell, performed by the open source code magic.

3.1 Numerical model

We use version 5.10 of the pseudo-spectral code magic⁷ designed for 3D (magneto-)hydrodynamical simulations of fluid motions in a spherical shell or full sphere (e.g. Christensen et al. 2001; Wicht 2002, for Boussinesq implementation and benchmarks). We employ Chebyshev polynomials in the radial direction and a spherical harmonic decomposition in the azimuthal and latitudinal directions. The code supports a number of mixed implicit/explicit time-stepping schemes, where the non-linear terms and the Coriolis force are treated explicitly in real (spatial) space to avoid inverting a large matrix (resulting from non-linear mode coupling), while the remaining linear terms are treated implicitly in spectral space (hence the name pseudo-spectral). The chosen time-stepping is based on an IMEX multistep method, specifically choosing the canonical combination of second-order Crank-Nicolson and Adams-Bashforth schemes to treat the implicit and explicit terms, respectively. The choice of Chebyshev polynomials evaluated at the Gauss-Lobatto collocation gridpoints, instead of low-order finite differences in radius, for example, allows for much better accuracy for smooth solutions for the same number of points, and it is also well adapted to resolve thin (viscous) layers near the boundaries, thanks to a denser grid of points close to the inner and outer radii. All of our simulations have benefited from the MKL and SHTns (Schaeffer 2013) libraries for fast Fourier, Legendre, and spherical harmonics transforms between real and spectral space, and vice versa, along with MPI parallelization on high-performance computing clusters, which significantly improves the speed of our calculations. Like in the code PARODY used in Paper I, the wavelike velocity is decomposed using poloidal–toroidal potentials, which automatically ensure that the flow is a solenoidal (i.e. divergence-free) field. Details about this decomposition and the spherical harmonic and Chebyshev representations can be found in the manual of the magic code, at https://magic-sph.github.io/numerics.html#secnumerics. Spectral projections of the unforced momentum equation for the poloidal and toroidal potentials, and the associated stress-free boundary conditions, in this formalism, are also described there.

We have modified magic subroutines to implement in real space the tidal forcing |$\boldsymbol f_\mathrm{t}$| and each of the non-linear terms involving the non-wavelike flow which contributes to the energy exchanges according to Section 2.2, namely the mixed non-linear terms |$(\boldsymbol u_\mathrm{nw}\cdot \boldsymbol \nabla)\boldsymbol u_\mathrm{w}$| and |$(\boldsymbol u_\mathrm{w}\cdot \boldsymbol \nabla)\boldsymbol u_\mathrm{nw}$|⁠. Since these changes are highly non-trivial (particularly for the mixed non-linear terms using the code variables), we ensured that we had done this correctly by comparing results obtained independently with an implementation for all of these terms (and subsets of these terms) computed using the spectral element code Nek5000 (which we also used for some of the simulations in Paper I). This detailed comparison is omitted here to save space. We have also modified the viscous dissipation routines to calculate the appropriate viscous dissipation and tidal power, and implemented routines to enable volume integrals (of products) of real-space quantities to be computed, so that, for example, we can evaluate each of the terms that contribute to equation (7).

In most of our simulations,⁸ we adopted a maximum spherical harmonic degree of l_max = 85 (meaning n_φ = 256 longitudinal grid points and n_θ = 128 latitudinal grid points), along with a number of radial grid points n_r = 97 (other spatial resolution choices will be specified later). This spatial resolution is sufficient to allow more than three orders of magnitude difference in the kinetic energy spectrum between the low- and high-order harmonic degrees and azimuthal wavenumbers, while maintaining a fairly good speed for the simulations.⁹ We also used a time-step of δt = 10⁻² in most of our simulations, and checked that this value is low enough to ensure both accuracy (by checking results for convergence using smaller δt) and stability according to the Courant–Friedrichs–Lewy condition. Lastly, we checked that strictly enforcing conservation of the total (equatorial and axial) angular momentum at every time-step, which is an option in magic, does not change our results with only wavelike-wavelike non-linearities.

3.2 Surface versus body forcing

The main difference in our model with Paper I, apart from the use of a different numerical code, is the way inertial waves are forced. They imposed an outer radial velocity at the surface r = R to model the non-wavelike tidal flow and how it forces inertial waves, which we will refer to as ‘surface forcing’ in the following. In that case, the impenetrability condition only applies at the inner radius r = αR. When studying the linear tidal response, surface forcing and the body forcing that we use are equivalent in the asymptotic limit of small tidal forcing frequencies compared to the dynamical frequency |$\omega _\mathrm{d}=\sqrt{g/R}$|⁠, and they are approximately similar for astrophysically relevant small but non-zero values of ω/ω_d. This regime is valid for slowly rotating stars or planets for which centrifugal effects can be neglected (e.g. Ogilvie 2009). Within this limit, both forcing mechanisms can be related by equating the amplitude of the non-wavelike velocity in both models, i.e.

$$\begin{eqnarray} \sqrt{\frac{32\pi }{15}}A\approx \omega C_\mathrm{t}, \end{eqnarray}$$

(25)

where A is the forcing amplitude in Paper I and C_t was defined above. For the purposes of recovering the results of Paper I, our tidal forcing amplitude C_t must be scaled with frequency, since they chose to fix A instead.

We now discuss non-linear simulations performed for three initial tidal forcing frequencies ω = 1.05, 1.1, and 1.15, which were those initially analysed in the linear regime by Ogilvie (2009) and in the prior non-linear simulations of Paper I. These three cases were chosen because they exhibited drastically different behaviours of the dissipation rate as a function of Ekman number, and had correspondingly different flow structures, despite having similar tidal frequencies. The tidal forcing amplitude C_t has been chosen to satisfy equation (25) with |$A=10^{-2}/\sqrt{32\pi /15}\approx 0.00386$| (e.g. for comparison with Fig. 6 of Paper I), but which is slightly different in our simulations for the three tidal frequencies (but C_t ≈ 0.009). Fig. 1 shows each of the terms appearing in the energy balance equation (7) over time, for the three non-linear simulations with different frequencies and E = 10⁻⁵.

Figure 1.

Contributions to the energy balance in equation (7) as a function of time for three different forcing frequencies ω, all with E = 10⁻⁵ and C_t = 10⁻²/ω. A running average with a period 2π/ω has been performed to smooth the curves. All y-axis quantities have been re-scaled by a factor A² × 32π/15 ≈ 10⁻⁴ by taking A = 0.00386 so as to provide a comparison with Fig. 6 of Paper I.

First, we checked that the energy balance as described in equation (7) is well satisfied in these simulations. Moreover, in all three cases, the energy transfer term I_w-nw is negligible compared to the viscous dissipation and tidal power, while the term I_nw-w sustaining the unphysical non-wavelike flux through the boundaries contributes significantly to the energy balance for ω = 1.05 and 1.1. For the ω = 1.15 case, the energy injected by the tidal forcing is mainly viscously dissipated, without important transfers of energy between wavelike and non-wavelike flows. These observations are consistent with our discussion of how the terms in the energy balance scale in the last paragraph of Section 2.2, by noting that highly dissipative shear layers are excited in linear simulations for ω = 1.05 and 1.1 unlike ω = 1.15 (as described in Ogilvie 2009), making the scaling laws most relevant for the first two frequencies there. However, one may wonder whether the ‘physical’ transfer term I_w-nw is always negligible for our parameters (as seen in Fig. 1) or not. In Fig. 2, we investigate how the omission of I_w-nw impacts the energy balance for different tidal forcing amplitudes C_t and Ekman numbers E at ω = 1.1. In the left-hand panel, one can observe that this term is important for very high-tidal forcing amplitudes C_t ≳ 0.1 (the energy transfer has been checked beforehand), but for moderate-to-low C_t, typically C_t ≲ 0.05, the energy balance without I_w-nw remains almost unchanged for these parameters (similar results are found for ω = 1.05). Changes in Ekman number (right-hand panel), from E = 10⁻⁴ to E = 2 × 10⁻⁶, show that I_w-nw is also playing a negligible role in the energy balance for C_t ≈ 0.009.

Figure 2.

Comparison between the left-hand and right-hand sides of the energy balance equation (7) for ω = 1.1 without the transfer term I_w-nw involving wavelike correlations and gradients of the non-wavelike flow. Left: for different amplitudes of the tidal forcing C_t at fixed Ekman number E = 10⁻⁵ (not re-scaled). Right: for different values of the Ekman number E at fixed tidal forcing amplitude C_t ≈ 0.009 (re-scaled as in Fig. 1).

Since the non-wavelike flux through the boundaries described by I_nw-w in these simulations is artificial, we wish to remove the contribution from this term from our simulations. To do so, we can remove both of the non-linear terms involving the non-wavelike tide in equation (6), which is justified to represent the dominant balance of terms in the astrophysical regime based on the physical scaling arguments we outlined in Section 2.2. However, this approach is only justified in our simulations if I_w-nw also does not contribute significantly to the energy balance (since we also remove this term in our approach). From the above discussion, this is likely to be valid for C_t ≲ 0.05 and E ∼ 10⁻⁶ − 10⁻⁴. It is thus of interest to compare simulations which feature all of the non-linear terms with ones that include only the wavelike/wavelike non-linearity. This is done in Fig. 3, displaying the kinetic energy and dissipation over time for the same three tidal forcing frequencies as in Fig. 1. Linear predictions computed with an independent spherical spectral code (similar to e.g. Ogilvie 2009) are also displayed in both panels as the horizontal dashed lines for each frequency. When all non-linear terms are included, we recover the kinetic energy and dissipation shown in Fig. 6 of Paper I after appropriate re-scaling by a factor A² × 32π/15, except that our model does not include the non-wavelike tide dissipation |$\rho \nu \Delta \boldsymbol u_\mathrm{nw}$| which is included in their model. However, in all three simulations, this term has been checked to be at least one order of magnitude smaller than the term |$\rho \nu \Delta \boldsymbol u_\mathrm{w}$|⁠. Our quantitative agreement with Paper I (for the times shown in their Fig. 6) verifies that our implementation of tidal forcing and non-linear terms is correct.

Figure 3.

Kinetic energy K (left) and dissipation D_ν (right) over time for three different forcing frequencies (colour), comparing linear calculations (horizontal dotted lines) with simulations including wavelike/wavelike (w-w, solid lines) non-linearities, or all non-linear terms (all nl, dashed lines). Kinetic energy and dissipation have been re-scaled as in Fig. 1. The Ekman number is set to E = 10⁻⁵ and the forcing amplitude C_t = 10⁻²/ω.

The difference in kinetic energy and dissipation between linear and full non-linear simulations keeps increasing with time, without a steady state being reached at the end of the simulation (here for Ωt = 5000) when all non-linear terms are included, contrary to the case with only wavelike/wavelike non-linear terms in which K and D_ν are closer to the corresponding linear predictions. This is likely to be due to the evolution of the total angular momentum |$\boldsymbol L$| in the full non-linear case [primarily due to unphysical non-wavelike fluxes through the outer boundary, as demonstrated in equation (22)]. The vertical component of angular momentum is shown in the left-hand panel of Fig. 4, which shows that the angular momentum is increasing with time, as would be expected due to tidal spin-synchronization, while |$\boldsymbol L$| is conserved in wavelike/wavelike non-linear simulations (since we are not then modelling long-term tidal evolution, just ‘instantaneous transfers’ in a snapshot of the system).

Figure 4.

Left: vertical component of the total angular momentum in the rotating frame, L_z, as a function of time Ωt in simulations, including all non-linearities (w-w, nw-w, and w-nw) for three initial forcing frequencies ω. Right: non-linear dissipation D_ν (solid lines) as a function of normalized frequency ω/Ω* in simulations from Ωt = 20 to 5000, with all (bottom) or only wavelike/wavelike (top) non-linear terms (blue pale bullets indicate the final dissipation when an average steady state is reached). The colour of each line indicates the kinetic energy in the differential rotation E_dr. For comparison, the frequency-dependent linear dissipation for a uniformly rotating body is plotted as the dashed black lines and the impulsive energy injected |$\hat{E}$| as the horizontal dotted line. All the dissipation rates are re-scaled by a factor |$C_\mathrm{t}^2$| (contrary to previous figures showing dissipation where a different scaling is used to match with Paper I) and the bottom curves are shifted downwards by a factor 10⁻² for ease of visualization. The tidal power term computed from the final steady-state attained in linear simulations with a background cylindrical zonal flow are also added as blue stars.

Furthermore, the fact that the difference between linear and non-linear K and D_ν is more pronounced for ω = 1.05 and 1.1 compared to the ω = 1.15 case seems to be linked to the strong focusing of waves into localized shear layers for the two lowest frequencies, as described in Ogilvie (2009, as opposed to the more diffuse ‘space-filling’ distribution of inertial wave beams for ω = 1.15). While for the linear ω = 1.1 case, the excitation of a periodic wave attractor has been identified, the linear ω = 1.05 case rather features large-scale structures hidden beneath shear layers (see Lin & Ogilvie 2021, and possibly visible in the top-left panel of Fig. 5). These have been shown to be potentially related to global modes in full sphere geometry, explaining the resonant peak in linear dissipation rates¹⁰.

Figure 5.

Radial component of the wavelike velocity in the meridional plane for the same three frequencies as in Fig. 6. Top: in linear simulations computed with LSB. Bottom: in non-linear simulations computed with magic where |$\varphi \equiv \omega t/m\, [2\pi ]-\pi /2$| (the π/2 factor is necessary to find the correct phase for comparison). The extrema of the colourbars are chosen to be the same values as those in the above panels.

The non-linear self-interactions of inertial waves in such localized shear layers favour the emergence of strong zonal flows through the non-uniform deposition of angular momentum in the shell (Tilgner 2007, Paper I). These zonal flows can be identified by calculating the energy of the differential rotation

$$\begin{eqnarray} E_\mathrm{dr}=\frac{\rho }{2}\int \left[\langle u_\varphi \rangle _\varphi -\delta \Omega r\sin \theta \right]^2\, \mathrm{d}V, \end{eqnarray}$$

(26)

where 〈u_φ〉_φ is the azimuthally averaged azimuthal component of the flow, and

$$\begin{eqnarray} \delta \Omega =\frac{1}{V}\int \frac{u_\varphi }{r\sin \theta }\, \mathrm{d}V, \end{eqnarray}$$

(27)

is the newly created volume-averaged rotation rate of the fluid in the rotating frame. The latter is superimposed upon the mean rotation rate Ω such that the volume-averaged rotation rate in the inertial frame is Ω* = Ω + δΩ. In this new differentially rotating environment, the propagation and dissipation of inertial waves are significantly modified, as was also shown in the linear studies of Baruteau & Rieutord (2013) and Guenel et al. (2016b, a), which could be primarily responsible for explaining the discrepancies we find in tidal dissipation rates between linear and non-linear regimes (this will be further discussed in the next section).

The right-hand panel of Fig. 4 shows the tidal dissipation rate as a function of normalized frequency ω/Ω* from non-linear simulations that include all non-linear terms (bottom curves, shifted downwards by a factor 10⁻²) compared with simulations including only the wavelike/wavelike non-linearity (upper curves). Since δΩ evolves with time, particularly in the full non-linear simulations, it means that the simulations also explore a range of tidal frequencies as they progress. We also identify E_dr by the colour of each line, which shows that differential rotation develops in each of these simulations. Especially, the removal of mixed non-linear terms does not prevent the development of zonal flows that are only driven by the wavelike/wavelike non-linearity in that case, as is also shown in Fig. 6. In that figure, we display snapshots in the meridional (r, θ) plane of azimuthally averaged zonal flows for the three frequencies. The energy in the differential rotation is similar or even greater on long time-scales for ω = 1.05 and 1.1, compared with the cases where all non-linearities are included, though the change in the rotation rate δΩ is much less important, as shown in both panels of Fig. 7. Note that non-zero values of δΩ can still occur in simulations with only wavelike/wavelike non-linearities despite the conservation of L_z because the flow is differentially rotating. We have checked that simulations, including only wavelike/wavelike non-linearities, conserve L_z to an accuracy better than 10⁻¹⁰.

Figure 6.

Azimuthally averaged zonal flows 〈u_φ〉_φ for three different frequencies when a steady state is reached in wavelike/wavelike non-linear simulations. The colour maps are invariant with respect to the equator and the rotation axis in each of these simulations.

Figure 7.

Left: energy in the differential rotation E_dr versus time for three initial forcing frequencies and for all (dotted lines) or only wavelike/wavelike (w-w) non-linearities (solid lines). Right: evolution of the mean rotation rate of the fluid normalized by the initial uniform rotation δΩ/Ω versus time for the same non-linear simulations.

In the right-hand panel of Fig. 4, we have also added the normalized constant |$\hat{E}/C_\mathrm{t}^2$| (horizontal dotted lines), which is the energy transfer in an impulsive interaction (equation (112; Ogilvie 2013)

$$\begin{eqnarray} \frac{\hat{E}}{C_\mathrm{t}^2}=\frac{1}{2\pi }\int _{-\infty }^{+\infty }P_\mathrm{t}(\omega)\, \frac{\mathrm{d}\omega }{\omega ^2}=\frac{20}{189}\frac{\alpha ^5}{1-\alpha ^5}\rho R^5\Omega ^2, \end{eqnarray}$$

(28)

with ρ, R, and Ω fixed to unity in our simulations. This quantity gives a measure of a typical dissipation rate over one-time unit (in units of Ω⁻¹) and is close to a (unweighted) mean value of the linear frequency-dependent tidal dissipation across most of spectrum. It is related to the frequency-averaged tidal dissipation of inertial waves, which is commonly applied to study dynamical evolution of exoplanetary systems (e.g. Mathis 2015; Bolmont & Mathis 2016), given by¹¹:

$$\begin{eqnarray} \Lambda =\int _{-\infty }^{+\infty }\mathrm{Im}\left[K^2_2(\omega)\right]\, \frac{\mathrm{d}\omega }{\omega }=\frac{16\pi }{63}\epsilon _\Omega ^2(1+k_2)^2\frac{\alpha ^5}{1-\alpha ^5}, \end{eqnarray}$$

(29)

where |$\epsilon _\Omega ^2=\Omega ^2/\omega _d^2$|⁠. The integrand is the imaginary part of the quadrupolar potential Love number, related to the impulsive energy transfer as |$\Lambda =(12\pi \epsilon ^2_\Omega /(5\epsilon ^2))\, \hat{E}$|⁠.

Final dissipation rates for wavelike/wavelike non-linear simulations shown in the right-hand panel of Fig. 4 (bullets) tend to be closer to the mean tidal dissipation proxy given by equation (28), compared to the highly resonant frequency-dependent linear dissipation (dashed lines), which can be one order of magnitude (or more) lower or higher than the mean value. At the beginning of all non-linear simulations, the energy injected into the differential rotation increases as the dissipation approaches the value predicted by linear theory. Then, the establishment of zonal flows (in particular) causes the dissipation to deviate from its linear prediction, even more or so in the simulations where all non-linearities are included since |$\boldsymbol L$| is not conserved (see the left-hand panel of Fig. 4). As expected, for ω = 1.15 where strong shear layers are not activated and which corresponds to a lower level of dissipation in the linear regime, zonal flows are hardly generated (note the order of magnitude of difference for 〈u_φ〉_φ in Fig. 6, right-hand panel, compared to the other frequencies shown in the left-hand and middle panels). There are no appreciable departures in the tidal dissipation rates from the linear regime in this case. From these three cases, it would appear that the closer we are to a ‘resonant peak’ in which the dissipation is enhanced, the differential rotation excited is stronger and the gap between the linear and wavelike/wavelike non-linear regimes (upper curves in the right-hand panel of Fig. 4) is larger. As shown in Fig. 6, but also in Paper I, Tilgner (2007), Morize et al. (2010), or Subbotin & Shiryaeva (2021), zonal flows produced by non-linear self-interaction of non-axisymmetric inertial waves always take the form of nested and axisymmetric cylinders with different angular velocity profiles depending upon cylindrical radius inside the shell. The radial distribution, magnitude, and direction of the flow at a given location (prograde or retrograde, i.e. red or blue in meridional snapshots of azimuthal velocities, respectively) can be very different depending on the tidal frequency. For instance, high-amplitude prograde flows can develop near the equator or near the rotation axis, as we show in the left-hand and middle panels of Fig. 6, respectively, even though the tidal forcing frequencies are quite similar. While it is difficult to straightforwardly predict the final magnitude and structure of these columnar flows, it is clear that the strongest zonal flows preferentially develop at the points of reflection between the inertial wave beams and the outer boundary or rotation axis. These are locations characterized by large amplitudes in the velocity components or dissipation, as we can see from the meridional cuts of the radial component of the linear velocity in Fig. 5 (top panels) and the dissipation in Fig. 8 based on complementary linear calculations (computed using the spectral code LSB, ``Linear Solver Builder'', Valdettaro et al. 2007). Powerful shear layers emerge from the critical latitude at |$\theta =\arcsin (\omega /(2\Omega))$| tangent to the inner sphere, approximately at latitudes of 32°, 33°, and 34°, respectively, from the left-hand to the right-hand panels of Figs 5 and 8. These are recurring features of forced inertial waves in 3D rotating spherical shells (e.g. Rieutord & Valdettaro 2010, 2018). In meridional cuts displaying the non-linear radial velocity for ω = 1.05 and 1.1 (left-hand and middle bottom panels of Fig. 5), we observe that shear layers are attenuated compared to their linear counterparts (upper panels), presumably due to the development of zonal flows. Though cylindrical zonal flows are the dominant patterns when looking at the φ −average of the azimuthal velocity in non-linear simulations, the imprints of shear layers are still visible in snapshots in the three panels of Fig. 6 when comparing with the same panels showing linear dissipation in Fig. 8. The key role of shear layers and critical latitudes in the establishment of these geostrophic flows have already been pointed out for example in the experimental study of Subbotin & Shiryaeva (2021).

Figure 8.

Viscous dissipation in the meridional plane in linear calculations computed with LSB for three different tidal forcing frequencies ω = 2ω_f. Note the amplitude is arbitrary as this is from a linear calculation.

From now on and throughout the rest of this paper, we analyse simulations with wavelike/wavelike non-linearities only.

3.3 Effects of the zonal flows on tidal dissipation

In this section, we demonstrate that the development of differential rotation can explain the discrepancies between linear and wavelike/wavelike non-linear tidal dissipation rates in many of these simulations, as we observe in the right-hand panels of Figs 3 and 4. At late times in our non-linear simulations (several thousand rotation times), these simulations reach an approximate steady state for K, D_ν, E_dr and δΩ/Ω, in which zonal flows are fully developed and hardly evolve any more. To further analyse these simulations, we extract a φ − and z −average of the azimuthal component of the velocity 〈u_φ〉_{φ, z} when the simulation reaches such a steady state (as done e.g. in Fig. 9, showing this quantity for the three frequencies studied earlier). We then investigate the energy exchanges in new linear simulations of tidally forced inertial waves, where we apply the final zonal flow from non-linear simulations as a ‘background flow’. This allows us to determine whether the most important non-linear effect in our simulations is the generation of differential rotation and the back-reaction of this on the tidally forced waves, as opposed to triadic wave interactions (or parametric instabilities) involving inertial waves, for example. We adopt a similar approach to the linear calculations of Baruteau & Rieutord (2013) with initial cylindrical differential rotation, except that inertial waves are tidally forced in our model and we perform these as an initial value problem using simulations (while they studied free inertial modes as an eigenvalue problem).

Figure 9.

Vertically and azimuthally averaged wavelike azimuthal velocity 〈u_φ〉_{φ, z} (‘output’) at the end of the non-linear simulations (Ωt = 5000) against cylindrical radius for three different forcing frequencies. The fitting curves (‘fit’) are degree 10 polynomials computed using a least-squares method.

In our model, the azimuthal background flow is |$\boldsymbol U=\boldsymbol r\wedge \boldsymbol \Omega =s\Omega _s(s)\, \boldsymbol e_\varphi$|⁠, where the rotation Ω_s(s) is a function of the cylindrical radius¹²s = rsin θ. The rotation can be decomposed into a constant mean rotation Ω plus an s-dependent small departure δΩ_s, such that Ω_s(s) = Ω + δΩ_s(s). Since the numerical simulations are in the corotating frame, the correct initial zonal flow to choose is given by δΩ_s = 〈u_φ〉_{φ, z}/s. Unlike in Baruteau & Rieutord (2013) where a quadratic function of s is used to describe the cylindrical rotation (chosen to maintain matrix sparsity as much as possible), the rotation profile obtained is not a simple function of s in our models (see Fig. 9). We have made the choice to reconstruct the zonal flow 〈u_φ〉_{φ, z} as a degree n polynomial such that |$\langle u_\varphi \rangle _{\varphi ,z}=\sum _{i=0}^{n}a_is^i$| (with n = 10 − 20) and where a_i are fitted constants. In the frame rotating at the rate Ω, the linearized momentum equation for inertial waves excited by the effective tidal forcing |$\boldsymbol f_\mathrm{t}$|⁠, and propagating in our background zonal flow, reads (see also Baruteau & Rieutord 2013):

$$\begin{eqnarray} (\partial _t+\delta \Omega _s\partial _\varphi)\boldsymbol u_\mathrm{w}&+&2\Omega _s\boldsymbol e_z\wedge \boldsymbol u_\mathrm{w}+s\boldsymbol e_\varphi (\boldsymbol u_\mathrm{w}\cdot \boldsymbol \nabla)\delta \Omega _s=-\frac{\boldsymbol \nabla p_\mathrm{w}}{\rho } \\ &+&\nu \Delta \boldsymbol u_\mathrm{w}+\boldsymbol f_\mathrm{t}. \end{eqnarray}$$

(30)

The terms on the left-hand side of equation (30, except for |$\partial _t\boldsymbol u_\mathrm{w}$|⁠) come from the linear development of the inertial terms |$(\boldsymbol U+\boldsymbol u_\mathrm{w})\cdot \boldsymbol \nabla (\boldsymbol U+\boldsymbol u_\mathrm{w})$| in the corotating frame, remembering that |$\boldsymbol U=s(\Omega +\delta \Omega _s)\boldsymbol e_\varphi$| and keeping the first-order perturbation only. Thus, these new terms directly originate from the initial differential rotation imposed here.

We can also derive an energy balance as we did in equation (7) here for linear inertial waves in a cylindrically, differentially rotating background flow (to avoid any confusion, we add the subscript dr):

$$\begin{eqnarray} \partial _tK_\mathrm{dr}= I_\mathrm{w-dr}-D_{\nu , \mathrm{dr}}+P_{\mathrm{t},\mathrm{dr}}. \end{eqnarray}$$

(31)

We have also introduced the energy transfer term (involving Reynolds stresses) between the wavelike and zonal flows:

$$\begin{eqnarray} \begin{aligned} I_\mathrm{w-dr}&=-\langle \rho s u_\varphi (\boldsymbol u_\mathrm{w}\cdot \boldsymbol \nabla)\delta \Omega _s\rangle . \end{aligned} \end{eqnarray}$$

(32)

Since we also have |$I_\mathrm{w-dr}=\langle \rho \delta \Omega (\boldsymbol u_\mathrm{w}\cdot \boldsymbol \nabla)s u_\varphi \rangle$| using the divergence-free constraint and the boundary conditions on the wavelike flow, a negative I_w-dr means energy is transferred from inertial waves to the zonal flow, while a positive I_w-dr means energy is extracted from the zonal flow towards the waves.

In the right-hand panel of Fig. 4, we have added the tidal power (featured by stars) from the new linear simulations solving equation (30) with magic, with the initial zonal shear flow coming from the end of the associated non-linear simulations with the same tidal forcing frequencies ω = 1.05, 1.1, and 1.15 (see Fig. 9). Interestingly, the dissipation rates at the end of the non-linear simulations (when an overall steady state is reached, given by blue circles) match very well with the linear tidal power terms for all three cases. For these simulations, it means that the lower rates of non-linear dissipation compared to linear predictions are fully explained by the setting up of zonal flows inside the shell, and of their effects on the waves. We emphasize that the presence of the transfer term I_w-dr in the energy balance equation (32) means that all of the energy injected by the tidal flows is not entirely dissipated by viscosity, but can also be transferred to (or extracted from) zonal flows.

We have further explored the range of frequencies for which inertial waves can be excited by the tidal forcing, between −1 and 1 in Figs 10 and 11 (with our choice of time unit, we can excite inertial waves between −2 and 2). For a tidal amplitude of C_t = 10⁻², there are no significant departures in non-linear dissipation rates from linear predictions (dashed line), though the simulations having an initial tidal forcing frequency close to a (linear) resonant peak of dissipation have more energy inside the differential rotation and a dissipation rate that differs more from the linear one, as was already observed in Fig. 4 (right-hand panel). Moreover, the tidal power term P_{t, dr} agrees quite nicely in all cases with the final non-linear dissipation. For a tidal amplitude forcing of C_t = 5 × 10⁻², this is no longer systematically true, with simulations having a notable mismatch between D_ν and P_{t, dr} for example, are ω = 0.192, 0.2, and 0.7, or even no steady linear tidal power P_{t, dr}, which still wildly oscillates at late times for ω = 0.8 and 0.9, as we can see in the left-hand panel of Fig. 12. Some of these cases are investigated in further detail in the following.

Figure 10.

Dissipation D_ν as a function of the normalized frequency ω/Ω*. Non-linear dissipation with wavelike/wavelike non-linearities (reddish curves and pale blue bullets indicating the final dissipation) is computed from Ωt = 20 to more than 8000 (when an average steady state is reached) for a tidal amplitude C_t = 10⁻² and initial tidal forcing frequency ω in the range [− 1, 0]. The non-linear dissipation is rescaled by a factor |$C_\mathrm{t}^2$| to compare with the linear predictions for a uniformly rotating body of frequency-dependent (dashed lines). The colourbar indicates the kinetic energy in the differential rotation E_dr. The tidal power rates P_{t, dr} computed from the final steady-state reached in linear simulations with an initial cylindrical zonal flow (corresponding to the one attained in non-linear simulations) are also added as blue stars.

Figure 11.

Same as in Fig. 10 but for initial tidal frequencies spanning the range [0, 1] and for two tidal forcing amplitudes C_t. The simulations for the highest C_t are shifted upwards by a factor 10³ for ease of visualization (upper curves).

Figure 12.

Left: tidal power P_{t, dr} and dissipation D_{ν, dr} versus time Ωt in three linear simulations with a background cylindrical zonal flow. Right: natural logarithm of the kinetic energy K_dr for the same three simulations. The meridional plane snapshots of the kinetic energy in Fig. 16 are at times indicated by the black dots. Black dashed lines at late times are linear fits with slopes ∼1.6 × 10⁻² for ω = 0.8, and ∼8.1 × 10⁻³ for ω = 0.9.

We show in Fig. 13 the tidal power P_{t, dr} and the dissipation D_{ν, dr} for the linear simulations with background zonal flows for which P_{t, dr} differs by more than |$5{{\ \rm per\ cent}}$| from D_{ν, dr}, for the two tidal forcing amplitudes C_t = 10⁻² and 5 × 10⁻² (left-hand and right-hand panels, respectively). In all cases presented in these panels, we observe that the amplitude of the tidal power input is greater than the viscous dissipation, confirming that part of the energy injected by tides is redirected to the zonal flows via the energy transfer term I_w-dr. The difference between P_{t, dr} and D_{ν, dr}, i.e. the value of I_w-dr, is even more important for higher tidal forcing amplitudes (right-hand versus left-hand panel), so for stronger zonal flows.¹³ The maximum differences for steady cases are for ω = 0.824 (more than |$50{{\ \rm per\ cent}}$|⁠), as is shown in the left-hand panel of Fig. 12. In the linear simulations with ω = 0.8 and ω = 0.9, P_{t, dr} > D_{ν, dr} is not satisfied at late times, as soon as the kinetic energy and tidal power start to diverge, implying that energy is being extracted from the zonal flow due to the positive sign of I_{w, dr}.

Figure 13.

Tidal power P_{t, dr} and dissipation D_{ν, dr} versus time Ωt in linear simulations (for tidal forcing frequencies ω of Figs 10 and 11) for which the difference between P_{t, dr} and D_{ν, dr} is larger than |$5{{\ \rm per\ cent}}$| (except for ω = 0.8, 0.824, and 0.9 with C_t = 5 × 10⁻² shown in the left-hand panel of Fig. 12). Left: the tidal forcing amplitude is set to C_t = 10⁻²Right: C_t = 5 × 10⁻².

In fact, these two last cases exhibit an instability; the kinetic energy K_dr grows exponentially soon after the beginning of the simulation, as we can see in the left-hand panel of Fig. 12. We have verified that this behaviour is not due to unphysical numerical effects by varying the spatial and temporal resolution (as well as by using a different time-stepping scheme). We have also run again these linear simulations but with a random noise in the velocity field, but without tidal forcing (namely C_t = 0). We found that such a random noise is not able to destabilize these flows and produce an exponential growth of the kinetic energy. We also point out that no instabilities leading to a diverging kinetic energy are found in the associated non-linear simulations for ω = 0.8 and ω = 0.9. Non-linearities may have a stabilizing effect here by slightly modifying the zonal flows while this is not possible in these linear simulations.

To understand the nature of this instability, we consider first whether our zonal flows could be centrifugally unstable. In the inviscid limit, cylindrical differentially rotating flows depending only on cylindrical radius (albeit which is an imperfect approximation of our flow) can be destabilized by axisymmetric perturbations if the specific angular momentum decreases outwards in magnitude from the rotation axis, or in other words if the Rayleigh discriminant, defined here as

$$\begin{eqnarray} \Phi (s)=\frac{1}{s^3}\frac{\mathrm{d}(s^2\Omega _s)^2}{\mathrm{d}s}, \end{eqnarray}$$

(33)

is negative inside the shell (the condition is necessary and sufficient from the work of Synge 1933; Billant & Gallaire 2005, generalizes this result to local non-axisymmetric perturbations). In Fig. 14, the Rayleigh discriminant Φ is plotted versus s using the zonal flows extracted from the non-linear simulations of Fig. 11 for C_t = 5 × 10⁻². For these columnar flows, Φ is always positive inside the shell, meaning that they are centrifugally stable to any inviscid axisymmetric perturbations.

Figure 14.

Plot of the Rayleigh discriminant Φ versus the cylindrical radius s for the same initial forcing frequencies as in Fig. 11 for a tidal forcing amplitude C_t = 5 × 10⁻². Bullets indicate corotation resonances [equation (35)], while ticks indicate inner critical latitudes and crosses outer critical latitudes, all determined by equation (38).

However, it is unclear whether non-axisymmetric perturbations (like m = 2 inertial waves) are able to destabilize these zonal flows. Considering an inertial mode of complex frequency ω = ω_r + iω_i of real part ω_r and growth (damping) rate ω_i, if this is positive (negative), the temporal dependence of its velocity scales as exp (− iωt), and thus its kinetic energy is proportional to exp (2ω_it). Hence, the ‘linear’ slopes for ω = 0.8 and ω = 0.9 in the plot showing log (K_dr) (right-hand panel of Fig. 12) after the early transient phase, would correspond to twice the growth rate of the triggered unstable mode, here ω_i ≃ 8 × 10⁻³ for ω = 0.8 and ω_i ≃ 4 × 10⁻³ for ω = 0.9. For comparison, the maximum shear rate s∂_sΩ_s for these cases is of the order of 10⁻¹. The knowledge of unstable eigenmodes in the shell with our particular differential rotation profiles could help us to know whether the growth of the kinetic energy in the linear simulations with 0.8 and 0.9 is related to their presence or not. In the analytical study of Dandoy et al. (in preparation), the authors looked at the destabilization of a columnar convective vortex located near the rotation axis that is interacting with (tidal) inertial waves. They found that when vortices are centrifugally unstable under the Rayleigh criterion [equation (33)], tidal inertial waves can trigger the most unstable (non-axisymmetric) mode of the vortex and its destabilization. Since the Rayleigh criterion for stability is satisfied here, such a complex analysis is beyond the scope of this paper and should be postponed to a future study.

Another possible cause of the sharp increase of the kinetic energy observed in some linear simulations with zonal flows is the presence of corotation resonances inside the shell. For the initial m = 2 inertial wave perturbation of frequency ω in the corotating frame (i.e. rotating with Ω), a corotation resonance (or ‘critical layer’) happens when the Doppler-shifted frequency of the wave in the zonal flow frame¹⁴

$$\begin{eqnarray} \sigma =\omega -m\delta \Omega _s, \end{eqnarray}$$

(34)

vanishes, namely at a cylindrical radius s = s_c, when:

$$\begin{eqnarray} \delta \Omega _s(s_\mathrm{c})=\omega /2. \end{eqnarray}$$

(35)

One should note that the definition of the Doppler-shifted frequency equation (34, and thus of the corotation resonance) differs from the one in Baruteau & Rieutord (2013) for example, due to the different sign convention for the wavelike part of the solution [equation (4)] and the different reference frames used (they worked in an inertial frame). In Fig. 15, we show the different zonal flows δΩ_s produced in non-linear simulations for the same tidal forcing frequencies and amplitudes as in Fig 11. We first notice that some flows exhibit one corotation resonance at a critical cylinder s_c satisfying equation (35), usually close to the rotation axis. These corotation resonances are more likely to be present in strong zonal flows, arising when the tidal forcing amplitude is high and the (non-)linear dissipation is strong (e.g. which can be seen by comparing the left-hand and right-hand panels of Figs 11 and 15).

Figure 15.

Zonal flows δΩ_s produced in the non-linear simulations presented in Fig. 11 (left-hand and right-hand panels) in terms of the cylindrical radius s for C_t = 10⁻² (dashed lines) and C_t = 5 × 10⁻² (solid lines). For each zonal flow (of different colour and line-style), the corotation resonance is indicated by a bullet when it exists [i.e. when it satisfies equation (35)].

When present, the corotation resonance seems to have an effect on the zonal flow profile by perturbing it in most cases, producing a small bump around s = s_c. When approaching a corotation resonance, a wave can exchange (deposit or extract) energy (or angular momentum) with the zonal flow, which is in turn modified. It has been analytically and numerically studied in great detail for (inertia-)gravity waves in (rotating) stratified shear flows (e.g. Booker & Bretherton 1967; Grimshaw 1975; Lindzen & Barker 1985; Barker & Ogilvie 2010; Barker 2011; Alvan, Mathis & Decressin 2013; Su, Lecoanet & Lai 2020) and also for inertial waves in a differentially rotating shell or box (Baruteau & Rieutord 2013; Guenel et al. 2016a, b; Astoul et al. 2021). For the simulations for which a corotation resonance exists for both C_t = 10⁻² and C_t = 5 × 10⁻², the location of the resonance s_c appears to move away from the rotation axis for higher tidal forcing amplitudes. As we can see in Baruteau & Rieutord (2013) for cylindrical differential rotation, free inertial waves are often very strongly damped at the corotation resonance. We can thus expect that the successive deposition of energy from inertial waves emitted from the critical latitudes at these corotation points can efficiently contribute to modifying the zonal flow at this location. This process can potentially move the location of the corotation resonance, as has been observed for example in Barker & Ogilvie (2010) and Su et al. (2020) in their local simulations of gravity waves breaking non-linearly.

The role of corotation resonances is clear from Fig. 16, which shows the azimuthally averaged kinetic energy for three linear simulations containing corotation resonances close to the poles at early times. In each of the three panels, inertial waves are emitted from the inner and possibly outer critical latitudes, whose expressions (detailed at the end of the section) are modified by differential rotation. When approaching the corotation resonance at the critical cylinder, shear layers bend as a result of the differential rotation. In the left-hand and right-hand panels (for which the kinetic energy diverge in these linear simulations at later times, see Fig. 12), the kinetic energy is concentrated along the critical cylinder on both sides, while in the middle panel (for which the kinetic energy reaches a steady state at late times) the critical cylinder seems to act like an absorbing barrier for incoming inertial waves which do not cross the corotation resonance.

Figure 16.

Azimuthal average of the wavelike kinetic energy 〈u²〉_φ in the meridional plane in linear simulations with initial zonal flows for the same three tidal frequencies as in Fig. 12 at early times (indicated by black dots on that figure). The critical cylinder corresponding to the corotation resonance is indicated by a vertical dashed line and the critical latitudes θ defined by equation (38) on the inner and outer shells are featured by the inclined black ticks.

Since we observe instabilities inside the simulations (left-hand and right-hand panels of Fig. 16), we first investigate the possibility of an instability similar to the ‘Papaloizou–Pringle instability’ (Papaloizou & Pringle 1984, 1985, 1987), which occurs in Rayleigh-stable (i.e. with Φ > 0, just like here) Keplerian discs, and has been observed e.g. in Barker et al. (2016). For this kind of instability to occur, the inertial wave needs to be sent multiple times into the corotation region to exhibit amplification. Since in our cases, the differential rotation is not strong enough to confine the waves in a portion of the domain¹⁵ (i.e. for turning points to exist in the differential equation governing inertial wave propagation, where it changes character from hyperbolic/wavelike to elliptic/evanescent), it is the spherical boundaries and the rotation axis that could reflect the waves back towards the corotation region. This has the potential to lead to instability if the waves can ‘over-reflect’ and be amplified at corotation.

This phenomenon is also similar to the over-reflection mechanism described in the review of Lindzen (1988) in the case of internal gravity waves in shear flows and leading to shear instability (see also e.g. Jones 1968; Lindzen & Tung 1978; Lindzen & Barker 1985; Baines & Mitsudera 1994; Harnik & Heifetz 2007; Alvan et al. 2013, for Rossby waves). Several conditions on the ‘wave geometry’ are required for such an instability to happen. The bounded domain where the waves propagate must be divided into three regions: two ‘propagative’ regions with wavelike solutions where the Richardson number¹⁶ Ri is greater than 1/4, surrounding a region containing the critical layer where Ri < 1/4 (a necessary condition for instability, Miles & Howard 1964). Under suitable boundary conditions for the waves to be reflected back to the critical cylinder and interfere constructively with waves on the other side of corotation (i.e. with wave quantization on the phase speed or frequency), a shear instability can occur after multiple amplifications of the waves at corotation.

In Astoul et al. (2021), an analogous criterion to the Richardson criterion (also known as the Miles–Howard theorem, and establishing a necessary condition for instability when Ri ≤ 1/4) but for inertial waves has been derived in the context of a local model of cylindrical differential rotation. They demonstrate that over-reflection and over-transmission are also possible for inertial waves in a similar three layer model, which opens the way to possible shear instabilities (though their study did not investigate these instabilities directly). In Appendix  A, we demonstrate that this criterion, based on an inviscid analysis of free inertial waves with cylindrical rotation, may predict the stable regime for all zonal flows investigated in this study. In such case, we expect waves to be strongly damped at corotation, which would rule out shear-type instabilities. However, the derivation of this criterion in Appendix  A and Astoul et al. (2021) does not take into account viscous dissipation and tidal forcing, which could modify the criterion. For gravito-inertial waves, viscosity is able to destabilize stratified shear flows and allow instabilities even if the Richardson number is greater than 1/4 (e.g. Miller & Lindzen 1988). For inertial waves, Baruteau & Rieutord (2013) and Guenel et al. (2016a) numerically observed unstable modes triggered below a certain viscosity threshold for shellular and conical differential rotation. Tidal forcing could also modify this criterion, as highlighted in Appendix  A, but a more detailed analysis is postponed to a future study. In either case (whether it is a shear or ‘Papaloizou–Pringle’-like instability), the mechanism leading to non-axisymmetric instability looks similar: incoming inertial waves launched from the critical latitude(s) (or elsewhere) propagate towards the corotation resonance, partly tunnel through, partly (over-)reflect and are partly damped (Narayan, Goldreich & Goodman 1987; Goldreich & Narayan 1985; Lindzen 1988). The transmitted wave reflects from the rotation axis, and the reflected wave reflects from outer boundary (or core) and back. These waves potentially sustain a positive growth of the kinetic energy where they meet near the corotation resonance (near s = s_c as in Fig. 16).

An important point to stress about the linear simulations involving these kind of instabilities is that the meridional circulation observed in the associated non-linear simulations are not negligible, as assumed in equation (30). In the snapshots showing the φ −average of the non-linear velocity components for ω = 0.8 in Fig. 17, the azimuthal average of the radial velocity 〈u_r〉_φ is not negligible anymore compared to the zonal flow 〈u_φ〉_φ, especially close to the corotation resonance, though it is still weaker in amplitude. This also demonstrates that radial (vertical) angular momentum flux is efficiently deposited at the corotation resonance and contributes to the creation of a strong radial meridional flow there.¹⁷ We also observed strong meridional (and especially radial) flows close to the corotation resonance in the non-linear simulations for ω = 0.7 and ω = 0.9 (not shown in this study). Whether these meridional flows are responsible for the departure of the non-linear dissipation from linear predictions in Fig. 11 is an open question. We emphasize, however, that no strong meridional flows have been observed for the non-linear simulations for ω = 0.192 and ω = 0.2 where a departure in the dissipation rate is also observed. In these simulations, the corotation resonance is extremely close to the rotation axis (see Fig. 15) and does not appear to play an important role in the redistribution and dissipation of inertial wave kinetic energy, as shown in the meridional snapshots and averages of u_r and u_φ displayed in Fig. 18 for ω = 0.2. This is presumably because the inner critical latitude is close to the equatorial plane, from which the waves propagate in almost vertical shear layers and deposit their energy preferentially close to the equator.

Figure 17.

Azimuthal average in one quadrant of the three components of the wavelike velocity |$\boldsymbol u_\mathrm{w}$| at the end of the non-linear simulation with ω = 0.8. As in Fig. 16, the ticks indicate the critical latitudes and the dashed line the corotation resonance.

Figure 18.

Snapshots of one quadrant for the non-linear simulation with ω = 0.2 and C_t = 5 × 10⁻². Critical latitudes and the corotation resonance are indicated as in Fig. 16 (the inner critical is almost in the equatorial plane); crosses are the locations of the points where the Fourier transform of u_r has been performed (in the left-hand panel of Fig. 19). Left: φ-slice of the radial component of the wavelike velocity u_r at early time. Middle: φ-average of the azimuthal component of the wavelike velocity 〈u_φ〉_φ at the same time. Right: φ-average of the azimuthal component of the wavelike velocity 〈u_φ〉_φ at a later time.

We finally compute the location of the critical latitudes that differ from the uniform rotation case when the zonal flows are strong. The critical latitudes θ for which rays (inertial wave shear layers) are tangent to the shell’s inner core and outer surface satisfy the equation in the first quadrant in the meridional plane

$$\begin{eqnarray} \frac{\mathrm{d}z}{\mathrm{d}s}=-\frac{s}{z}, \end{eqnarray}$$

(36)

at |$r=\sqrt{s^2+z^2}=\alpha R$| or R, respectively. Moreover, following Baruteau & Rieutord [2013, equation (3.8)], the paths of characteristics for inviscid and free inertial waves are defined here by:

$$\begin{eqnarray} \frac{\mathrm{d}z}{\mathrm{d}s}=\pm \sqrt{\frac{\kappa ^2(s)}{\sigma ^2(s)}-1}, \end{eqnarray}$$

(37)

with |$\kappa =\sqrt{\Phi }$|⁠, the epicyclic frequency (commonly used in astrophysical discs). Combining equations (36) and (37), we obtain the expression for the critical latitudes in the first quadrant:

$$\begin{eqnarray} \sin \theta =\frac{\sigma (s)}{\kappa (s)}, \end{eqnarray}$$

(38)

at r = αR for the inner shell and r = R for the outer shell, with θ measured from the equator, where z = αRsin θ, s = αRcos θ describes the inner critical latitude (and z = Rsin θ, s = Rcos θ the outer critical latitude). One should note that for weak amplitudes of the zonal flow and shear, σ ≃ ω, κ ≃ 2Ω, i.e. Φ ≃ 4 in our time unit, as we can observe in Figs 14 and 15 (it concerns notably the outer critical latitudes), and equation (38) tends towards the classical definition of the critical latitude with uniform rotation sin θ = ω/2Ω. For strong enough zonal flows, the difference between equation (38) and the latter definition is quite important (note also that inner and outer critical latitudes are the same with uniform rotation). For example, for a tidal forcing frequency ω = 0.8, θ ≃ 23.6° with uniform rotation, while the inner critical latitude is θ ≃ 39.6° and the outer critical latitude is θ ≃ 20.4° using equation (38) derived with the zonal flows at this frequency (note that the outer critical latitude is closer to the critical latitude with uniform rotation for the reason explained above). These inner and outer critical latitudes nicely match the locations of the high amplitudes in kinetic energy in Fig. 16 (black ticks across inner and outer surfaces show the predicted θ’s).

Figure 19.

Left: Fourier transform of the radial component of the non-linear wavelike velocity u_r for two points in the shell indicated by crosses in Fig. 19 and for the time range t ≤ 6000. Right: Fourier transform of the poloidal velocity potential W for the whole time range t ≤ 15 000, with azimuthal order m and degree l indicated in colour. Dominant frequencies excited are indicated by dashed grey lines and labelled.

3.4 Parametric instability of inertial waves

In this section, we focus on the non-linear simulation with ω = 0.2 and C_t = 5 × 10⁻². After an early transient phase, the kinetic energy of this simulation reaches a steady state, and then an instability sets in producing oscillations with several different periods for t ≳ 10 000, as we can see on the left-hand panel of Fig. 20. To understand what is happening here, we analyse the azimuthal kinetic energy spectrum in the right-hand panel of Fig. 20 as a function of azimuthal order m at different times (before and after the oscillations start), which we also compare with the C_t = 10⁻² simulation. When C_t = 10⁻², m = 2 (of the initial tidally forced waves) is the dominant component in the energy spectrum, but when C_t = 5 × 10⁻², the m = 0 component is of the same order of magnitude (both at early and late times), implying a strong zonal flow, and m = 1 is also strong at late times.

Figure 20.

Left: kinetic energy versus time Ωt in two non-linear simulations with different forcing amplitudes C_t with ω = 0.2. Right: normalized (to the maximum value for each simulation) kinetic energy spectrum in terms of the azimuthal order m + 1 for ω = 0.2 for two tidal amplitudes C_t, at different times. The m⁻³ (red) and m⁻² (black) scaling laws, found for example in numerical simulations of Barker & Lithwick (2013, for elliptical instability with strong geostrophic flows) and Le Reun, Favier & Le Bars (2020, for ‘inertial wave turbulence’), have been plotted for reference, although neither is a particularly good fit to our data. Note that odd m values of the kinetic energy for the orange and blue curves are of the order of numerical errors or lower (less than 10⁻¹²), so their contributions are negligible and mostly not shown.

We examine in further detail the Fourier transform of the non-linear radial velocity u_r (after applying a Hamming window function) at two points inside the shell (indicated by crosses in Fig. 18): one in the equatorial plane where the amplitude of the zonal flow is strong, the other in a shear layer directly emerging from the inner critical latitude. This is done in the left-hand panel of Fig. 19 at the early-time indicated. As expected, the contribution of the tidal forcing at ω = 0.2 is dominant for the point inside the shear layer. Interestingly, superharmonics of the frequency ω are also excited before the instability is triggered, with particularly high power in the region where the zonal flow is strong (point 1) compared to inside the shear layer (point 2). These oscillations have a frequency nω and azimuthal order mn = 2n with n ≤ 10, a positive integer which must satisfy nω ≤ 2 in our time units, corresponding to the upper limit of propagation of inertial waves. We also observe superharmonics in the simulation for C_t = 10⁻². Superharmonics have been clearly reported for gravity waves (e.g. Baker & Sutherland 2020; Boury, Peacock & Odier 2021; Ivanov, Chernov & Barker 2022), but have not been observed in many prior studies of inertial waves, though they have been predicted analytically, and observed, for instance in Barik et al. (2018). One reason could be the limited frequency range allowed for inertial waves, imposing |ω| ≪ 2Ω to observe multiple higher harmonics of the main frequency ω.

In the right-hand panel of Fig. 19, we perform a Fourier transform (applying again a Hamming window) of the poloidal velocity potential W(l, m) at the mid-shell for the whole time range to explore the properties of the waves (frequency, m and degree l) excited after the instability becomes important for t ≥ 10 000. In addition to the primary wave with a frequency ω and its (self-excited) superharmonics, we also observe a strong peak at the frequency ω₁ ≃ 0.13 with wavenumbers (l, m) = (1, 1), along with a number of even l + m modes (due to symmetry; with odd l + m being excited for the toroidal velocity potential). This provides clear evidence for triadic resonances (theorized and observed for inertial waves in e.g. Kerswell 1999; Barik et al. 2018) since all of the excited modes can be recovered using only linear combinations of ω and ω₁, and their corresponding azimuthal orders and degrees. In particular, we have

$$\begin{eqnarray} \begin{aligned} &\omega _\mathrm{a}\pm \omega _\mathrm{b}=\omega _\mathrm{c}, &m_\mathrm{a}\pm m_\mathrm{b}=m_\mathrm{c},\ \text{ and }\ l_\mathrm{a}\pm l_\mathrm{b}=l_\mathrm{c}, \end{aligned} \end{eqnarray}$$

(39)

with for example a parent primary wave ω_a = ω with (l_a, m_a) = (2, 2), interacting with two secondary/daughter waves, one with ω_b = ω₁ with (l_b, m_b) = (1, 1), and another with ω_c = ω − ω₁ = ω₂ with (l_c, m_c) = (l_a − l_b, m_a − m_b) = (1, 1). In this case, we have a parametric instability of the primary wave involving these three components, as is shown in the right-hand panel of Fig. 19. It should be noted that there are also superharmonics of these daughter waves (for example with frequencies 2ω₂ and 2ω₁).

The exponential growth of the parametric instability involving these m = 1 modes is clear from Fig. 21, which shows the time-series of the logarithm of the poloidal velocity potential for (l, m) = (1, 1), (2,2), and (3,3). Indeed, the dashed green curved with linear slope ∼0.033 indicates that the secondary waves with wavenumbers (1,1) and (3,3) are exponentially growing at the expense of the primary wave by draining its energy, until a saturation level is reached after Ωt ≳ 10 500. Thanks to viscous dissipation, the simulation reaches a quasi-steady state, similarly to Cui & Latter (2022), where an inertial mode also triggers wave–wave interactions and parametric instabilities in their local simulations of protoplanetary discs.

Figure 21.

Time series of the poloidal velocity potential log (W) (base e) for different modes (note that the same (l, m) component can relate to different frequencies as in Fig. 19). The dashed green line illustrates the exponential growth of the modes of components (l, m) = (1, 1) and (l, m) = (3, 3). This indicates that the primary tidal wave is unstable to a slowly growing parametric instability.

Parametric instabilities typically transfer energy to waves excited close to reflection points, with lower frequencies and shorter wavelengths, as observed in non-linear simulations of inertial wave attractors in 2D Cartesian geometry (Jouve & Ogilvie 2014). It is possible that the mode of frequency ω₁ is related to this phenomenon. As the Ekman number is lower (for which non-linearities are expected to be enhanced), the frequency of the secondary could become closer to exact subharmonics with ω/2 for the fastest growing modes (though this has not always been found in related problems in spherical geometry, e.g. Lin, Marti & Noir 2015, where ω/3 and 2ω/3 have been obtained instead for moderate Ekman numbers, similar to our case), as obtained in the inviscid theoretical parametric subharmonic instability study in the aforementioned paper, and as also expected for the related elliptical instability (e.g. Kerswell 2002; Barker 2016, which is usually analysed for simple elliptical flows). This kind of instability was not observed in Paper I, perhaps because of the strong zonal flows that dominated their simulations.

It is clear from the energy spectrum in the right-hand panel of Fig. 20, that kinetic energy is cascaded towards high values of azimuthal wavenumbers m (when comparing blue and orange curves at different tidal forcing amplitudes). This is initially due to superharmonics of the main frequency ω, that correspond with even m (orange), and then to odd values of m (green curve) presumably due to the triadic resonances with secondary waves with m = 1, which subsequently initiate an inertial wave turbulent cascade (e.g. Le Reun et al. 2017).

Finally, we emphasize that undertaking the eigenvalue problem for this differentially rotating shell¹⁸ could help us to understand which modes are likely to be excited and cause the particular instability we observe.

3.5 Scaling laws for tidally generated differential rotation: effects of the viscosity

In this section, we analyse how the energy inside the tidally generated differential rotation, E_dr, varies as we vary the viscosity/Ekman number E. In Fig. 22, we display E_dr in terms of E when an overall steady state is reached in our non-linear simulations for the three frequencies studied in Section 3.2. For the simulations with E ≤ 5 × 10⁻⁶, the spatial resolution is increased to l_max = 170 (so n_φ = 512) and n_r = 161.

Figure 22.

Energy inside the differential rotation E_dr as a function of the Ekman number E for various non-linear simulations with three tidal forcing frequencies ω (with C_t ≈ 9 × 10⁻³). Simulations are run until an overall steady state is reached and the values taken are the mean in the last 3000 rotation times (and errorbars correspond to minimum and maximum values in that range, being non-negligible only for one point). Some scaling laws have been added as dashed lines for comparison.

Looking at E_dr for ω = 1.05 and 1.1 (for which the linear dissipation is close to a resonant peak), we observe that E_dr ∝ E⁻², while for ω = 1.15 (for which the linear dissipation is at a trough), the scaling is less steep, closer to E_dr ∝ E⁻¹. At low Ekman numbers, such as E ≤ 5 × 10⁻⁶, there is a break in the E⁻² law for ω = 1.05. At such low viscosities, the amplitude of prograde zonal flows starts to become large (for E = 5 × 10⁻⁶) or even dominant (for E = 2 × 10⁻⁶) close to the rotation axis for the ω = 1.05 and 1.15 cases, contrary to the ‘equatorially’ dominant flow we found for the same frequencies in Fig. 6 for E = 10⁻⁵. This is presumably explained by the development of a corotation resonance close to the pole, since the amplitude of the zonal flow is anticorrelated with the viscosity. For the ω = 1.1 case, since the zonal flow is already strong near the pole, the presence of the corotation resonance may have amplified it, which is what we suspect for E = 5 × 10⁻⁶.

It is difficult to determine theoretically a universal scaling law, given the strongly frequency-dependent predictions from linear theory. However, Tilgner (2007) argued for an upper threshold E_dr < E⁻²ℓ⁻², which predicts E_dr ∝ E⁻² as we have observed if ℓ = O(1), but is steeper than E⁻² if we naively use ℓ ∝ E^1/3 here (or indeed any positive power of E). Following similar lines to Tilgner’s analysis, who distinguished a momentum equation for the m = 2 tidal inertial waves |$\boldsymbol u_1$| (his equation 3) and one for the axisymmetric flow |$\bar{\boldsymbol u}_2$| generated by their non-linear self-interaction (equation 5), we can derive |$\mathrm{E}\hat{\boldsymbol \varphi }\cdot \boldsymbol \nabla ^2\bar{\boldsymbol u}_2\sim \hat{\boldsymbol \varphi }\cdot (\boldsymbol u_1 \cdot \boldsymbol \nabla) \boldsymbol {u}_1$|⁠, with |$\hat{\boldsymbol \varphi }$| the unit vector in the azimuthal direction. If we assume the tidal waves to have a typical (transverse) length-scale ℓ and the zonal flow to have a radial scale L, this would suggest the zonal flow magnitude to scale as |$\bar{u}_2\sim u_1^2 L^2/(\mathrm{E}\ell)$|⁠, in terms of the typical velocity magnitude of the waves. Hence, |$E_\mathrm{dr}\sim \langle \bar{u}_2^2\rangle \sim \langle u_1^4L^4/(\mathrm{E}^2 \ell ^2)\rangle$|⁠, which predicts E_dr ∝ E⁻² if u₁, L, ℓ, and the volume-filling fraction (VFF) of the waves combine to produce an O(1) number in the above scaling. On the other hand, if we assume u₁ ∼ E^−1/6, L ∼ ℓ ∼ E^1/3 and VFF also scales as E^1/3, we find the less steep scaling E_dr ∝ E^−5/3 (reducing to E_dr ∝ E⁻² under the same assumptions if VFF is O(1)). If we instead focus on flows generated near the critical latitude, and consider u₁ ∼ E^−1/5, L ∼ ℓ ∼ E^2/5 and VFF scaling like E^2/5, we would obtain E_dr ∝ E^−8/5. Different assumptions would lead to alternative scalings, but these are hard to justify a priori.

The Ekman number scalings we have observed are very different from the zonal flows generated by tidal forcing in a full sphere in the experiments (with a deformable no-slip boundary) of, for example, Morize et al. (2010), where the zonal velocity scales as ϵ²E^−3/10, or for those produced by libration-driven inertial waves, as studied in Cébron et al. (2021) and Lin & Noir (2021), where it instead scales as ϵ²E⁰ or ϵ²E^−1/10, respectively (though these also result from self-interaction, so the zonal velocity also scales as ϵ²). We underline that the best-fitting laws here seem quite different from the ones emerging in Paper I in E^−3/2 and E^−1/2, presumably because of the non-wavelike non-linearities involved in their work. It is difficult to predict whether the scaling laws we have observed would still hold for solar-like and Jupiter-like values of the Ekman number, particularly because the corotation resonance at δΩ_s = ω/2 near the rotation axis starts to trigger instabilities for E = 10⁻⁶ for the three frequencies studied in this section. These instabilities appear to flatten the dependence, and could even cause E_dr to become independent of E below a critical value. Further work is required to understand these results theoretically and to explore a wider range of parameters.

4 IMPORTANCE OF NON-LINEAR EFFECTS FOR EXOPLANETARY SYSTEMS

In parallel with this theoretical modelling, it is important to consider which astrophysical systems are likely to be affected by these non-linear effects. Using the time and length-scales defined in Section 2.2, the importance of wavelike non-linearities can be quantified by:

$$\begin{eqnarray} \frac{(\boldsymbol u_\mathrm{w}\cdot \boldsymbol \nabla)\boldsymbol u_\mathrm{w}}{\partial _t\boldsymbol u_\mathrm{w}}\sim \frac{u_\mathrm{w}}{\omega \ell }\sim C_\mathrm{t}E^{-\alpha -\beta }, \end{eqnarray}$$

(40)

where −α − β = −1/2 or −5/12 for shear layers, and −α − β = −3/5 for the critical latitudes. Note that as hypothesized in Paper I, wave breaking is expected when the right hand side is O(1), which is more likely to happen at the critical latitudes for small E, perhaps preventing the propagation of shear layers from these locations (see also Goodman & Lackner 2009). We confirm that non-linear effects are important for both shear layers and critical latitudes according to equation (40) in most of our simulations (i.e. for E ≲ 10⁻⁵ and C_t ≳ 10⁻²), in which the RHS is ≳ O(1).

An estimation of equation (40) in compact exoplanetary systems can be made using typical values for the stellar or planetary tidal amplitude ϵ (depending on whether the primary body in which tidal flows are investigated is the star or the planet), which we display in Fig. 23. We recall that C_t and ϵ differ only by k₂, which is around 1/2 for giant gaseous planets, like Jupiter and Saturn, and also possibly for Hot Jupiters (e.g. Guillot et al. 2022). For a n = 1 non-rotating polytrope, a similar value is found, but can vary with rotation to between 0.3 and 1.2 (Dewberry & Lai 2022). The planetary tidal amplitude parameter varies in Fig. 23 between 10⁻⁶ and 10⁻¹, with maximum values reached for the Hot Jupiters WASP-12 b and WASP-19 b at ϵ ≈ 5 × 10⁻². Using a typical value for the kinematic Ekman number of about E = 10⁻¹⁸, and C_t ≈ 3ϵ/2, non-linear effects for wavelike tidal flows are predicted to matter for every gaseous planet according to equation (40) for any of our choices of exponent −α − β.

Figure 23.

Stellar and planetary tidal amplitudes ϵ in terms of the orbital period P of the planet in days, and the mean density |$\bar{\rho }=M_1/(4/3\pi R^3)$| of the star or planet in colour. Only compact systems satisfying P < 10 d, and M_p/M_⋆ > 10⁻⁴ (bodies with |$\rho \gt 100\, \mathrm{g\,cm^-3}$| have also been removed) have been selected (a significant fraction of these are Hot Jupiter systems), and radii, masses, periods, and semimajor axes have been taken from the online data base http://exoplanet.eu/.

The stellar tidal amplitude parameter varies between 10⁻⁸ and 10⁻⁴, with a maximum value for the host star WASP-19, in which ϵ ≈ 2 × 10⁻⁴. Assuming a Sun-like kinematic Ekman number of about E = 10⁻¹² (probably even lower at the base of the convective envelope), non-linear effects are significant for only a fraction of low-mass stars (preferentially those with ϵ ≳ 10⁻⁵), depending on the value of −α − β. We point out that the results in our paper may be more relevant to fast-rotating bodies like Hot Jupiters, which possibly exhibit cylindrical-like differential rotation in their upper atmospheres (e.g. Gastine, Wicht & Aurnou 2013), than to host stars for which we generally expect a more conical-like differential rotation in the convective envelope (e.g. Benomar et al. 2018), like in the Sun. Although extrapolating E_dr for much lower Ekman numbers in Fig. 22 suggests strong shears and thus a much higher impact on the tidal dissipation (than observed e.g. in Fig. 11), these statements must be tempered by the unknown action of turbulent convective motions on tidal inertial waves, which potentially drive even stronger zonal flows, and the presence of magnetism and fluid instabilities that could mitigate them.

If we assume that the action of convection on tidal flows can be modelled by an eddy viscosity that behaves just like a microscopic viscosity, as done for instance to study Rossby and inertial waves in the Sun in Gizon, Fournier & Albekioni (2020) and Fournier et al. (2022), and which is in the order of 10⁻⁴ − 10⁻⁵ from mixing length theory, major differences in outcome are to be expected for the scaling law equation (40), and tidal dissipation estimates from our simulations may be more directly applicable. With a turbulent eddy viscosity in the order of 10⁻⁵, like in most of our simulations, our numerical results first suggest that for all low-mass host stars and Hot Jupiters with ϵ ≲ 10⁻² (see Figs 4, 10, and 11 for C_t = 10⁻²), non-linear effects have a limited effect on the tidal dissipation, modifying the frequency-dependent (and presumably frequency-averaged) tidal dissipation rates by a factor of unity only. However, for Hot Jupiters having ϵ ≳ 10⁻² (depending on their Love number), and especially for ultra Hot Jupiters like WASP-12 b, WASP-19 b, WASP-121 b, or HIP 65A b having C_t ≳ 5 × 10⁻², non-linearities could have a much greater impact on tidal dissipation rates and the generation of zonal flows, depending on the forcing frequencies in these systems, the exact value of the Ekman number inside the convective envelope, and also the aspect ratio of the shell. Further exploration, and scanning of the frequency range allowed for inertial waves, could allow us to determine if the frequency-averaged non-linear dissipation is modified in comparison to its linear analogue (Ogilvie 2013).

5 CONCLUSIONS AND DISCUSSION

In this paper, we have investigated the impact of non-linearities on tidal flows in an incompressible and adiabatic spherical shell, which models the convective envelope of a low-mass star or a giant gaseous planet in a compact system. In our study, we have performed a suite of 3D hydrodynamical direct numerical simulations in spherical shell geometry for a range of tidal frequencies and amplitudes, and fluid viscosities (Ekman numbers), building upon the prior study of (Paper I, Favier et al. 2014). We show in particular that non-linear tidal effects are likely to be important in the convective envelopes of planets and stars in Hot Jupiter systems.

Unlike in the non-linear simulations of Paper I where the flow was forced from the outer surface through an imposed radial velocity, here we use a tidal body forcing to excite inertial waves, resulting from the residual action of the Coriolis acceleration on the equilibrium (non-wavelike) tide. We consider this to be a more realistic way to tidally excite inertial waves, and it has the advantage that it allows us to clearly identify and assess the contribution of wavelike (⁠|$\boldsymbol u_\mathrm{w}\cdot \boldsymbol \nabla \boldsymbol u_\mathrm{w}$|⁠), mixed (⁠|$\boldsymbol u_\mathrm{w}\cdot \boldsymbol \nabla \boldsymbol u_\mathrm{nw}+\boldsymbol u_\mathrm{nw}\cdot \boldsymbol \nabla \boldsymbol u_\mathrm{w}$|⁠), and non-wavelike (⁠|$\boldsymbol u_\mathrm{nw}\cdot \boldsymbol \nabla \boldsymbol u_\mathrm{nw}$|⁠) non-linear terms in the energy and momentum balances. We demonstrate that the mixed non-linearities generate a non-physical radial flux through the boundaries due to our spherical (and not tidally elliptically deformed) boundaries, which is responsible for the unrealistic angular momentum evolution leading to the de-synchronization of the body observed for some frequencies in Paper I. By removing these non-linear terms, which is justified by physical scaling arguments for astrophysical parameter regimes, the angular momentum is instead conserved in this model. Our model thus probes the instantaneous energy transfers between tidal inertial waves and zonal flows on times much shorter than long tidal evolutionary time-scales.

Like in Paper I and Tilgner (2007), we report the development of strong axisymmetric zonal flows describing cylindrical differential rotation in the shell. These flows are generated by the non-linear self-interaction of tidal inertial waves in various locations in the shell: inside shear layers, and developing around critical latitudes, at the points of reflection between the wave beams and the boundaries or the rotation axis, and finally near a corotation resonance. We observe the formation of corotation resonances (also called critical layers) in our simulations, namely critical cylinders where the angular velocity matches the angular pattern speed of a wave (as in the 2D linear calculations of Baruteau & Rieutord 2013), which form preferentially close to the rotation axis if the tidal forcing is strong or the viscosity is small. This is explained by the fact that the zonal flows in these cases are stronger, since the energy inside the differential rotation varies as E_dr ∝ ϵ⁴E^−γ, where ϵ is the tidal amplitude and γ is a positive number (typically γ ∈ [1, 2], though we find indications that γ → 0 as E → 0) depending on the tidal frequency. When present, these corotation resonances strongly modify the zonal flow profile and the tidal angular momentum exchanges in the system.

We also find that the non-linear tidal dissipation rates depart from linear predictions, with a large discrepancy when the tidal frequency is close to a resonant peak of enhanced dissipation according to linear theory. These cases are also correlated with stronger energy inside the differential rotation. However, unlike Paper I, the discrepancy between linear and non-linear tidal dissipation rates is in general somewhat less important in our simulations (we typically find less than one order of magnitude differences between the two). We demonstrate, by injecting the zonal flows resulting from non-linear simulations into new tidally forced linear simulations as a ‘background flow’, that the development of cylindrical differential rotation is responsible for this discrepancy between non-linear and linear dissipation rates in many of our simulations. This is especially the case for moderate values of the tidal forcing C_t = 10⁻² and Ekman number E = 10⁻⁵.

For higher values of the tidal forcing amplitude, or lower values of the Ekman number, we observe the emergence of complex wave/wave and wave/zonal flow interactions that we try to characterize. In particular, we identify parametric instabilities involving triadic resonances between inertial waves, which are reminiscent of the elliptical instability (e.g. Kerswell 2002; Barker 2016; Le Reun et al. 2017, albeit for a primary tidal wave with a more complicated flow than one with simple elliptical streamlines). We also identify the appearance of corotation resonances, which may trigger (secondary) shear instabilities in some cases, and which lead to strong absorption of waves in other cases. These interactions contribute significantly to change tidal dissipation rates, either by absorbing or amplifying waves at corotation resonances (see also Astoul et al. 2021), or through generating daughter waves that dissipate on smaller length-scales (e.g. Jouve & Ogilvie 2014), and thus redistribute energy inside the shell. Analytical investigation of non-linear three-wave interactions in Kerswell (1999) and Barik et al. (2018) notably emphasize three different outcomes of non-linear self-interaction of inertial waves that have all been observed in our model: a geostrophic m = 0 mode, which is a quasi-axisymmetric zonal flow, superharmonics, as we observe for low frequencies compared to 2Ω, and triadic resonances between two daughter waves and a primary tidal wave.

We have demonstrated that non-linear effects are likely to play an important role on tidal flows, and in modifying tidal dissipation rates, particularly when considering realistic values of atomic viscosities. If turbulent convective motions can be modelled as an effective viscosity on tidal inertial waves, and are thus the dominant contribution for the viscous dissipation term, the latter statement may only be true for the envelopes of ultra short-period Hot Jupiters. However, these statements have to be qualified since various different key processes have not been investigated in this study, including the influence of magnetism and proper turbulent convective motions, that could shape different amplitudes and profiles of the differentially rotating background flow, the effect of compressibility, and also the size of the convective shell. The latter is likely to be important because the linear tidal dissipation prediction with uniform rotation depends strongly on this parameter, with a frequency-averaged dissipation depending on α⁵/(1 − α⁵). Although this implies that tidal dissipation is higher for a thin convective shell, Baruteau & Rieutord (2013) have shown that even for a tiny core, differential rotation is able to retain shear layers, which potentially make a huge difference in tidal dissipation rates. There remains much to be explored in this problem.

ACKNOWLEDGEMENTS

This research has been supported by STFC (Science and Technology Facilities Council) grants ST/R00059X/1, ST/S000275/1 and ST/W000873/1. Simulations were undertaken using the magic software on ARC4, part of the HPC (High Performance Computing) facilities at the University of Leeds, and using the DiRAC (Distributed Research utilising Advanced Computing) Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS (Business, Energy & Industrial Strategy) capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. We are grateful to T. Gastine for his help and valuable advice on the use of magic. We would like to thank M. Rieutord, R. V. Valdetarro, and C. Baruteau for allowing us to use the linear spectral code LSB. We also thank C. Baruteau, J. Park and S. Mathis for fruitful exchanges about the paper, and the referee for a prompt, helpful, and constructive report. This research has made use of data obtained from or tools provided by the portal exoplanet.eu of The Extrasolar Planets Encyclopaedia, and of NASA’s Astrophysics Data System Bibliographic Services.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

1

e.g. from the The Extrasolar Planets Encyclopaedia data base http://exoplanet.eu/readme/guideline.

2

The characteristic frequency associated with buoyancy forces.

3

Due to the incompressibility, the definitions of the equilibrium tide made by Zahn (1966a) and Goldreich & Nicholson (1989) are equivalent, since the non-wavelike flow as defined by equation (4) hereafter is both divergence-free and curl-free.

4

For example E = 10⁻¹¹ − 10⁻¹² for the Sun (e.g. Charbonneau 2013), while E = 10⁻¹⁸ − 10⁻¹⁹ for Jupiter (e.g. Gastine et al. 2014).

5

It has been introduced in the analytical studies of Zahn (1989 ,1966a, b) and Goldreich & Nicholson (1977) using mixing length theory, and supported by the numerical analysis of Duguid, Barker & Jones (2020a, b, in local box models) and Vidal & Barker (2020a, b, in global spherical models).

6

It is the latitude where the shear layer emerges from the surface of the inner boundary.

7

Available at https://magic-sph.github.io/.

8

Typically for simulations with E ≳ 10⁻⁵ and A ≲ 5 × 10⁻².

9

On the order of a few hours with our choice of time-step with 96 cores to reach several thousand rotation times.

10

It is, however, not clear to what extent this relation holds in our non-linear simulation for the lowest frequency, where the zonal flow leads to both a detuning effect by shifting the frequency of the wave away from the resonant linear mode, but also presumably changing the location of the resonant peak itself (e.g. Guenel et al. 2016a).

11

For a fully homogeneous body, k₂ = 3/2, and one recovers equation (113) of Ogilvie (2013).

12

Any meridional circulation is neglected since 〈u_r〉_φ and 〈u_θ〉_φ have been found to be negligible in most of the non-linear simulations, or at least they are weaker (with such cases to be emphasized in the following), compared to 〈u_φ〉_φ.

13

This might be partly because inertial waves are primarily damped near a corotation resonance, as we investigate later in this section.

14

A wave with a frequency |$\tilde{\omega }$| in the inertial frame has a frequency in the Ω-rotating frame of |$\omega =\tilde{\omega }-m\Omega$|⁠, and in the zonal flow frame, it is |$\sigma =\tilde{\omega }-m\Omega _s$| which gives equation (34).

15

The condition for turning surfaces to exist is not met for the zonal flows presented here i.e. where σ² = Φ [see equation (37)], and hence we only observe what Baruteau & Rieutord (2013) refer to as ‘D modes’.

16

It is the ratio of the squared Brunt–Vaäsälä frequency with the squared shear rate.

17

This also has been derived and observed in the pre-liminary study of V. Skoutnev et al., as part of the Kavli Summer Program in Astrophysics 2021 (https://kspa.soe.ucsc.edu/sites/default/files/KSPA_VSkoutnev.pdf).

18

Note that the eigenvalues are known to drift from the uniformly rotating case (e.g. Barik et al. 2018) depending on the strength of the shear and also with the Ekman number (e.g. Bekki, Cameron & Gizon 2022).

19

Equation (A1) can also be derived from Baruteau & Rieutord [2013, equation (3.5)] by taking the leading-order term in s − s_c at the corotation resonance at the critical cylinder s_c.

APPENDIX A: INERTIAL WAVES NEAR COROTATION

We recall that the linear momentum equation (30) for inertial waves with a background zonal flow is analogous to equation (2.1) of Baruteau & Rieutord (2013) with cylindrical differential rotation. The differences are that they omit tidal forcing and work in an inertial frame. By assuming a plane wave in the vertical and azimuthal directions, we can derive a similar second-order ordinary differential equation (ODE) for the free (unforced) pressure perturbation p_w in the inviscid limit (as in their Appendix  A). Under the short-wavelength approximation, and in the vicinity of the corotation resonance at the critical cylinder s_c, this ODE reads¹⁹

$$\begin{eqnarray} \frac{\partial ^2p_\mathrm{w}}{\partial s^2}+\frac{\alpha _c^2(1+\mathrm{Ro}_\mathrm{c})}{\mathrm{Ro}_\mathrm{c}^2(s-s_\mathrm{c})^2}p_\mathrm{w}=0, \end{eqnarray}$$

(A1)

with |$\mathrm{Ro}_\mathrm{c}=\cfrac{s_\mathrm{c}\partial _s\Omega _s(s_\mathrm{c})}{2\Omega _s(s_\mathrm{c})}$|⁠, which is a kind of Rossby number, and the ratio of the vertical to azimuthal wavenumbers α_c = s_ck_z/m. We should point out that in our simulations the singularity at the critical cylinder in equation (A1) is regularized by viscosity, which is not taken into account here. Equation (A1) is analogous to the second-order equation derived in Astoul et al. [2021, equation (68)] for free inertial waves near a corotation resonance with cylindrical differential rotation in a local model (the different sign in the numerator results from the opposite local orientation of the box). Two regimes emerge depending on the value of

$$\begin{eqnarray} \mathcal {R}=\frac{\alpha _\mathrm{c}^2(1+\mathrm{Ro}_\mathrm{c})}{\mathrm{Ro}_\mathrm{c}^2}, \end{eqnarray}$$

(A2)

in equation (A1), which plays an equivalent role to the Richardson number for internal gravity waves. Similarly to the Miles–Howard theorem, if |$\mathcal {R}\gt 1/4$|⁠, solutions of equation (A1) take the form of wavelike solutions, which are preferentially strongly damped when reaching the critical cylinder s_c, while if |$\mathcal {R}\lt 1/4$|⁠, equation (A1) admits hyperbolic solutions, meaning that waves can be amplified. The inequality |$\mathcal {R}\lt 1/4$| is also a necessary but not sufficient condition to trigger shear instabilities, that may depend on the shear flow profile, the boundary conditions, and the properties of the incident waves.

The main difficulty to evaluate |$\mathcal {R}$| in our simulations is to estimate the vertical wavenumber k_z, which should be taken close enough to s_c but in a region where equation (A1) admits wavelike solutions (i.e. when |$\mathcal {R}\gt 1/4$|⁠). The wavenumber transverse to a shear layer scales as E^−1/3 (e.g. Rieutord & Valdettaro 2018), thus the vertical wavenumber k_z close to the inner critical latitude θ scales as E^−1/3sin θ. Of course, this is only valid close to the inner critical latitude and does not account for the bending of the shear layer as it approaches the corotation resonance, where a WKBJ analysis (like in Baruteau & Rieutord 2013; Guenel et al. 2016a; Astoul et al. 2021) predicts the vertical wavenumber to go to zero, or the transverse wavenumber k_s to go to infinity (in which case the phase and the group velocities go to zero). This estimate is based on the shear layer scaling and may no longer apply close to s_c, since the wave beam is enlarged there, as seen in Fig. 16. We have, however, evaluated |$\mathcal {R}$| given by equation (A2) for all of the ‘background’ zonal flows in our linear simulations. In all cases, we find |$\mathcal {R}\gt 1/4$|⁠, leading us to two opposite conclusions:

• −either our zonal flows are stable under the |$\mathcal {R}$|-criterion, and the instabilities we have observed for ω = 0.8 and 0.9 have nothing to do with local shear instabilities, and are possibly more similar to the ‘Papaloizou–Pringle instability’.

• −the criterion on |$\mathcal {R}$| discussed here is no longer valid with viscous dissipation and tidal forcing, and must be modified.

Regarding our second hypothesis, the presence of tidal forcing presumably modifies equation (A1) close to s_c, based on the appendix of Astoul et al. [2021, equation (A.9) is the equivalent ODE for forced inertial waves]. This term could thus modify the expression for |$\mathcal {R}$| in equation (A2), and hence also change the stability criterion.