© The Authors 2024

1 Introduction

Tidal interactions between Jupiter and its Galilean satellites are recognised as influential factors in both the system’s orbital evolution and the internal dynamics (e.g. Lainey et al. 2009). Traditionally, the tidal response of a gaseous star or planet such as Jupiter is treated using the concept of equilibrium tide (Zahn 1966a, 1989; Remus et al. 2012), where hydrostatic deformation exhibits a minor phase lag in response to the dissipative processes caused by tidal forcing. However, the observed strong tidal dissipation in Jupiter (Lainey et al. 2009) and the gravitational perturbations recently measured by the Juno spacecraft (Durante et al. 2020) cannot be fully explained by the equilibrium tide alone. In fact, the Juno spacecraft has not only enhanced our understanding of Jupiter’s tidal dynamics, but it has also allowed us to delve deeper into the gravitational perturbation and tidal dissipation phenomena associated with the planet. On the one hand, it has acquired precise measurements of Jupiter’s tidal Love numbers, k_ℓm, which quantitatively characterise the planet’s response to tidal forcing represented by spherical harmonics of degree ℓ and order m. By analysing the real part of these Love numbers, we gain valuable information about the gravitational perturbations experienced by Jupiter.

On the other hand, the imaginary part of the Love numbers provides us with insights into the processes of tidal dissipation occurring within the planet. Recently, Durante et al. (2020) measured the Love number value for the dominant tidal component k₂₂ = 0.565 ± 0.018 (3σ uncertainty). This is lower than the theoretical hydrostatic value of $k_{22}^{(\text{eq})}=0.589$ as stated by Wahl et al. (2020), indicating a difference of approximately ∆k₂₂ ≈ −4%. Wahl et al. (2020) noted that the influence of the interior structure on k_ℓm is negligible when considering models that accurately reproduce the zonal harmonics J₂, J₄, and J₆, which have already been measured with high precision by Juno. This discrepancy between the observed and the computed values of k₂₂ could potentially be attributed to the influence of dynamical tides (Zahn 1975; Ogilvie & Lin 2004). Indeed, the conventional concept of the equilibrium tide does not satisfy the full equation of motion because the acceleration of the fluid is neglected (Zahn 1966b). Hence, a comprehensive understanding of the planet’s tidal response requires the inclusion of corrections. These corrections introduce wave-like motions within the planet and depend on both the tidal frequency and the internal structure (Ogilvie 2014). The dynamical (wave-like) tide offers additional channels for tidal dissipation and generates supplementary gravitational perturbations, surpassing the effects solely attributed to the hydrostatic deformation (Idini & Stevenson 2021; Lai 2021; Lin 2023).

The detection of gravitational signatures from dynamical tides can provide valuable insights into the interior structure of Jupiter, as it is influenced by both the tidal frequency and the internal structure. Idini & Stevenson (2021), Lai (2021), Idini & Stevenson (2022a,b), Dewberry & Lai (2022), Dewberry (2023), Lin (2023) found that ∆k₂₂ can be largely attributed to the Coriolis effect on the fundamental modes (f-modes). Additionally, Idini & Stevenson (2022a) proposed that resonant locking with a gravity mode in an extended diluted core could explain a ∆k₄₂ ≈ −11% difference between the observed and computed values of k₄₂. This finding supports the existence of a diluted core in Jupiter, a possibility also suggested by Stevenson (1985) and by Juno’s measurements of gravitational moments (Wahl et al. 2017; Militzer et al. 2022). The Coriolis force plays a crucial role in Jupiter’s tidal responses because tidal frequencies of its Galilean satellites are comparable with the planet’s spin frequency. Including the Coriolis force introduces inertial waves in (magnetised) convective regions (Rieutord & Valdettaro 1997, 2010; Ogilvie & Lin 2004; Ogilvie 2009, 2013; Baruteau & Rieutord 2013; Guenel et al. 2016a,b; Mathis et al. 2016; Wei 2016, 2018; Lin & Ogilvie 2018; Astoul et al. 2019) and a combination of gravity waves and inertial waves (known as gravito-inertial waves) in stably stratified regions (Dintrans et al. 1999; Dintrans & Rieutord 2000; Mathis 2009; Auclair Desrotour et al. 2015; André et al. 2019; Pontin et al. 2023).

The dissipation of dynamical tides occurs through various friction mechanisms, including turbulent friction in convective layers and heat diffusion in stably stratified regions (e.g. Ogilvie 2014; Mathis 2019; Duguid et al. 2020; Vidal & Barker 2020; de Vries et al. 2023). The rate of tidal dissipation in convective and stably stratified regions of planets has significant implications for the evolution of planet-moon systems. In the context of the Jupiter and Saturn systems, our understanding of tidal evolution has undergone a remarkable transformation. Both planets exhibit tidal dissipation that is one or several orders of magnitude stronger than previous predictions based on moon-formation scenarios (Goldreich & Soter 1966). This intense dissipation is essential to explaining their rapid orbital migration, a phenomenon that came to light through precise astrometric measurements (Lainey et al. 2009, 2012, 2017, 2020). For instance, Lainey et al. (2009) fitted a dynamical model, including parameterised tidal dissipation, to astrometric observations from 1891 to 2007 of the Galilean satellites. They found that the tidal dissipation is k₂₂/Q = (1.1 ± 0.2) × 10⁻⁵ (where Q is the quality factor that evaluates the ratio between the maximum energy stored in the tidal distortion and the energy dissipated during an orbital period), for the asynchronous tide due to Io.

Giant gas planets have traditionally been modelled with a three-layer model. This model entails a central rocky or icy core enveloped by a convective layer comprising metallic hydrogen and helium, which is further encompassed by an outer layer consisting of molecular hydrogen and helium (Stevenson 1982; Guillot et al. 1994). While this model serves as a reference for gas giant planets, uncertainties persist regarding the specific size of each region and the precise characteristics of the transitions between them. Recent studies have been diverging from the conventional standard model and delving into alternative interior structures. Specifically, Leconte & Chabrier (2012) proposed a model with a gradient of entropy and heavy elements throughout an entirely semi-convective (convective regions that are well mixed and separated by thin interfaces with stable stratification, creating a staircase-like structure in the entropy profile) planet, suppressing the need for a compact core in Jupiter. Stevenson (1985), Wahl et al. (2017), Debras & Chabrier (2019) investigated the possibility of incorporating stable stratification into their models, which takes the form of a substantial but diffuse core that extends beyond what was previously believed to be the convective zone. Therefore, significant portions of giant planet interiors are expected to exhibit an unstable entropy gradient, which competes with stable composition gradients. This competition can result in the emergence of double-diffusive convection, leading to the formation of semi-convective layers (e.g. Garaud 2018). Although these layers cannot be directly observed, their formation introduces distinct physics compared to traditional adiabatic models. Consequently, they have a profound impact on the behaviour and subsequent evolution of the system (e.g. Debras & Chabrier 2019). We thus go beyond the classical models of tides that invoke inertial waves in the deep convective envelope (Ogilvie & Lin 2004), viscoelastic tides in the rocky or icy core (Remus et al. 2012, 2015), or the combination of both (Guenel et al. 2014), and instead move towards models that consider gravito-inertial waves propagating in giant planets’ interiors, where both convective and stably stratified layers co-exist (André et al. 2017, 2019; Pontin et al. 2020, 2023; Lin 2023; Dewberry 2023).

In this study, we developed a method to calculate the dissipation of the dynamical tidal response of a self-gravitating, rotating fluid body composed of alternating convective layers and stably stratified layers and which takes into account the viscous, thermal, and chemical dissipation processes. This is the first time that global models incorporate the consideration of all three dissipation mechanisms, as opposed to solely focusing on viscosity in previous models. The Coriolis force is fully taken into account, but the centrifugal force is neglected as a first step. This method allows us to compute the imaginary part of the tidal Love number for a given planetary interior model. We focused especially on the latest Jupiter interior model constrained by Juno data and calculated by Debras & Chabrier (2019).

The paper is structured as follows. In Sect. 2, we describe how we derived the model that allows us to study the dissipation of tidally forced waves. We provide a detailed explanation for the separation of equilibrium tides and dynamical tides. Additionally, we derived the energy equation of tidal flows within this framework. In Sect. 3, we present the Jupiter interior model used in this study. In Sect. 4, we outline our methodology for calculating the equilibrium tide in the adiabatic case. Then, we focus on computing the dynamical (wave-like) tide from 2D linear pseudo-spectral numerical simulations, which allowed us to derive the associated dissipation. We then present the novel results obtained concerning tidal waves in Jupiter in Sect. 5. Specifically, we cover the simultaneous inclusion of inertial waves in convection zones and gravito-inertial waves in stably stratified zones along with the evaluation of the dissipation resulting from the different dissipative processes. Finally, in Sect. 6, we summarise the key findings and implications of our study.

2 Modelling tidally forced waves in giant planet interiors

We studied the linear excitation of (gravito-)inertial waves using an external tidal body forcing F_forcing. These waves are subject to dissipative processes, namely viscosity, thermal diffusion, and molecular diffusion (we assume that these diffusivities are uniform (cf. Sect. 3.2)).

2.1 Governing equations

We begin by writing the system of dynamical equations formed by the following set of equations. First, we write the continuity equation: $$D_t\rho +\rho \nabla \cdot V=0,$$(1)

where ρ is the density, V is the velocity field, and D_t = ∂_t + (V⋅∇) is the Lagrangian derivative. Then, we introduce the momentum equation: $$\rho \frac{DV}{Dt}=-\nabla P-\rho \nabla \text{Φ}+\rho v\left({\nabla }^{2}V+\frac{1}{3}\nabla \nabla \cdot V\right)+F_{\text{forcing }},$$(2)

where P is the pressure, Φ the gravitational potential, v is the kinematic viscosity assumed to be constant, and F_forcing the tidal forcing. We adopt here the Stokes hypothesis, where the bulk viscosity is neglected. We also introduce the heat (energy) equation: $$\rho T{D}_{t}s=k{\nabla }^{2}T,$$(3)

where T is the temperature, k is the thermal conductivity, and s is the specific entropy such as $\text{d}s:=c_p\left(\frac{\text{d}T}{T}-{\nabla }_{\text{ad}}\frac{\text{d}P}{P}\right)$; with ${\nabla }_{\text{ad}}:={\left(\frac{\text{d}\ln T}{ \text{d}\ln P}\right)}_S$ being the adiabatic temperature gradient and c_p the specific heat capacity. We neglect the viscous heating term here and suppose k is constant. The chemical composition equation can be written as $$D_t\mu =D_{\mu }{\nabla }^2\mu ,$$(4)

where µ is the molecular weight and D_µ is the molecular diffusion supposed constant. The Poisson equation reads $${\nabla }^2\text{Φ}=4\pi G\rho ,$$(5)

where G the universal gravitational constant. Finally, the general differential form of the equation of state (Kippenhahn & Weigert 1994) is defined by $$\frac{\partial \rho }{\rho }=\alpha \frac{\partial P}{P}-\delta \frac{\partial T}{T}+\varphi \frac{\partial \mu }{\mu },$$(6)

with $$\alpha :={\left(\frac{\partial \ln \rho }{\partial \ln P}\right)}_{T,\mu },\delta :=-{\left(\frac{\partial \ln \rho }{\partial \ln T}\right)}_{P,\mu },\varphi :={\left(\frac{\partial \ln \rho }{\partial \ln \mu }\right)}_{P,T}.$$(7)

2.2 Linearisation

We linearise the hydrodynamic system (Eqs. (1)–(6)) around the hydrostatic steady state. Each scalar field X := {P,ρ, Φ, T, µ} is expanded as the sum of its hydrostatic value X₀ and of the Eulerian perturbations associated with the tides X′: $$X(r,\theta ,\phi ,t)=X_0(r)+X\prime (r,\theta ,\phi ,t).$$(8)

We neglect the non-spherical character of the hydrostatic background here due to the deformation associated with the centrifugal acceleration, and the associated perturbation of the gravitational potential since gravito-inertial waves are only slightly affected by the deformation (e.g. Ballot et al. 2010; Dhouib et al. 2021). This implies that the background is independent of θ, so X₀ = X₀(r). We can write the velocity field, V, as the sum of the large-scale azimuthal velocity associated with the uniform rotation (as a first step, we neglect the differential rotation since Jupiter’s relative differential rotation is 4%; Guillot et al. 2018), Ω, and of the wave velocity, υ: $$V(r,\theta ,\phi ,t)=r\sin \theta \text{\hspace{0.17em}Ω\hspace{0.17em}\hspace{0.17em}}{e}_{\phi }+\upsilon (r,\theta ,\phi ,t),$$(9)

where t is time and (r, θ, φ) are the usual spherical coordinates with their associated unit vector basis (e_r, e_θ, e_φ). In this case, the linearised system ((1)–(6)) can be rewritten in the rotating frame as $$\frac{\partial {\rho }^{\prime }}{\partial t}+\nabla \cdot \left({\rho }_0\upsilon \right)=0,$$(10) $$\begin{array}{l}\frac{\partial \upsilon }{\partial t}+2\text{Ω}\times \upsilon =-\nabla W-\frac{\nabla {\rho }_0}{{\rho }_0}{W}^{\prime }+\frac{{\rho }^{\prime }}{{\rho }_0}{g}_0-\nabla {\text{Φ}}^{\prime }\\ +v\left({\nabla }^2\upsilon +\frac{1}{3}\nabla \nabla \cdot \upsilon \right)+F_{\text{forcing }},\end{array}$$(11) $$\frac{\partial T^{\prime }}{\partial t}+\upsilon \cdot \nabla T_0-{\nabla }_{\text{ad}}\frac{T_0}{P_0}\left(\frac{\partial P^{\prime }}{\partial t}+\upsilon \cdot \nabla P_0\right)=\kappa {\nabla }^2{T}^{\prime },$$(12) $$\frac{\partial {\mu }^{\prime }}{\partial t}+\upsilon \cdot \nabla {\mu }_0=D_{\mu }{\nabla }^2{\mu }^{\prime },$$(13) $${\nabla }^2{\text{Φ}}^{\prime }=4\pi G{\rho }^{\prime },$$(14) $$\frac{{\rho }^{\prime }}{{\rho }_0}=\alpha \frac{P^{\prime }}{P_0}-\delta \frac{T^{\prime }}{T_0}+\varphi \frac{{\mu }^{\prime }}{{\mu }_0},$$(15)

where W′ = P′/ρ₀ is the normalised pressure, κ = k/ρ₀c_p is the thermal diffusivity supposed constant, and 𝑔₀ = −∇Φ₀ = ∇P₀/ρ₀ = −𝑔₀e_r is the gravitational acceleration.

2.3 Non-wave-like and wave-like tides

We decompose the fluctuations associated with the tides into non-wave-like and wave-like parts: $$Y=Y^{\text{nw}}+Y^{\text{w}},$$(16)

with Y := {υ_r, υ_θ, υ_φ, X′}, where Y^nw is the non-wave-like (equilibrium) tide that satisfies the hydrostatic equilibrium (Zahn 1966a, 1989) and Y^w the wave-like (dynamical) tide that describes the propagation of waves (Zahn 1975; Ogilvie & Lin 2004).

2.3.1 Non-wave-like part

We assume that the non-wave-like part is adiabatic (α ≈ 1/Γ_ı and κ = 0), where Γ₁ = (∂ ln P₀/∂ ln ρ₀)_s is the first adiabatic exponent, and non-dissipative (v = D_µ = 0). The planet is assumed to be tidally forced by a single potential component, F_forcing = −∇Ψ, where $$\text{Ψ}(r,\theta ,\phi ,t)={\text{Ψ}}_{\ell }(r)Y_{\ell }^m(\theta ,\phi ){\text{e}}^{-\text{i}{\omega }_{0}t}\text{, with }{\text{Ψ}}_{\ell }(r)=A{\left(\frac{r}{R}\right)}^{\ell }\text{, }$$(17)

where $$A\propto \frac{G{R}^{\ell }{M}_{\text{satellite }}}{a^{\ell +1}}$$(18)

is the forcing amplitude, M_sateilite is the mass of the satellite, and a is the semi-major axis. In our linear numerical calculations, we used a normalised value of A, so we set A = 1.

$Y_{\ell }^m(\theta ,\phi )$ is an orthonormalised spherical harmonic of degree ℓ and order m, and ω₀ = nΩ_orbital is the tidal frequency in the inertial frame centred on the planet (n labels temporal harmonics of the orbital motion and Ω_orbital denotes the orbital frequency). We only considered the dominant quadrupolar tidal component ℓ = m = 2.

If we suppose that the adiabatic equilibrium tide is stationary in the frame rotating with the fluid inside the planet (Remus et al. 2012), we simplify the linearised heat (Eq. (12)) and chemical composition (Eq. (13)) equations by neglecting ∂P^nw/∂t. In that case, the system of equations that describes the non-wave-like tides can be written as $${\rho }^{\text{nw}}+\nabla \cdot \left({\rho }_0{\xi }^{\text{nw}}\right)=0,$$(19) $$-\nabla W^{\text{nw}}-\frac{\nabla {\rho }_0}{{\rho }_0}{W}^{\text{nw}}+\frac{{\rho }^{\text{nw}}}{{\rho }_0}{g}_0-\nabla {\text{Φ}}^{\text{nw}}-\nabla \text{Ψ}=0,$$(20) $$T^{\text{nw}}+{\xi }_r^{\text{nw}}{T}_0\frac{N_{\text{t}}^2}{g_0\delta }=0$$(21) $${\mu }^{\text{nw}}-{\xi }_r^{\text{nw}}{\mu }_0\frac{N_{\mu }^2}{g_0\varphi }=0,$$(22) $${\nabla }^2{\text{Φ}}^{\text{nw}}=4\pi G{\rho }^{\text{nw}},$$(23) $$\frac{{\rho }^{\text{nw}}}{{\rho }_0}=\frac{1}{{\Gamma }_1}\frac{P^{\text{nw}}}{P_0}-\delta \frac{T^{\text{nw}}}{T_0}+\varphi \frac{{\mu }^{\text{nw}}}{{\mu }_0},$$(24)

where ξ^nw is the displacement defined as υ^nw = ∂ξ^nw/∂t (it is customary to consider the displacement ξ^nw instead of the velocity υ^nw in the calculation of the non-wave-like tide since it is a deformation induced by mass redistribution), $$N_{\text{t}}^2=-g_0\delta \frac{\text{d}\ln P_0}{ \text{d}r}\left({\nabla }_{\text{ad}}-\frac{\text{d}\ln T_0}{ \text{d}\ln P_0}\right)$$(25)

is the thermal Brunt-Väisälä frequency squared, and $$N_{\mu }^2=-g_0\varphi \frac{\text{d}\ln {\mu }_0}{ \text{d}r}$$(26)

is the compositional Brunt-Väisälä frequency squared. The sum of these two qualities gives us the total Brunt–Väisälä frequency squared: $$N^2=N_{\text{t}}^2+N_{\mu }^2=-g_0\left(\frac{1}{{\rho }_0}\frac{\text{d}{\rho }_0}{ \text{d}r}-\frac{1}{{\text{Γ}}_1{P}_0}\frac{\text{d}{P}_0}{ \text{d}r}\right).$$(27)

If we write Eq. (20) as $$-\nabla \left(W^{\text{nw}}+{\text{Φ}}^{\text{nw}}+\text{Ψ}\right)+\left({\rho }^{\text{nw}}-\frac{\text{d}{\rho }_0}{ \text{d}{P}_0}{P}^{\text{nw}}\right)\frac{\nabla P_0}{{\rho }_0^2}=0,$$(28)

we deduce that (Ogilvie 2014) $$P^{\text{nw}}=-{\rho }_0\left({\text{Φ}}^{\text{nw}}+\text{Ψ}\right),$$(29) $${\rho }^{\text{nw}}=\frac{\text{d}{\rho }_0}{ \text{d}{P}_0}{P}^{\text{nw}}.$$(30)

The non-wave-like gravitational potential is obtained by solving the Poisson equation (Eq. (23)), which can be rewritten as $$\frac{1}{r^2}\frac{\text{d}}{\text{d}r}\left(r^2\frac{\text{d}{\text{Φ}}_{\ell }^{\text{nw}}}{\text{d}r}\right)-\frac{\ell (\ell +1)}{r^2}{\text{Φ}}_{\ell }^{\text{nw}}+4\pi G\frac{\text{d}{\rho }_0}{ \text{d}{P}_0}{\rho }_0\left({\text{Φ}}_{\ell }^{\text{nw}}+{\text{Ψ}}_{\ell }\right)=0,$$(31)

where $${\text{Φ}}^{\text{nw}}(r,\theta ,\phi ,t)={\text{Φ}}_{\ell }^{\text{nw}}(r)Y_{\ell }^m(\theta ,\phi ){\text{e}}^{-\text{i}{\omega }_{0}t},$$(32)

and with the following boundary conditions to ensure its regularity near the centre (r = η, where η is the aspect ratio) and its continuity at the surface (r = 1) (e.g. Ogilvie 2009): $$\frac{\text{d}\ln {\text{Φ}}_{\ell }^{\text{nw}}}{\text{d}\ln r}=\ell \text{ at }r=\eta R,$$(33) $$\frac{\text{d}\ln {\text{Φ}}_{\ell }^{\text{nw}}}{\text{d}\ln r}=-(\ell +1)\text{ at }r=R.$$(34)

Using Eqs. (21) and (22) in Eq. (24) we obtain $$\frac{1}{{\text{Γ}}_1}\frac{P^{\text{nw}}}{P_0}-\frac{{\rho }^{\text{nw}}}{{\rho }_0}+\frac{N^2}{g_0}{\xi }_r^{\text{nw}}=0,$$(35)

where ${\xi }^{\text{nw}}={\xi }_r^{\text{nw}}{e}_r+{\xi }_h^{\text{nw}}$ such that ${\xi }_h^{\text{nw}}\cdot e_r=0$. From this equation we can derive the expression for the non-wave-like radial displacement, $$N^2\left({\xi }_{r,\ell }^{\text{nw}}+\frac{{\text{Φ}}_{\ell }^{\text{nw}}+{\text{Ψ}}_{\ell }}{g_0}\right)=0;$$(36)

so, when N² ≠ 0 we obtain $${\xi }_{r,\ell }^{\text{nw}}=-\frac{{\text{Φ}}_{\ell }^{\text{nw}}+{\text{Ψ}}_{\ell }}{g_0}.$$(37)

Subsequently, from Eq. (19) we can derive the non-wave-like horizontal displacement $${\xi }_{h,\ell }^{\text{nw}}=\frac{1}{\ell (\ell +1)}\left(2{\xi }_{r,\ell }^{\text{nw}}+r\frac{\text{d}{\xi }_{r,\ell }^{\text{nw}}}{\text{d}r}\right);$$(38)

therefore, we can deduce that $$\mathrm{div }{\xi }^{\text{nw}}=0$$(39)

This is the conventional equilibrium tide (Zahn 1966a, 1989; Remus et al. 2012). This solution applies not only to stably stratified zones, but also to convective regions, since N² is not equal to zero but slightly negative and the fact that we generally set N² in these zones to zero is only an approximation. Thus, Eqs. (37) and (38) may be applied to the whole fluid domain inside the planet. Terquem et al. (1998) and Goodman & Dickson (1998) argued that this equilibrium tide solution does not apply to con-vective regions, as they assumed that the convective zone is adiabatically stratified (N² = 0). A comparison between these two definitions of the non-wave-like tides performed by Barker (2020) highlights the fact that in the interface between convective zones and stably stratified zones, a discontinuity arises in the horizontal component of displacement. This situation poses a problem both physically, since the ellipsoidal deformation and the related displacement has no reason to be discontinuous at convective-radiative boundaries, and numerically when dealing with multi-zone problems. We can therefore use the solution derived by Zahn (1966a, 1989) and Remus et al. (2012) which applies in stably stratified zones, but also in convective regions, since N² is not strictly equal to zero, but slightly negative.

2.3.2 Wave-like part

To derive the wave-like part, we first assumed the Boussinesq approximation (Spiegel & Veronis 1960), which neglects the density variations except where they appear in the buoyancy term, so the acoustic waves are filtered out. This approximation is an essential first step for addressing such a complex problem where the eigenmodes at these low frequencies are generally singular and are regularised by diffusion processes (see Sect. 6 for the discussion on the use of the anelastic approximation instead of the Boussinesq one). In fact, calculating inertial and gravito-inertial waves in an internal structure model with multiple transition layers poses a challenge, particularly given the presence of the strong density gradients. Additionally, incorporating three diffusion processes with coefficients spanning several orders of magnitude, potentially reaching very low values, adds another layer of complexity that is numerically demanding. To manage these complexities effectively, it is necessary to start with a simplified model to control the physical processes before moving on to the following stages. This approach is crucial for acknowledging and addressing possible biases introduced during the analysis. In Sect. 6, we discuss the potential limitations of this first necessary study within the Boussinesq approximation and the requirements to go beyond it in a near future. Then, we assumed the Cowling approximation (Cowling 1941), which neglects the perturbations of the gravitational potential induced by the waves since the perturbations induced by the non-wave-like tides are dominant (e.g. Ogilvie & Lin 2004). The system of equations that describes the wave-like tides can thus be written as $$\nabla \cdot {\upsilon }^{\text{w}}=0,$$(40) $$\frac{\partial {\upsilon }^{\text{w}}}{\partial t}+2\text{Ω}\times {\upsilon }^{\text{w}}=-\nabla W^{\text{w}}+g_0\left(\delta \frac{T^{\text{w}}}{T_0}-\varphi \frac{{\mu }^{\text{w}}}{{\mu }_0}\right)e_r+v{\nabla }^2{\upsilon }^{\text{w}}+f^{\text{nw}},$$(41) $$\frac{\partial T^{\text{w}}}{\partial t}+{\text{υ}}_r^{\text{w}}{T}_0\frac{N_{\text{t}}^2}{g_0\delta }=\kappa {\nabla }^2{T}^{\text{w}},$$(42) $$\frac{\partial {\mu }^{\text{w}}}{\partial t}-{\text{υ}}_r^{\text{w}}{\mu }_0\frac{N_{\mu }^2}{g_0\varphi }=D_{\mu }{\nabla }^2{\mu }^{\text{w}},$$(43) $$\frac{{\rho }^{\text{w}}}{{\rho }_0}=-\delta \frac{T^{\text{w}}}{T_0}+\varphi \frac{{\mu }^{\text{w}}}{{\mu }_0},$$(44)

with $$f^{\text{nw}}=-\frac{\partial {\upsilon }^{\text{nw}}}{\partial t}-2\text{Ω}\times {\upsilon }^{\text{nw}},$$(45)

the forcing term that arises when solving the non-wave-like tides as a residual force, as the non-wave-like one does not satisfy the equation of motion due to the omission of inertial forces associated with this flow. This force encompasses the acceleration of the non-wave-like tide and the Coriolis acceleration applied to it and will force the gravito-inertial tidal waves (see also Ogilvie 2005; André et al. 2019).

We choose the planet’s radius R for the length scale and (2Ω)⁻¹ for the timescale (t = τ/2Ω). Therefore, we can define the normalised quantities as follows: υ^w = 2ΩRu, W^w = 4Ω²R²π, $g_0=4{\text{Ω}}^{2}R{g}_0^{*}$, ƒ^nw = 4Ω²Rƒ^*, and the normalised differential operator as R∇ = ∇*. We write the normalised temperature and chemical composition as T^W = T_iΘ, µ^w = µ_iM, $T_0=T_{\text{i}}{T}_0^{*}$ and ${\mu }_0={\mu }_{\text{i}}{\mu }_0^{*}$, where T_i and µ_i are the temperature and the molecular weight, respectively, at the inner boundary. So, we can write the normalised system as $${\nabla }^{*}\cdot u=0,$$(46) $$\frac{\partial u}{\partial \tau }+e_z\times u=-\nabla \Pi +\left(\frac{\delta }{T_0^{*}}\text{Θ}-\frac{\varphi }{{\mu }_0^{*}}M\right)g_0^{*}{e}_r+{\text{E∇}}^{*2}u+f^{*},$$(47) $$\frac{\partial \text{Θ}}{\partial \tau }+\frac{T_0^{*}{N}_{\text{t}}^{*2}}{g_0^{*}\delta }{u}_r=\frac{\text{E}}{\mathrm{Pr}}{\nabla }^{*2}\text{Θ},$$(48) $$\frac{\partial M}{\partial \tau }-\frac{{\mu }_0^{*}{N}_{\mu }^{*2}}{g_0^{*}\varphi }{u}_r=\frac{\text{E}}{\text{Sc}}{\nabla }^{*2}M,$$(49)

where we define the normalised Brunt–Väisälä frequencies: $$N_{\text{t}}^{*2}=\frac{N_{\text{t}}^2}{4{\text{Ω}}^2}\text{ and }{N}_{\mu }^{*2}=\frac{N_{\mu }^2}{4{\text{Ω}}^2}\text{. }$$(50)

These equations are governed by three dimensionless numbers. These are the Prandtl number, defined as the ratio of the kinematic viscosity (v) to the thermal diffusivity (κ), $$\mathrm{Pr}=\frac{v}{\kappa };$$(51)

the Schmidt number, defined as the ratio of the kinematic viscosity (v) to the molecular diffusivity (D_µ), $$\text{Sc}=\frac{v}{D_{\mu }};$$(52)

and the Ekman number, which compares the ratio between the viscous force and the Coriolis force, $$\text{E}=\frac{v}{2\text{Ω}{R}^2}.$$(53)

2.4 Energy equation

From the scalar product between $\overline{u}$ (where $\overline{\square }$ denotes the complex conjugate) and the momentum Eq. (47) and by using Eqs. (48) and (49), we obtain the energy equation $${\partial }_{\tau }\left(E_{\text{k}}+E_{\text{p},\text{th}}+E_{\text{p},\mu }\right)=P_{\text{acou}}+D_{\text{th}}+D_{\text{ch}}+D_{\text{visc}}+P_{\text{tide}},$$(54)

with $$E_{\text{k}}=\frac{1}{2}\Vert u{\Vert }^2$$(55)

being the specific kinetic energy, $$E_{\text{p},\text{th}}=\{\begin{array}{ll}\frac{1}{2}{\left(\frac{\delta g_0^{*}}{T_0^{*}{N}_{\text{t}}^{*}}\right)}^2{\left|\text{Θ}\right|}^2\hfill & \text{ if }{N}_{\text{t}}^2\ne 0\hfill \\ 0\hfill & \text{ if not }\hfill \end{array}$$(56)

the specific potential energy associated with thermal stratification, $$E_{\text{p},\mu }=\{\begin{array}{ll}\frac{1}{2}{\left(\frac{\varphi g_0^{*}}{{\mu }_0^{*}{N}_{\mu }^{*}}\right)}^2{\left|M\right|}^2\hfill & \text{ if }{N}_{\mu }^2\ne 0\hfill \\ 0\hfill & \text{ if not }\hfill \end{array}$$(57)

the specific potential energy associated with chemical stratification, $$P_{\text{acou }}=-{\nabla }^{*}\text{Π}\cdot \overline{u}$$(58)

the specific work of pressure forces that can be related to the acoustic flux ∇⋅(P′u), $$D_{\text{th}}=\{\begin{array}{ll}\frac{\text{E}}{\mathrm{Pr}}{\left(\frac{\delta g_0^{*}}{T_0^{*}{N}_{\text{t}}^{*}}\right)}^2\overline{\text{Θ}}{\nabla }^{*2}\text{Θ}\hfill & \text{ if }{N}_{\text{t}}^2\ne 0\hfill \\ 0\hfill & \text{ if not }\hfill \end{array}$$(59)

the specific power dissipated by thermal diffusion, $$D_{\text{ch}}=\{\begin{array}{ll}\frac{\text{E}}{\text{Sc}}{\left(\frac{\varphi g_0^{*}}{{\mu }_0^{*}{N}_{\mu }^{*}}\right)}^2\overline{M}{\nabla }^{*2}M\hfill & \text{ if }{N}_{\mu }^2\ne 0\hfill \\ 0\hfill & \text{ if not }\hfill \end{array}$$(60)

the specific power dissipated by chemical diffusion, $$D_{\text{visc }}=E{\nabla }^{*2}u\cdot \overline{u}$$(61)

the specific power dissipated by viscous friction, and $$P_{\text{tide }}=f^{*}\cdot \overline{u}$$(62)

the specific tidal power. Then, after spatial integration over the volume V, we obtain $${\partial }_t\left({\tilde{E}}_{\text{k}}+{\tilde{E}}_{\text{p},\text{th}}+{\tilde{E}}_{\text{p},\mu }\right)={\tilde{P}}_{\text{acou }}+{\tilde{D}}_{\text{th }}+{\tilde{D}}_{\text{ch}}+{\tilde{D}}_{\text{visc }}+{\tilde{P}}_{\text{tide }},$$(63)

where $\tilde{\square }=\frac{1}{2}\Re \left[{\displaystyle {\int }_V\square }\text{d}V\right]$. Here, we assume that the density is constant by adopting the Boussinesq approximation, allowing us to simplify it in this equation on both sides.

3 Interior model of Jupiter as revealed by Juno

3.1 Five-layer model

In order to assess the dissipation of dynamical tides on Jupiter, we need to prescribe the background state profiles of the different quantities and the Brunt–Väisälä frequencies. Recent observations have made significant advancements in our understanding of the internal structure of gas giant plants, such as Jupiter and Saturn, yet some degree of uncertainty remains. As we can see in Fig. 1, the internal structure model computed by Debras & Chabrier (2019) to reproduce Jupiter’s multipolar moments as measured by Juno assumes an extended diluted core of radius 0.69R treated as a stably stratified fluid layer (a similar zone is probably also present in Saturn; Mankovich & Fuller 2021), and a convective envelope that features a small stably stratified layer between 0.9R and 0.92R, possibly resulting from H-He immis-cibility (Debras & Chabrier 2019). The equation of state used to compute this model is the one derived in Chabrier & Debras (2021). As illustrated in Fig. 2, starting from its surface and moving towards the core, Jupiter is thought to exhibit the following layers:

• Gaseous envelope: this outermost layer is characterised by convective motion and differential rotation.

• Transitional stably stratified zone: this region is considered to be potentially semi-convective. It is also known as a double-diffusive zone, as proposed by Leconte & Chabrier (2012).

• Internal convective zone: this layer is magnetised and composed of metallic hydrogen and helium, and it is rotating quasi-uniformly since if the rotation were significantly different from solid rotation, the ohmic dissipation would become inexplicable (e.g. Liu et al. 2008; Guillot et al. 2018).

• Stably stratified zone: this layer is located closer to the core and may exhibit double diffusion convection or a diluted core structure due to stabiliing composition gradients (Leconte & Chabrier 2012, 2013; Wilson & Militzer 2012b,a; Wahl et al. 2013, 2017; González-Cataldo et al. 2014; Mazevet et al. 2015).

• Potential unstable solid core of size 1.4% of radius: made up of rock or ice. While this size may vary slightly based on different models, we are unable to construct a model with a substantial core (bigger than ~10%) that respects the Juno constraints (e.g. Debras & Chabrier 2019).

We note that as a first step, we did not take the differential rotation and magnetism into account in this study.

Fig. 1 Profiles of compositional $\left(N_{\mu }^2\right)$, thermal $\left(N_t^2\right)$, and total (N2) Brunt-Väisälä frequencies squared normalised by the inertial frequencies squared (4Ω2) as a function of the normalised radius r/R.

Fig. 2 Schematic of model of Jupiter’s interior used in this study.

3.2 Transport properties

Using the transport properties outlined in Stevenson & Salpeter (1977), we made calculations to determine the various molecular diffusivities within Jupiter and the associated dimensionless numbers (the expressions of these numbers are given in Appendix A). The radial profile of these dimensionless numbers (Prandtl, Schmidt and Ekman) is displayed in Fig. 3. We can see that the influence of viscous forces is generally small compared to the Coriolis acceleration. Consequently, the Ekman number, which characterises the ratio of viscous to Coriolis forces, becomes extremely low, typically of the order of 10⁻¹⁷ when assuming molecular viscosity. However, due to the challenges in accurately resolving shear layers at very low diffusivity levels, such Ekman number regimes cannot be reached numerically. The numerical strategy adopted is therefore to reach the lowest possible Ekman number values, hoping to have reached a regime where the scaling laws obtained will apply to the lower astrophysical values. Moreover, this Ekman number value is calculated with a molecular viscosity value, whereas inertial tidal waves could in reality be subject to a turbulent effective viscosity (Ogilvie & Lin 2004, 2007; Mathis et al. 2016; Duguid et al. 2020; Vidal & Barker 2020; de Vries et al. 2023), with larger values that would lead to a larger Ekman number. Indeed, by employing the non-rotating mixing-length theory, we can make a rough estimation of the turbulent effective eddy viscosity in convective regions resulting in an Ekman number of approximately 10^-7 (Guillot et al. 2004). When replacing the standard non-rotating mixing length theory by the rotating mixing length theory developed in Stevenson (1979), we end up (following Mathis et al. 2016) with a much smaller turbulent Ekman number of approximatively 10⁻¹⁵, which is closer to the microscopic value, because of the inhibition of convection by rapid rotation (Fuentes et al. 2023). On the other hand, the Prandtl number, which measures the ratio of viscosity to thermal diffusivity, is low within the planet (approximately 10⁻²), but it increases to around 1 near the surface (at r > 0.9R). Regarding the Schmidt number, which characterises the ratio of viscosity to molecular diffusivity, it remains close to unity throughout the planet. As a first step, we assumed in our simulations that these dimensionless numbers are constant. Then, we were able to study the impact of their variations by exploring the parameter domain.

Fig. 3 Prandtl, Schmidt, and Ekman number profiles as a function of the normalised radius based on Stevenson & Salpeter (1977).

4 Numerical resolution

Our attention in this paper is directed towards the ℓ = m = 2 component of the tide, as it is commonly considered to be the most prominent for a quasi-circular and coplanar two-body system (Mathis & Le Poncin-Lafitte 2009; Ogilvie 2014).

4.1 Solving the non-wave-like part

Our goal here is to find the non-wave-like displacement ξ^nw , as it is needed in the expression for the effective forcing driving dynamical tides (Eq. (45)). Therefore, we must first solve Poisson’s equation (Eq. (31)) numerically with boundary conditions (33) and (34) using the density and pressure background profiles (ρ₀ and P₀) that we computed using the structure model defined in Sect. (3.1). Then, we were able to easily compute the vertical and horizontal non-wave-like displacement using Eqs. (37) and (38). We can see, in Fig. 4, these quantities as a function of the normalised radius for A = 1 and ℓ = 2. We emphasise here that setting the forcing value to A = 1 gives very high displacement values, while a realistic forcing value (Eq. (18)) would give much lower values. Since we are dealing with the linear case, the choice of A has no impact on our results (ultimately, we want to calculate the Love number defined just after in Eq. (64), which is a ratio where A will be simplified).

Once the solution to Poisson’s equation (Eq. (31)) is obtained numerically, the Love number (Love 1911) is readily given by $$k_{\ell }^{\text{nw}}=\frac{{\text{Φ}}_{\ell }^{\text{nw}}(r=R)}{{\text{Ψ}}_{\ell }(r=R)}.$$(64)

We note that this includes only non-wave-like tides, and it is real since we do not take into account its dissipation. For ℓ = 2 we find that $k_2^{\text{nw}}=0.638$. However, it is important to acknowledge that our assumption of a spherical planet does not hold true for Jupiter, as its rapid rotation causes it to be flattened, leading to a discrepancy between our value and the calculated values by Wahl et al. (2020). In fact, they investigated the nonrotating case and found a value of $k_2^{\text{nw}}=0.536$. However, when taking the planet’s rotation into account, they determined that the value changes to $k_2^{\text{nw}}=0.589$ for m = 2. It is worth noting that the Love number becomes dependent on the azimuthal order m when considering the effects of rotation. Our approach takes into account rotation while neglecting flattening, whereas the method of Wahl et al. (2020) takes into account both rotation and induced deformation simultaneously. Since our main objective in this study is to understand wave-like tides, we computed the non-wave-like tides only to calculate the forcing term. Therefore, we can omit the flattening as a first step, considering that the forcing term will be only slightly modified. The discrepancy of the calculated equilibrium Love number in comparison to the observed value (k₂₂ = 0.565, Durante et al. 2020) suggests uncharacterised dynamical (wave-like) contribution due to tidal waves.

Fig. 4 Non-wave-like (equilibrium) gravitational potential (left) and corresponding radial and horizontal displacements (right) as a function of the normalised radius for A = 1 and ℓ = 2.

4.2 Solving the wave-like part

We expanded the velocity, temperature, molecular weight, and reduced pressure on spherical harmonics (Rieutord 1987) as $$\begin{array}{l}u={\displaystyle \sum _{\ell =0}^{+\infty }{\displaystyle \sum _{m=-\ell }^{\ell }(u_{\ell }^m(r)R_{\ell }^m(\theta ,\phi )+v_{\ell }^m(r)S_{\ell }^m(\theta ,\phi )}}\\ +w_{\ell }^m(r)T_{\ell }^m(\theta ,\phi )){\text{e}}^{-i\tilde{\omega }\tau },\end{array}$$(65) $$\text{Θ}={\displaystyle \sum _{\ell =0}^{+\infty }{\displaystyle \sum _{m=-\ell }^{\ell }{t}_{\ell }^m}}(r)Y_{\ell }^m(\theta ,\phi ){\text{e}}^{-i\tilde{\omega }\tau },$$(66) $$M={\displaystyle \sum _{\ell =0}^{+\infty }{\displaystyle \sum _{m=-\ell }^{\ell }{\mu }_{\ell }^m}}(r)Y_{\ell }^m(\theta ,\phi ){\text{e}}^{-i\tilde{\omega }\tau },$$(67) $$\text{Π}={\displaystyle \sum _{\ell =0}^{+\infty }{\displaystyle \sum _{m=-\ell }^{\ell }{p}_{\ell }^m}}(r)Y_{\ell }^m(\theta ,\phi ){\text{e}}^{-i\tilde{\omega }\tau },$$(68)

with $R_{\ell }^m=Y_{\ell }^m{e}_r,S_{\ell }^m=\nabla Y_{\ell }^m,T_{\ell }^m=\nabla \times R_{\ell }^m$, and $\tilde{\omega }=\omega /2\text{Ω}$ being the normalised frequency associated with the tidal forcing in the rotating frame. Then, the projection of the linearised dimensionless system (Sect. 4.2.1) and the associated boundary conditions (Sect. 4.2.2) are solved using the linear 2D pseudo-spectral Linear Solver Builder (LSB) code (Valdettarov 2007). These equations are discretised in the radial direction on the Gauss-Lobatto collocation nodes associated with Chebyshev polynomials. They are truncated at order N_r for the Cheby-shev basis and at order N_ℓ for the spherical harmonics basis. The governing equations, Eqs. (69) to (74), and adopted boundary conditions (Sect. 4.2.2) form a linear system of the form ℳX = F, which is solved on each point of the radial grid, given the value of the forcing frequency $\tilde{\omega }$ and azimuthal order m.

4.2.1 System of equations to solve

Using Eqs. (65)–(68) and the expressions of the operators in the spherical harmonics basis specified in Appendix B, we can rewrite the system (46)–(49) as $${\text{d}}_r{u}_{\ell }^m+\frac{2}{r}{u}_{\ell }^m-\ell (\ell +1)\frac{v_{\ell }^m}{r}=0,$$(69) $$\begin{array}{l}\text{E}{\Delta }_{\ell }{u}_{\ell }^m-\left(\frac{2\text{E}}{r^2}-i\tilde{\omega }\right)u_{\ell }^m+\left(im+\frac{2\text{E}}{r^2}\ell (\ell +1)\right)v_{\ell }^m\hfill \\ -{\beta }_{\ell -1}^{\ell }{w}_{\ell -1}^m-{\beta }_{\ell +1}^{\ell }{w}_{\ell +1}^m-{\text{d}}_r{p}_{\ell }^m+\frac{\delta g_0^{*}}{T_0^{*}}{t}_{\ell }^m-\frac{\varphi g_0^{*}}{{\mu }_0^{*}}{\mu }_{\ell }^m=-f_{\text{R}}^{\ell ,m},\hfill \end{array}$$(70) $$\begin{array}{l}\text{E}{\Delta }_{\ell }{v}_{\ell }^m+\left(i\tilde{\omega }+\frac{im}{\ell (\ell +1)}\right)v_{\ell }^m+\left(\frac{2\text{E}}{r^2}+\frac{im}{\ell (\ell +1)}\right)u_{\ell }^m\\ -{\gamma }_{\ell -1}^{\ell }{w}_{\ell -1}^m-{\gamma }_{\ell +1}^{\ell }{w}_{\ell +1}^m-\frac{p_{\ell }^m}{r}=-f_{\text{S}}^{\ell ,m},\end{array}$$(71) $$\begin{array}{l}\text{E}{\Delta }_{\ell }{w}_{\ell }^m+\left(i\tilde{\omega }+\frac{im}{\ell (\ell +1)}\right)w_{\ell }^m+{\gamma }_{\ell -1}^{\ell }{v}_{\ell -1}^m+{\gamma }_{\ell +1}^{\ell }{v}_{\ell +1}^m\\ -\frac{{\alpha }_{\ell -1}^{\ell }}{\ell }{u}_{\ell -1}^m+\frac{{\alpha }_{\ell +1}^{\ell }}{\ell +1}{u}_{\ell +1}^m=-f_{\text{T}}^{\ell ,m},\end{array}$$(72) $$\frac{\text{E}}{\mathrm{Pr}}{\Delta }_{\ell }{t}_{\ell }^m-u_{\ell }^m\frac{T_0^{*}{N}_{\text{t}}^{*2}}{g_0^{*}\delta }+i\tilde{\omega }{t}_{\ell }^m=0,$$(73) $$\frac{\text{E}}{\text{Sc}}{\Delta }_{\ell }{\mu }_{\ell }^m+u_{\ell }^m\frac{{\mu }_0^{*}{N}_{\mu }^{*2}}{g_0^{*}\varphi }+i\tilde{\omega }{\mu }_{\ell }^m=0,$$(74)

where $${\text{d}}_r\square =\frac{\text{d}\square }{\text{d}r},$$(75) $${\text{d}}_{r^2}^2\square =\frac{{\text{d}}^2\square }{{\text{d}}^2},$$(76) $${\Delta }_{\ell }\square ={\text{d}}_{r^2}^2\square +\frac{2}{r} {\text{d}}_r\square -\frac{\ell (\ell +1)}{r^2}\square ,$$(77)

and the coupling coefficients, which all depend on m, are given by $${\alpha }_{\ell -1}^{\ell }={\alpha }_{\ell }^{\ell -1}=\sqrt{\frac{{\ell }^2-m^2}{(2\ell -1)(2\ell +1)}},$$(78) $${\beta }_{\ell -1}^{\ell }=(\ell -1){\alpha }_{\ell -1}^{\ell },{\beta }_{\ell +1}^{\ell }=-(\ell +2){\alpha }_{\ell +1}^{\ell },$$(79) $${\gamma }_{\ell -1}^{\ell }=\frac{\ell -1}{\ell }{\alpha }_{\ell -1}^{\ell },{\gamma }_{\ell +1}^{\ell }=\frac{\ell +2}{\ell +1}{\alpha }_{\ell +1}^{\ell }.$$(80)

Since we only considered the dominant quadrupolar tidal component ℓ = m = 2, the forcing term can be written as $$f^{*}=f_{\text{R}}^{\ell ,m}(r)R_{\ell }^m+f_{\text{S}}^{\ell ,m}(r)S_{\ell }^m+f_{\text{T}}^{\ell ,m}(r)T_{\ell }^m,$$(81)

with $$f_{\text{R}}^{\ell ,m}=\left[{\tilde{\omega }}^2{\xi }_{r,\ell }^{\text{nw}}+\frac{m\tilde{\omega }}{\ell (\ell +1)}\left(2{\xi }_{r,\ell }^{\text{nw}}+r {\text{d}}_r{\xi }_{r,\ell }^{\text{nw}}\right)\right]\delta (\ell -2)\delta (m-2),$$(82) $$\begin{array}{l}{f}_{\text{S}}^{\ell ,m}=[\frac{m\tilde{\omega }}{{\ell }^2{(\ell +1)}^2}\left((\ell (\ell +1)+2){\xi }_{r,\ell }^{\text{nw}}+r {\text{d}}_r{\xi }_{r,\ell }^{\text{nw}}\right)\hfill \\ +\frac{{\tilde{\omega }}^2}{\ell (\ell +1)}\left(2{\xi }_{r,\ell }^{\text{nw}}+r {\text{rd}}_{}{\xi }_{r,\ell }^{\text{nw}}\right)]\delta (\ell -2)\delta (m-2),\hfill \end{array}$$(83) $$f_{\text{T}}^{\ell ,m}=-i\tilde{\omega }\frac{{\alpha }_{\ell -1}^{\ell }}{{\ell }^2}\left(-(\ell -2){\xi }_{r,\ell -1}^{\text{nw}}+r {\text{d}}_r{\xi }_{r,\ell -1}^{\text{nw}}\right)\delta (\ell -3)\delta (m-2),$$(84)

obtained by projecting Eq. (45) on the spherical harmonics basis.

Fig. 5 Meridional cut (z = r cos θ and s = r sin θ) of dimensionless kinetic energy (left), potential energy associated with chemical (middle), and potential energy associated thermal stratification (right) of the forced mode $\tilde{\omega }=-0.76$ for m = 2, E = 10−7, Pr = Sc = 1 with a spatial resolution of (Nr, Nℓ) = (300, 301). The trajectories of characteristics are represented with white curves, while the surfaces at which they undergo reflections are shown with brown curves.

Fig. 6 Spectral content of velocity field components (u, v, w), pressure (p), temperature (t), and molecular weight (µ) of the forced mode $\tilde{\omega }=-0.76$ for m = 2, E = 10−7, Pr = Sc = 1. Chebyshev and spherical harmonics coefficients are shown in the top and bottom panels, respectively.

4.2.2 Boundary conditions

Given our specific emphasis on (gravito-)inertial modes while excluding surface gravity modes, we were able to adopt the classical boundary conditions established in pioneering studies by Dintrans et al. (1999), Dintrans & Rieutord (2000), Valdettaro et al. (2007), Ogilvie & Lin (2004, 2007), Ogilvie (2005, 2009), Rieutord & Valdettaro (2010). Namely, we employed impenetrable and stress-free boundary conditions while assuming that the spheres bounding the fluid domain can absorb any flux of heat or chemical elements while remaining at constant temperature and molecular weight. That is to say, the radial functions must satisfy the following inner (r = η) and outer (r = 1) boundary conditions: $$u_{\ell }^m=0,$$(85) $$\frac{\text{d}{v}_{\ell }^m}{ \text{d}r}+\frac{u_{\ell }^m-v_{\ell }^m}{r}=0,$$(86) $$\frac{\text{d}}{\text{d}r}\left(\frac{w_{\ell }^m}{r}\right)=0,$$(87) $$t_{\ell }^m=0,$$(88) $${\mu }_{\ell }^m=0.$$(89)

5 Results

Having established the framework for our numerical work, we now present the numerical results. We discuss the basic properties of the forced waves before considering how the dissipative properties depend on the model’s key parameters and the implications for astrophysical tidal evolution. Our calculations are focused on the frequency range of $-1<\tilde{\omega }<-0.5$, which is directly relevant to the tidal frequencies of the Galilean moons. It is worth noting that the negative tidal frequency indicates that the tidal forcing is retrograde in the co-rotating frame with the planet, based on our convention.

Fig. 7 Dissipation spectra for m = 2, Pr = Sc = 1, and E = 10−6 (top), E = 10−7 (middle), and E = 10−8 (bottom) with a spatial resolution (Nr, Nℓ) = (300, 301). The total dissipation is represented by black crosses, and its viscous, thermal, and molecular contributions are represented by green points, red stars, and blue plus signs, respectively. The magenta dashed line indicates the viscous (total) dissipation in the case of a purely convective interior $\left(N_{\mu }^2=N_{\text{t}}^2=0\right)$.

Fig. 8 Same as middle panel of Fig. 7, but for Sc = 0.5 (top) and Sc = 6 (bottom).

5.1 Energies

In this section, we present our numerical results for the different types of energy examined in our study, with a specific focus on the forced mode $\tilde{\omega }=-0.76$ (the frequency exited on Jupiter by Io for m = n = 2). We adopt typical values for various dimensionless numbers, namely E = 10⁻⁷, Pr = Sc = 1. Furthermore, we set the aspect ratio to η = 0.014 (because of the solid core of size 1.4% of radius; see Sect. 3.1). The left panel of Fig. 5 illustrates the spatial distribution of the kinetic energy of this mode in a meridional quarter-plane since it is symmetrical with respect to the equator. Notably, we can observe two distinct types of modes. Firstly, there are gravito-inertial modes that exist within the inner stably stratified regions (and theoretically within the thin, outer, stably stratified layer; however, it is very thin, so it is not clear here). As expected, we find equatorial trapping of sub-inertial $\left(\tilde{\omega }<1\right)$ gravito-inertial modes (e.g. Dintrans et al. 1999; Dintrans & Rieutord 2000; Mathis 2009). Secondly, we have inertial modes present in the two convective zones, which are separated by the thin, stably stratified layer. These modes exhibit multiple reflections at the boundaries of their propagation zones, following specific trajectories known as attractors (Maas & Lam 1995). We note that the attractor starting from the critical latitude (e.g. Rieutord et al. 2001; Rieutord & Valdettaro 2018) in the inner convective zone and reflected at the pole, at the surface, at the equator, and at the interface with the innermost stably stratified zone seems to appear as well, regardless of the thin, intermediate, stably stratified region. We also represent the trajectories of characteristics with white curves in this figure, while the surfaces at which they undergo reflections are depicted by brown curves. These paths of characteristics are calculated based on the second-order partial differential equation satisfied by the pressure perturbation in the inviscid and short-wavelength approximations (see Mirouh et al. 2016 for the detailed derivation). The inclusion of these curves provides a valuable means of understanding the solutions to non-dissipative problems and validation of numerical calculations. We find that the patterns formed by the characteristics are in very good agreement with the numerical calculation, especially in the inner, stably stratified layers. In convective zones, the paths of characteristics follow straight lines that maintain a constant angle relative to the rotation axis (z-axis) in order to respect the inertial wave dispersion relation. In contrast, stably stratified regions introduce a distinct behaviour where the characteristics become curved, owing to the distortion caused by the presence of the stable stratification.

The middle and right panels of Fig. 5 reveal that the chemical and thermal energies primarily concentrate within the stably stratified regions, as the thermal and chemical Brunt-Väisälä frequencies approach zero within convective zones. Consequently, at the interfaces of the convective zones, both potential energies have a finite transition to zero.

To ensure the numerical convergence for this mode, we employed a spatial resolution of (N_r, N_ℓ) = (300, 301). This convergence can be appreciated by inspecting Fig. 6, where we display the spectral content of the velocity field components, pressure, temperature, and molecular weight for the same forced mode. In the top panel, we show the maximum Chebyshev coefficients C_k as a function of the Chebyshev order k, selecting the highest value among all the spherical harmonics coefficients corresponding to a given k. Likewise, the bottom panel displays the maximum spherical harmonics coefficients C_ℓ as a function of the spherical harmonic degree , considering the highest value among all Chebyshev coefficients. This spatial resolution (N_r, N_ℓ) = (300, 301) has proved to be sufficient up to values of E ≈ 10⁻⁸.

Fig. 9 Same as middle panel of Fig. 7, but for Pr = 0.5 (top) and Pr = 2 (bottom).

Fig. 10 Modified tidal quality factor as function of tidal frequency for m = 2, E = 10−7, and Pr = Sc = 1. Vertical dotted lines indicate the tidal frequencies for the four Galilean Moons of Jupiter (from right to left: Io, Europa, Ganymede, Callisto). The magenta dashed line indicates the values of these quantities in the case of a purely convective interior $\left(N_{\mu }^2=N_{\text{t}}^2=0\right)$. The dash-dotted orange line marks the observed value of this quantity due to Io.

Fig. 11 Imaginary part of Love number as function of tidal frequency for m = 2, Pr = Sc = 1, and E = 10−6 (top), E = 10−7 (middle), and E = 10−8 (bottom). Vertical dotted lines indicate the tidal frequencies for the four Galilean Moons of Jupiter (from right to left: Io, Europa, Ganymede, Callisto). The magenta dashed line indicates the values of these quantities in the case of a purely convective interior $\left(N_{\mu }^2=N_{\text{t}}^2=0\right)$. The dash-dotted orange line marks the observed value of this quantity due to Io.

Fig. 12 Same as left panel of Fig. 5, but for an interior model without the thin, external, stably stratified layer (two-zone model).

5.2 Dissipation spectra

We analysed three distinct forms of dissipation: viscous, thermal, and chemical. The total dissipation is defined as follows: $$\mathcal{D}(\tilde{\omega })={\tilde{D}}_{\text{th }}+{\tilde{D}}_{\text{ch }}+{\tilde{D}}_{\text{visc. }}$$(90)

Figure 7 shows the viscous $\left({\tilde{D}}_{\text{visc }}\right)$, thermal $\left({\tilde{D}}_{\text{th}}\right)$, molecular $\left({\tilde{D}}_{\text{ch}}\right)$, and total (𝒟) dissipation rates integrated over the volume as a function of the normalised forcing frequency $\left(\tilde{\omega }\right)$ for m = 2, E = {10⁻⁶, 10⁻⁷, 10⁻⁸}, and Pr = Sc = 1. We observe a significant frequency dependence, indicating a strong relationship between dissipation and the forcing frequency. Moreover, our analysis reveals that the dominant mechanism contributing to dissipation is viscosity, surpassing both thermal and chemical dissipations in magnitude. We ensure that the total energy is conserved $\mathcal{D}\approx {\tilde{P}}_{\text{tide }}$ to a given degree of confidence (maximum relative error of 5%). We note that, given our boundary conditions (Sect. 4.2.2), P_acou ≈ 0.

Following the comparison method adopted by André et al. (2019) in Cartesian coordinates, we also computed the dissipation spectra for the old vision of Jupiter’s interior, where there is a single purely convective zone extending from r = η = 0.014 to r = 1. In this scenario, the only form of dissipation present is viscous dissipation, as thermal and chemical dissipation are negligible due to $N_{\mu }^2=N_{\text{t}}^2=0$ The dissipation due to viscosity is represented by the magenta dashed line in Fig. 7. We can see that the spectra in this case exhibit a smooth profile, devoid of any pronounced peaks at specific frequencies, unlike the four-layer model. Additionally, it is worth noting that the dissipation in this case is significantly weaker, ranging from two to four orders of magnitude lower.

We also focused on the influence of Ekman number variations on dissipation spectra. Our results reveal that varying the Ekman number has a significant impact on the energy dissipation. Indeed, we find an increase in the number of peaks in the dissipation spectra as the Ekman number decreases. More specifically, the decrease in the Ekman number leads to lower viscosity, which results in higher and narrower resonance peaks associated with gravito-inertial modes, making the spectrum more complex, whereas all peaks are smoothed when a higher viscosity is used. This result is consistent with the predictions of Auclair Desrotour et al. (2015), which studied the dissipation of gravito-inertial waves by viscosity in a Cartesian box and thermal diffusion in a stably stratified medium. We also explored the impact of the Schmidt (Sc = {0.5,6}) and Prandtl (Pr = {0.5,2}) numbers in Figs. 8 and 9, respectively. We find that decreasing the Schmidt (Prandtl) number increases the molecular (thermal) dissipation. However, the total dissipation is not modified, since the viscous dissipation is dominant in this parameter regime.

5.3 Quality factor and Love number

To establish a connection between our numerical computations and observations, it was necessary to calculate the imaginary component of the Love number from the total dissipation. This calculation enabled us to conduct a comprehensive analysis by comparing our numerical models with actual observations and performing quantitative and qualitative comparisons.

The imaginary part of the Love number ℑ [k_ℓm] plays a crucial role in characterising the response of a celestial body to tidal forces, capturing the phase difference between the applied tidal forcing and the resulting response. It represents the transfer of energy and angular momentum within the system. By establishing a relation between the overall dissipation and the imaginary part of the Love number, we gain valuable insights into the evolution of the system. This relation can be expressed as in Ogilvie (2013): $$\mathfrak{J}\left[k_{\ell m}(\tilde{\omega })\right]=\frac{GM}{4R^3{\text{Ω}}^2}\frac{8\pi }{(2\ell +1)\tilde{\omega }}\mathcal{D};$$(91)

then, we can define the modified tidal quality factor as $$Q_{\ell m}^{\prime }(\tilde{\omega })=\mathrm{sign}(\tilde{\omega })\frac{3}{2\Im \left[k_{\ell m}(\tilde{\omega })\right]},$$(92)

which has the advantage of combining the tidal quality factor Q with the real part of the Love number ℜ[k_ℓm]: $$\tilde{J}\left[k_{\ell m}(\tilde{\omega })\right]=\mathrm{sign}(\tilde{\omega })\frac{\Re \left[k_{\ell m}(\tilde{\omega })\right]}{Q_{\ell m}(\tilde{\omega })}.$$(93)

Lainey et al. (2009) fitted a dynamical model including parameterised tidal dissipation to astrometric observations of the Galilean satellites from 1891 to 2007. They found that ℑ [k₂₂] = 1.1 × 10⁻⁵ $\left({Q^{\prime }}_{22}=-1.4\times {10}^5\right)$, for the asynchronous tide (ℓ = m = n = 2) due to Io $(\tilde{\omega }=-0.76)$.

Using Eqs. (91) and (92), we computed the imaginary part of the Love number and the modified tidal quality factor, and we represent them as a function of the normalised forcing frequency in Fig. 10 and the middle panel of Fig. 11 for m = 2, E = 10⁻⁷, and Pr = Sc = 1. We find a significant discrepancy between computed values of the imaginary part of the Love number (the modified tidal quality factor) due to Io and the observed ones, differing by roughly two orders of magnitude (the computed imaginary part of the Love number is 8.8 × 10⁻⁴, whereas the observed one is 1.1 × 10⁻⁵). Consequently, our calculations tend to overestimate the amplitude of tidal dissipation. Conversely, when examining the purely convective model, we observe an underestimation of tidal dissipation by approximately one and a half orders of magnitude.

We investigated the influence of varying Ekman number on the imaginary part of the Love number. The outcomes are illustrated in Fig. 11 for E = {10⁻⁶, 10⁻⁷, 10⁻⁸}. We find that decreasing the Ekman number impacts the imaginary part of the Love number (the dissipation) by increasing the number of peaks. However, for the frequency associated with the forcing imposed by Io, the impact is small because we are not on a resonance (a peak), and the imaginary part of the Love number does not vary significantly. We find the same result for the Schmidt and Prandtl numbers. Their impact on the total dissipation is very weak; therefore, they do not influence the imaginary part of the Love number.

Eventually, we find that stable stratification plays a crucial role in explaining the high dissipation. This conclusion was also highlighted by André et al. (2019), who investigated tidal dissipation in a rotating semi-convective region with a Cartesian box model. In addition, Lin (2023) and Dewberry (2023) have also studied tidal responses in some simplified scenarios conceivable for Jupiter’s interior with stably stratified layers, taking only the viscous diffusion into account. Our results confirm their results taking into account the three possible diffusion mechanisms, which are dominated by the viscous one, and more realistic internal structure models for Jupiter.

Fig. 13 Same as middle panel of Fig. 11, but for η = 0.1.

Fig. 14 Same as Fig. 1, but for a structure model with a smaller internal stably stratified layer.

5.4 Impact of the external stably stratified layer: Four-zone versus two-zone models

In Fig. 12, we observe the distribution of kinetic energy in a two-zone interior model. We can see that the internal part of the model follows a gravito-inertial pattern, while the external zone exhibits a single inertial mode. This discrepancy, in comparison with the left panel of Fig. 5, arises due to the absence of the external stably stratified layers, which theoretically facilitate wave reflection and the formation of two distinct inertial modes. Nevertheless, we can see that the attractor’s presence remains consistent, independently of the presence of the narrow intermediate stably stratified region. Furthermore, we find that the impact of this zone on the dissipation is very weak.

5.5 Impact of the size of solid and diluted cores

In order to study the impact of the size of the diluted and solid cores on the total dissipation, we first carried out a test with the five-layer model (Sect. 3.1), but with a bigger solid core (smaller diluted core) of 10% instead of 1.4%. As shown in Fig. 13, we find that the magnitude of imaginary part has slightly decreased (the imaginary part of the Love number due to Io is 5.4 × 10⁻⁴) and that the position of the peaks is only slightly modified. We performed another test with the two-layer model, but this time with a solid core of size 15% instead of 1.4%. We find that the position of the peaks changes and the dissipation due to Io increases by less than half an order of magnitude (with this model, the imaginary part of the love number due to Io is 2.5 × 10⁻³). Afterwards, we used another structure model that also satisfies Juno constraints and uses the equation of state of Chabrier & Debras (2021) where we reduce the size of the internal stably stratified layer (diluted core). As we can see in Fig. 14, the size of the diluted is reduced, and it is now localised between 48% and 56% of the radius (the outer stably stratified layer is in the same position). We can also see that the stratification in this model is stronger (N² /(2Ω)² ≈ 10 instead of ≈ 4). With this model, we find, as we can see in Fig. 15, that the imaginary part of the love number due to Io is closer to the observed value (with this model the imaginary part of the Love number due to Io is 1.4 × 10⁻⁴), but the gap remains significant (approximately one order of magnitude).

Fig. 15 Same as middle panel of Fig. 11, but for the structure model represented in Fig. 14.

6 Discussion and conclusions

We developed a numerical method that enables the calculation of the forced dynamic tidal response of an incompressible, non-magnetised, uniformly rotating fluid body. The Coriolis force is fully accounted for in our calculations. However, we do not consider centrifugal distortion, which allows us to solve the problem using spherical geometry. We take into consideration various types of dissipations such as fluid viscosity, thermal dissipation, and molecular diffusivity. By incorporating these dissipation mechanisms, we compute, using 2D numerical simulations, the total dissipation and determine the imaginary part of the tidal Love numbers for a given complex planetary interior model. In this study, we examine the dynamical tides in the latest Jupiter interior model (Sect. 3) and specifically investigate the quadrupolar tidal components (ℓ = m = 2). Our focus is on the frequency range that corresponds to the tidal frequencies associated with Jupiter’s Galilean moons. We considered a multi-layer model with alternating convective and stably stratified regions, which enables a more comprehensive and realistic representation of the physical processes occurring within giant gaseous planets’ interiors, in particular the dissipation of dynamical tides. We find that the presence of stably stratified regions plays a significant role in explaining the strong dissipation observed on Jupiter when compared to the case of a sole convective envelope. In this framework, we find that the dissipation depends on the chosen internal structure, in particular the size of the diluted core. In fact, with a large diluted core (around 68% of the radius), we find a discrepancy between the calculated and observed dissipation of two orders of magnitude due to Io, whereas with a smaller stably stratified inner layer (around 8% of the radius) the discrepancy becomes smaller (one order of magnitude). This may provide constraints on the size of the diluted core in the future. Our analysis also reveals that, in the chosen parameter regime in which the kinematic viscosity and thermal and molecular diffusivities are uniform and equal (the realistic variation in transport coefficients vary by several orders of magnitude and their ratios are potentially different from 1, depending on the considered region), the dominant mechanism contributing to dissipation is viscosity, surpassing both thermal and chemical dissipations in magnitude. Furthermore, it is important to note that our model is not limited to Jupiter, but it can also be applied to other giant planets such as Saturn and exoplanets.

There are several caveats that should be carefully considered in future studies in order to ensure accurate quantitative comparisons with high-precision observations, it is crucial to incorporate the relevant missing physical processes in a self-consistent manner. First, neglecting the influence of centrifugal effects may limit the accuracy of our solutions. Particularly for high-degree tidal components, the impact of centrifugal forces becomes increasingly significant (Dewberry 2023). Second, while adopting the Boussinesq approximation to investigate dynamical tides simplifies the system of equations to solve, it is important to acknowledge its limitations. These limitations are particularly significant when the Lamb frequency, which characterises the acoustic modes, approaches a comparable magnitude to the excited mode frequencies near the surface. Clearly, an important follow-up of this work would be to go from the Boussinesq approximation to the anelastic approximation and take into account density stratification. The outcomes of using the more realistic anelastic approximation are not expected to completely deviate from those obtained with the Boussinesq approximation; in fact, both approximations yield the same attractors of characteristics. This comes about because in the anelastic approximation the velocity υ is replaced by the specific linear momentum ρυ in the system of equations. This means that the momentum vector satisfies the same set of equations as the velocity vector does in the Boussinesq case (Dintrans & Rieutord 2000). In this respect, using simple polytropic models, the work of Ogilvie (2013) gives a first exploration of the effects of density variations of the background on tidal dissipation. The Boussinesq approximation may overestimate the tidal dissipation that could explain why the computed dissipation in our work is too large when compared to the observations. This will be carefully evaluated in a follow-up work where we shall use the anelastic approximation.

Finally, differential rotation can play a crucial role in the dynamics of the outer regions of gas giant planets. In addition, the ionised inner region, characterised by the presence of a magnetised gas, can exhibit significant effects due to ohmic dissipation and induced magnetic torques. We know that the presence of differential rotation in a convective zone is strongly dependent on electrical conductivity (e.g. Guillot et al. 2018; Galanti et al. 2019). Therefore, an interesting follow-up of this work would be to undertake a study to understand the profound impact of both differential rotation (Mathis 2009; Baruteau & Rieutord 2013; Mirouh et al. 2016; Guenel et al. 2016a,b; Dewberry et al. 2021) and magnetic fields (Rogers & MacGregor 2010; Mathis & de Brye 2011; Barker & Lithwick 2014; Wei 2016, 2018; Lin & Ogilvie 2018) on wave propagation and dissipation in gas giant planets.

Acknowledgements

We thank the referee for her/his positive and constructive report, which has allowed us to improve the quality of our article. H. Dhouib and S. Mathis acknowledge support from the CNES PLATO grant at CEA/DAp and from PNP (CNRS/INSU). A. Astoul acknowledges support from the Science and Technology Facilities Council (STFC) grant ST/S000275/1, as well as the Leverhulme Trust for early career grant. We are also very grateful to Lorenzo Valdettaro for his kind help and support with the LSB code. Some preliminary calculations have been performed thanks to HPC resources from CALMIP super-computing centre (Grants P16024 and 2022-P0107). S. Mathis and M. Rieutord acknowledge support from the European Research Council through HORIZON ERC SyG Grant 4D-STAR 101071505.

Appendix A Dimensionless number expressions in Jupiter

The expressions for the different diffusivities and associated dimensionless numbers are derived in Stevenson & Salpeter (1977) for metallic and molecular phases. The region with R < 0.9 exhibits a metallic phase, whereas the region with R > 0.9 is characterised by a molecular phase. We added a smooth transition between both zones. We reiterate the expression for these numbers here: $$\mathrm{Pr}=\{\begin{array}{ll}\frac{4}{0.3}{T}^{-1/2}{\rho }^{-1/3}\hfill & \text{ if }R<0.9\hfill \\ 1\hfill & \text{ if not }\hfill \end{array},$$(A.1) $$\text{Sc}=\{\begin{array}{ll}\frac{4}{3}{10}^6{T}^{-2}{\rho }^{2/3}\hfill & \text{ if }R<0.9\hfill \\ \frac{1}{4}{10}^{6.5}{T}^{-2}{\rho }^{5/6}\hfill & \text{ if not }\hfill \end{array}\text{, }$$(A.2) $$\text{E}=\{\begin{array}{ll}2\times {10}^{-5}\frac{T^{-1/2}}{\text{Ω}{R}^2}\hfill & \text{ if }R<0.9\hfill \\ 5\times {10}^{-5}\frac{T^{-1/2}}{\text{Ω}{R}^2}\hfill & \text{ if not }\hfill \end{array}.$$(A.3)

All quantities must be expressed in SI units.

Appendix B Expression of some operators in the spherical harmonics basis

The spherical vector harmonics form a complete family, and the orthogonality relations ensure that any sufficiently regular vector field u can be uniquely expanded over the spherical vector harmonics (Rieutord 1987): $$u=\left[u_{\ell }^m,v_{\ell }^m,w_{\ell }^m\right].$$(B.1)

Appendix B.1 General case

The divergence of the vector field u $$\mathbf{\nabla }\cdot u=\frac{1}{r^2}{\partial }_r\left(r^2{u}_{\ell }^m\right)-\frac{\ell (\ell +1)}{r}{v}_{\ell }^m.$$(B.2)

The curl of the vector field u is expressed as follows: $$\mathbf{\nabla }\times u=\left[\ell (\ell +1)\frac{w_{\ell }^m}{r},\frac{1}{r}{\partial }_r\left(r{w}_{\ell }^m\right),\frac{u_{\ell }^m}{r}-\frac{1}{r}{\partial }_r\left(r{v}_{\ell }^m\right)\right].$$(B.3)

Appendix B.2 Particular case : a vector field with zero divergence

If ∇⋅u = 0., the curl of the vector field u can be rewritten as $$\mathbf{\nabla }\times u=\left[\ell (\ell +1)\frac{w_{\ell }^m}{r},\frac{1}{r}{\partial }_r\left(r{w}_{\ell }^m\right),-\frac{{\text{Δ}}_{\ell }\left(r{u}_{\ell }^m\right)}{\ell (\ell +1)}\right],$$(B.4)

with ${\text{Δ}}_{\ell }\square ={\partial }_{r^2}^2\square +\frac{2}{r}{\partial }_r\square -\frac{\ell (\ell +1)}{r^2}\square$.

The Laplacian of the vector field u is expressed as follows: $${\mathbf{\nabla }}^{2}u=\left[\frac{1}{r}{\text{Δ}}_{\ell }\left(r{u}_{\ell }^m\right),\frac{1}{r}{\partial }_r\left(\frac{D_{\ell }{u}_{\ell }^m}{\ell (\ell +1)}\right),{\text{Δ}}_{\ell }{w}_{\ell }^m\right],$$(B.5)

with $D_{\ell }\square ={\partial }_r^2\left(r^2\square \right)-\ell (\ell +1)\square$.

The vector product between the unit vector e_z and the vector field u is given by $$\begin{array}{l}{\mathit{e}}_{\mathbf{z}}\times \mathit{u}=[(\ell -1){\alpha }_{\ell -1}^{\ell }{w}_{\ell -1}^m-(\ell +2){\alpha }_{\ell +1}^{\ell }{w}_{\ell +1}^m-im{v}_{\ell }^m,\\ \frac{\ell -1}{\ell }{\alpha }_{\ell -1}^{\ell }{w}_{\ell -1}^m+\frac{\ell +2}{\ell +1}{\alpha }_{\ell +1}^{\ell }{w}_{\ell +1}^m-\frac{im}{\ell (\ell +1)}\left(u_{\ell }^m+v_{\ell }^m\right),\\ -\frac{{\alpha }_{\ell -1}^{\ell }}{{\ell }^2}{r}^{\ell -1}{\partial }_r\left(\frac{u_{\ell -1}^m}{r^{\ell -2}}\right)-\frac{{\alpha }_{\ell +1}^{\ell }}{{(\ell +1)}^2}{r}^{-\ell -2}{\partial }_r\left(r^{\ell +3}{u}_{\ell +1}^m\right)-\\ \frac{im}{\ell (\ell +1)}{w}_{\ell }^m]\end{array}$$(B.6)

References

All Figures

	Fig. 1 Profiles of compositional $\left(N_{\mu }^2\right)$, thermal $\left(N_t^2\right)$, and total (N2) Brunt-Väisälä frequencies squared normalised by the inertial frequencies squared (4Ω2) as a function of the normalised radius r/R.
In the text

	Fig. 2 Schematic of model of Jupiter’s interior used in this study.
In the text

	Fig. 3 Prandtl, Schmidt, and Ekman number profiles as a function of the normalised radius based on Stevenson & Salpeter (1977).
In the text

	Fig. 4 Non-wave-like (equilibrium) gravitational potential (left) and corresponding radial and horizontal displacements (right) as a function of the normalised radius for A = 1 and ℓ = 2.
In the text

Fig. 5 Meridional cut (z = r cos θ and s = r sin θ) of dimensionless kinetic energy (left), potential energy associated with chemical (middle), and potential energy associated thermal stratification (right) of the forced mode $\tilde{\omega }=-0.76$ for m = 2, E = 10−7, Pr = Sc = 1 with a spatial resolution of (Nr, Nℓ) = (300, 301). The trajectories of characteristics are represented with white curves, while the surfaces at which they undergo reflections are shown with brown curves.

In the text