© The Authors 2023

1 Introduction

Tidal interactions between Jupiter and the Galilean satellites play an important role in the orbital evolution of the system and the internal dynamics of the moons (Lainey et al. 2009). The highly active volcanic eruptions on Io are believed to be due to strong tides raised by Jupiter (Peale et al. 1979). Meanwhile, tides are also raised in Jupiter by its moons, which are probably dominated by Io (Gavrilov & Zharkov 1977). The tidal response of a gaseous body such as Jupiter is conventionally treated as a hydrostatic deformation, which acquires a small phase lag with respect to the tidal forcing due to dissipative processes. This is known as the equilibrium tide. However, the equilibrium tide alone does not suffice to account for 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).

In fact, the equilibrium tide does not satisfy the momentum equation of tidal flows and thus corrections have to be made to fully account for the tidal response of Jupiter. The corrections to the equilibrium tide are collectively referred to as the dynamical tide, which usually involves wave-like motions in the planet and depends on the tidal frequency as well as the interior structure (Ogilvie 2014). The dynamical tide may provide extra channels of tidal dissipation and produce additional gravitational perturbations in addition to the hydrostatic deformation. The Juno spacecraft has obtained highly accurate tidal Love numbers, k_lm (Durante et al. 2020), which quantitatively characterize the tidal response of Jupiter to a tidal forcing component represented in spherical harmonics of degree l and order m. The observed tidal Love numbers by Juno exhibit non-negligible discrepancies with respect to the theoretically calculated hydrostatic values (Wahl et al. 2020), suggesting that the dynamical tide has to be considered to explain the observed tidal response. Specifically, Juno observations found ∆k₂₂ ≈ −4% for the dominant tidal component l = 2 and m = 2 and ∆k₄₂ ≈ −11% for the high-degree tidal component l = 4 and m = 2, where $\text{Δ}{k}_{lm}=\left(k_{lm}-k_{lm}^{\left(hs\right)}\right)/k_{lm}^{\left(hs\right)}$ represents the fractional correction to the hydrostatic value $k_{lm}^{\left(hs\right)}$ (Wahl et al. 2020; Idini & Stevenson 2021, 2022a).

As the dynamical tides are sensitive to the tidal frequency and the interior structure, the detected gravitational signatures of dynamical tides may provide important constraints on Jupiter’s interior (Lai 2021; Idini & Stevenson 2021, 2022b; Dewberry & Lai 2022). Recent studies (Idini & Stevenson 2021; Lai 2021) have revealed that the discrepancy in k₂₂ can be mainly attributed to the Coriolis effect on the fundamental modes (f-modes). More recently, Idini & Stevenson (2022b) proposed that the resonant locking with a gravity mode in an extended dilute core can explain ∆k₄₂ ≈ −11%. This provides an independent constraint on the existence of a dilute core in Jupiter, which has also been suggested by the Juno measurements of gravitational moments of Jupiter (Wahl et al. 2017; Militzer et al. 2022). However, the tidal constraint on the existence of a dilute core is uncertain. The calculation of the tidal response in Idini & Stevenson (2022b) inadequately treated the rotational (Coriolis) effect, which plays an important role in Jupiter’s tidal responses because the tidal frequencies of Galilean satellites are comparable to the spin frequency of Jupiter. Including the Coriolis force introduces inertial waves in the neutrally buoyant regions (Ogilvie & Lin 2004; Wu 2005a) as well as mixed gravity waves and inertial waves (gravito-inertial waves) in the stably stratified region (Dintrans et al. 1999; Xu & Lai 2017). The mechanism proposed by Idini & Stevenson (2022b) is also struggling to reconcile both the real part (relevant to the gravitational perturbation) and imaginary part (relevant to the tidal dissipation) of the tidal Love numbers.

For this study, we developed a method to directly calculate the tidal response of a fully compressible, self-gravitating, rotating, and viscous fluid body. The Coriolis force is fully taken into account, but the centrifugal force is neglected, which allows us to numerically solve the problem in spherical geometry using a pseudo-spectral method based on spherical harmonic expansions (Ogilvie & Lin 2004; Lin & Ogilvie 2017). As we directly solve the tidally forced problem with explicit viscosity, we could simultaneously obtain the real and imaginary parts of the tidal Love number for a given planetary interior model. Our approach is different from recent studies on dynamical tides of Jupiter (Lai 2021; Idini & Stevenson 2022b; Dewberry & Lai 2022). They obtained the eigen modes of the inviscid fluid body first and then calculated the tidal Love number (only the real part) through projecting the tidal force onto each eigen modes. We consider three nominal interior models of Jupiter to investigate the dependence of the tidal response on the tidal frequency and the interior structure. We focus on the effect of a compact rigid core, an extended dilute core, and a thin stably stratified layer in the outer region on tidal responses. All of the simplified models can explain the observed ∆k₂₂ ≈ −4% as previous studies have shown. However, these simplified models cannot account for the observed ∆k₄₂ ≈ −11%. Resonances with gravito-inertial modes in an extended dilute core near the tidal frequency of Io can produce non-negligible dynamical correction to k₄₂, but it is insufficient to explain the Juno observation based on our simplified model.

2 Tidal model

We consider linear tidal responses of a rotating gaseous planet to a tidal potential component of ${\Psi }_l^m=𝒜{\left(r/R\right)}^l{Y}_l^m\left(\theta ,\varphi \right){\text{e}}^{-\text{i}\omega t}$, where 𝒜 is the tidal amplitude, R is the radius of the planet, $Y_l^m\left(\theta ,\varphi \right)$ represents spherical harmonics, and ω is the tidal frequency. The resulting tides of the planet produce an external gravitational potential perturbation $\text{Φ}\prime =ℬ{\left(R/r\right)}^{l+1}{Y}_l^m\left(\theta ,\varphi \right){\text{e}}^{-\text{i}\omega t}$ (and probably other spherical harmonic components). The ratio $K_l^m\left(\omega \right)=ℬ/𝒜$ defines the tidal Love number, which depends on the tidal frequency. The tidal Love number $K_l^m$ is a complex number because there exists a phase lag between the forcing and the gravitational perturbations due to dissipative processes (Ogilvie 2014). While the real part $k_{lm}=\mathrm{Re}\left[K_l^m\right]$ measures the in-phase gravitational perturbations with the tidal forcing, the imaginary part $\mathrm{Im}\left[K_l^m\right]$ quantifies the out-of-phase tidal response and is related to the dissipation rate. The ratio between the real and imaginary parts is related to the tidal quality factor $$Q=\mathrm{sgn}\left(\omega \right)\frac{k_{lm}}{\mathrm{Im}\left[K_l^m\right]},$$(1)

where sgn(ω) = ±1 is the sign function. Because the phase lag is generally very small, that is Q ≫ 1, the magnitude of the imaginary part is typically much smaller than the real part. For this study, we develope an approach to directly and simultaneously calculate the real and imaginary parts of the tidal Love number for a given planetary model.

2.1 Linearized equations

For a compressible, self-gravitating, and rotating fluid body which may contain a rigid core of radius R_i, linear perturbations to a tidal potential Ψ ∝ e^−iωt in the rotating frame are described by the following equations (e.g., Ogilvie & Lin 2004): $$-i\omega u^{\prime }=-2\Omega \times u^{\prime }-\frac{1}{{\rho }_0}\nabla P^{\prime }+\frac{{\rho }^{\prime }}{{\rho }_0^2}\nabla P_0-\nabla {\text{Φ}}^{\prime }-\nabla \text{ψ+}{f}_v,$$(2) $$-i\omega {\rho }^{\prime }+\nabla \cdot \left({\rho }_0{u}^{\prime }\right)=0,$$(3) $$-i\omega \left(\frac{P^{\prime }}{\text{Γ}{P}_0}-\frac{{\rho }^{\prime }}{{\rho }_0}\right)+u^{\prime }\cdot \left(\frac{1}{\text{Γ}}\nabla \ln P_0-\nabla \ln {\rho }_0\right)=0$$(4) $${\nabla }^2\text{Φ}\prime =4\pi G{\rho }^{\prime },$$(5)

where u is the velocity, Ω the rotation rate, ρ the density, P the pressure, Γ the adiabatic index, and G the gravitational constant. In the above equations, the subscript ₀ denotes physical quantities in the hydrostatic state (without tidal potential) and the notations with the prime represent Eulerian perturbations induced by the tidal forcing. In the momentum Eq. (2), we have explicitly included a viscous force f_ν defined as $$f_v=\frac{1}{{\rho }_0}\nabla \cdot \left(2\mu S\right),$$(6)

where µ is the dynamic shear viscosity (we neglected the bulk viscosity) and S is the strain-rate tensor: $$S=\frac{1}{2}\left[\nabla u^{\prime }+{\left(\nabla u^{\prime }\right)}^T\right]-\frac{1}{3}\left(\nabla \cdot u^{\prime }\right)I.$$(7)

We have included the viscous force in the momentum equation, but we neglected the viscous heating in the energy equation, that is the density and pressure perturbations were treated as adiabatic.

For this study, we have fully taken the Coriolis force into account, but neglected the centrifugal distortion for numerical convenience. The centrifugal effect can be measured by ϵ = Ω/ω_dyn, that is the ratio between the spin frequency Ω and the dynamical frequency ω_dyn = (GM/R³)^1/2. This ratio is not particularly small for Jupiter (ϵ = 0.288). Indeed, the centrifugal distortion of Jupiter has non-negligible contributions to the total Love number k_lm, especially for the high-degree Love number k₄₂ because the tidal response at l = m = 2 can produce a gravitational perturbation at l = 4 and m = 2 in an oblate figure (Idini & Stevenson 2022a). For the hydrostatic $k_{42}^{\left(\text{hs}\right)}$ of Jupiter due to Io, 93% of the total value is actually contributed by the centrifugal coupling with k₂₂ and only the remaining 7% is produced by the tidal forcing at l = 4 and m = 2 (Wahl et al. 2020; Idini & Stevenson 2022a). In this paper we do not aim to directly fit the k_lm observed by Juno, but rather focus on the fractional corrections ∆k_lm by the dynamical tides. In terms of the fractional correction ∆k_lm, the centrifugal contribution to ∆k₂₂ can be neglected in leading order (Lai 2021). However, the centrifugal contribution to ∆k₄₂ cannot be neglected even in leading order because $k_{42}^{\left(\text{hs}\right)}$ is mostly contributed to by the centrifugal coupling with k₂₂. This complicates the comparison between the calculated ∆k₄₂ in a spherical figure and the observation. Nevertheless, the calculated ∆k₄₂ in a spherical figure can be multiplied by the factor 0.07 to account for Jupiter’s centrifugal coupling effect for qualitative comparisons with the observation (Idini & Stevenson 2022b). Such a comparison would assume that the tidally excited internal modes are not significantly modified by the centrifugal deformation.

By neglecting the centrifugal deformation, the unperturbed basic state is spherically symmetric, that is to say it depends on the radius r only. Given the density ρ₀(r) and pressure P₀(r) profiles of the unperturbed state, the radial gravitational acceleration (inward) g(r) and the Brunt-Väisälä frequency N(r) are then determined by $$g=\frac{{\text{dΦ}}_0}{\text{d}r}=-\frac{1}{{\rho }_0}\frac{\text{d}{P}_0}{\text{d}r},$$(8) $$N^2=g\left(\frac{1}{\text{Γ}}\frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}{P}_0}{\text{d}r}-\frac{\text{d\hspace{0.17em}}\ln {\rho }_0}{\text{d}r}\right).$$(9)

2.2 Numerical method

In order to obtain the complex Love numbers, we numerically solved Eqs. (2)–(5) using a pseudo-spectral method for the prescribed basic states, which are subject to the relevant boundary conditions. The numerical scheme is based on the method used in previous studies (Ogilvie & Lin 2004; Lin & Ogilvie 2017), but we extended the method to solve the full set of linearized Eqs. (2)–(5) without making a low-frequency approximation (Ogilvie 2013). By introducing h′ = P′/ρ₀ and eliminating the density perturbation ρ′, Eqs. (2)–(5) can be reduced to the following equations: $$\begin{array}{ll}-i\omega {\rho }_0{u}^{\prime }=\hfill & -2{\rho }_0\Omega \times u^{\prime }-\nabla \left({\rho }_0{h}^{\prime }\right)+g{\nabla }^2{\phi }^{\prime }/\left(4\pi G\right)\hfill \\ \hfill & -{\rho }_0\nabla {\phi }^{\prime }-\nabla \text{ψ+}\nabla \cdot \left(2\mu S\right),\hfill \end{array}$$(10) $$-i\omega h^{\prime }=-c_s^2\left(N^2{u^{\prime }}_r/g+\nabla \cdot \left({\rho }_0{u}^{\prime }\right)/{\rho }_0\right),$$(11) $$-i\omega {\nabla }^2{\phi }^{\prime }=-4\pi G\nabla \cdot \left({\rho }_0{u}^{\prime }\right),$$(12)

where ${u^{\prime }}_r$ is the radial velocity perturbation and $c_s^2=\text{Γ}{P}_0/{\rho }_0$ is the square of the adiabatic sound speed.

We imposed boundary conditions including the regularity of the gravitational perturbations, zero radial velocity on the rigid inner boundary, and vanishing Lagrange pressure perturbation at the surface, that is $\delta P=P^{\prime }+{u^{\prime }}_r/\left(-i\omega \right)\nabla P_0=0$. In terms of h′ and ${u^{\prime }}_r$, the last boundary condition can be written as (Dewberry et al. 2021) $$\left(-i\omega h^{\prime }-g{u^{\prime }}_r\right)|{}_{r=R}=0.$$(13)

As the viscous force was included, additional boundary conditions were required to complete the boundary value problem. We used the so-called stress-free conditions, in other words the tangential stresses vanish at both boundaries.

For a given tidal potential ${\text{Ψ}}_l^m=𝒜{\left(r/R\right)}^l{Y}_l^m\left(\theta ,\varphi \right){\text{e}}^{-\text{i}\omega t}$, the tidal perturbations (including both equilibrium and dynamical tides) u′, h′ and Φ′ can be expanded as $$u^{\prime }={\displaystyle \sum _{n=m}^L{u}_n^m}\left(r\right)R_n^m+{\displaystyle \sum _{n=m}^L{\upsilon }_n^m}\left(r\right)S_n^m+{\displaystyle \sum _{n=m}^L{w}_n^m}\left(r\right)T_n^m,$$(14) $$h^{\prime }={\displaystyle \sum _{n=m}^L{h}_n^m}\left(r\right)Y_n^m\left(\theta ,\varphi \right),$$(15) $$\text{Φ}\prime ={\displaystyle \sum _{n=m}^L{\text{Φ}}_n^m}\left(r\right)Y_n^m\left(\theta ,\varphi \right),$$(16)

where $R_n^m,S_n^m,T_n^m$ are vector spherical harmonics $$R_n^m=Y_n^m\left(\theta ,\varphi \right)\widehat{r},S_n^m=r\nabla Y_n^m\left(\theta ,\varphi \right),T_n^m=r\nabla \times R_n^m.$$(17)

As the basic state is axisymmetric, the perturbations involve spherical harmonics with the same order m as the tidal potential ${\text{Ψ}}_l^m$, but the Coriolis force would couple all spherical harmonics with degree n ≥ m. For numerical calculations, we had to make a truncation at certain degree L. Substituting expansions of Eqs. (14)–(16) into Eqs. (10)–(12) and projecting onto spherical harmonics, we ended up with a set of ordinary differential equations (ODEs) involving $u_n^m\left(r\right),{\upsilon }_n^m\left(r\right),w_n^m\left(r\right),h_n^m\left(r\right)$, and ${\text{Φ}}_n^m\left(r\right)$. For the radial dependence, we used Chebyshev collocation on N_r Gauss-Lobatto nodes (Rieutord et al. 2001). The boundary conditions were applied through replacing the ODEs with the corresponding boundary conditions on the boundary nodes. The regularity of gravitational perturbations requires $$r\frac{{\text{dΦ}}_n^m}{\text{d}r}+\left(n+1\right){\text{Φ}}_n^m=0\text{\hspace{1em}at\hspace{0.17em}}r=R,$$(18) $$r\frac{{\text{dΦ}}_n^m}{\text{d}r}-n{\text{Φ}}_n^m=0\text{\hspace{1em}at\hspace{0.17em}}r=R_i.$$(19)

The vanishing Lagrangian pressure perturbation at the surface and zero radial velocity at the rigid inner boundary give $$-i\omega h_n^m=g{u}_n^m\text{\hspace{1em}at\hspace{0.17em}}r=R,$$(20) $$u_n^m=0\text{\hspace{1em}at\hspace{0.17em}}r=R_i.$$(21)

The stress-free boundary condition is given as (Ogilvie 2009) $$u_n^m+r\frac{\text{d}{\upsilon }_n^m}{\text{d}r}-{\upsilon }_n^m=0,r\frac{\text{d}{w}_n^m}{\text{d}r}-w_n^m=0$$(22)

at both boundaries.

Using the numerical discretization described above, the boundary value problem becomes a linear system involving a large complex block-tridiagonal matrix. The solution of the linear system was obtained using the standard direct solver. We used typical truncations of L = 200 and N_r = 100 for this study.

Once the solution of the linear system is obtained numerically, the complex tidal Love number is readily given by $$K_l^m={\text{Φ}}_l^m\left(r=R\right)$$(23)

for the tidal potential component ${\text{Ψ}}_l^m=𝒜{\left(r/R\right)}^l{Y}_l^m\left(\theta ,\varphi \right){\text{e}}^{-\text{i}\omega t}$ (we simply set 𝒜 = 1 for the linear tidal response). We note that the solution includes both the equilibrium and dynamical tides. For the real part of Love numbers, of particular interest is the fractional correction of dynamical tides $$\text{Δ}{k}_{lm}=\left(k_{lm}-k_{lm}^{\left(hs\right)}\right)/k_{lm}^{\left(hs\right)},$$(24)

where $k_{lm}^{\left(hs\right)}$ is the hydrostatic value and it was calculated by setting ω = 0. As our calculations neglect the centrifugal effect which significantly influences the high-degree Love number k₄₂, the calculated value of ∆k₄₂ should be multiplied by the factor 0.07 when compared with the observation as we have discussed in Sect. 2.1.

We can also calculate the tidal dissipation rate D_ν from the velocity perturbations $$D_v={\displaystyle {\int }_{V}2\mu S^2}\text{d}V,$$(25)

where the integral was taken over the fluid domain. The dissipation rate is related to the imaginary part of the tidal Love number (Ogilvie 2014) $$D_v=\frac{\left(2l+1\right)R{\mathcal{A}}^2}{8\pi G}\omega \mathrm{Im}\left[K_l^m\right],$$(26)

which can be used as an independent validation of the numerical code. The above relation is satisfied to a high degree of accuracy for all of numerical calculations presented in this paper.

Fig. 1 Three nominal models of Jupiter’s interior used in this study. The top panel shows the schematic models and the bottom panel shows the density (normalized by the density at the center) and the Brunt-Väisälä frequency (normalized by the dynamical frequency) as a function of the radius. The blue shadow in the bottom panel indicates solid regions. (a) A compact rigid core model; (b) an extended dilute core model; and (c) a compact rigid core and an outer stable layer model.

2.3 Interior models

In order to solve Eqs. (10)–(12), we need to prescribe basic state profiles ρ₀(r), g(r), and N²(r) to model Jupiter’s interior. Our understanding of Jupiter’s interior has been significantly improved by Juno observations (Stevenson 2020), yet some degrees of uncertainty remain. In this study, we do not aim to build a realistic model of Jupiter’s interior, but focus on the fractional contributions of dynamical tides to the tidal Love number for different possible scenarios of Jupiter’s interior. We consider three nominal interior models (Fig. 1) based on a polytrope of index 1, which is a good leading order approximation for Jupiter (Stevenson 2020).

For all of models used in this study, the unperturbed density and gravity follow a hydrostatic polytrope of index 1, $${\rho }_0=\frac{\pi M}{4R^3}\frac{\sin kr}{kr},$$(27) $$g=\frac{GM}{r^2}\left[\sin \left(kr\right)-kr\cos \left(kr\right)\right],$$(28)

where k = π/R. The first model consists of a small rigid core of radius 0.25R and an isentropic fluid envelope, that is Γ = 2 and N² = 0 in the fluid region (Fig. 1a).

The second model assumes an extended dilute core of radius 0.7R and an isentropic envelope (Fig. 1b). The dilute core is treated as a stably stratified fluid layer with the Brunt-Väisälä frequency given by $$\frac{N^2}{{\omega }_{\text{dyn}}^2}={\tilde{N}}_2\sin \left(\frac{\pi r}{R_c}\right),$$(29)

where Ñ₂ = 0.25 and R_c = 0.7 for this model. As we fixed the density and pressure profiles to that of a polytrope, the stratification was effectively realized by adjusting the adiabatic index (Γ > 2) in the dilute core (Lai 2021). This model is similar to the one used in Idini & Stevenson (2022b), but they adjusted the density profile to model the stable stratification in the dilute core while fixing the adiabatic index Γ = 2.

The third model is based on the model in Fig. 1a, but we further added a stably stratified layer between 0.8R and 0.9R (Fig. 1c), possibly resulting from H-He immiscibility (Debras & Chabrier 2019; Stevenson et al. 2022). The Brunt-Väisälä frequency in the top stable layer is prescribed as $$\frac{N^2}{{\omega }_{\text{dyn}}^2}={\tilde{N}}_2\frac{1}{\left[1+{\text{e}}^{-100\left(r-0.8\right)}\right]\left[1+{\text{e}}^{100\left(r-0.9\right)}\right]}.$$(30)

The degree of stratification of this layer remains uncertain, but it is estimated that typical values of $N^2/{\omega }_{\text{dyn}}^2$ would be roughly between 0.1 and 0.8 for Jupiter (Christensen et al. 2020; Gastine & Wicht 2021). Here we set a moderate value Ñ₂ = 0.5.

We note that an interior model with the coexistence of a dilute core and a top stable layer is also possible (Debras & Chabrier 2019). As this kind of model involves two different stably stratified layers, it would be difficult to characterize the role of the top stable layer on tides. We considered only the combination of a compact rigid core and a top stable layer for simplicity.

In all of these models, we set the total mass M, the radius R, and the spin rate Ω such that the ratio $ϵ=\text{Ω}/\sqrt{GM/R^3}=0.288$, corresponding the value of Jupiter. Our calculations also require the fluid viscosity, which is difficult to estimate in detail for giant planets. We simply assumed the dynamic viscosity µ is proportional to the background density ρ₀, so the kinematic viscosity ν = µ/ρ₀ is constant. The viscosity can be measured by the dimensionless number Ek = ν/(ΩR²), known as the Ekman number. We set Ek = 10⁻⁶ for most of the calculations (unless otherwise specified), roughly corresponding to the effective viscosity based on mixing-length theory (Guillot et al. 2004).

As we have mentioned that we do not aim to construct a realistic interior model for Jupiter in this study. These simplified models were designed to investigate the effects of a compact rigid core, an extended dilute core, and a top stable layer on the tidal responses of Jupiter. Nevertheless, the fractional corrections ∆k_lm and the tidal quality factor Q for these simplified models can be used to make some qualitative comparisons with the observations (Lai 2021; Idini & Stevenson 2021, 2022b).

Fig. 2 Complex Love number as a function of the tidal frequency for a full isentropic polytrope of index 1. Top panel: fractional correction ∆klm of the real part of the Love numbers. The fractional correction ∆k42 (orange curve in the top panel) was multiplied by 0.07. Bottom panel: minus imaginary part $-\mathrm{Im}\left[K_l^m\right]$, which is equivalent to klm/Ql. Vertical dashed lines indicate tidal frequencies of four Galilean moons of Jupiter (from right to left: Io, Europa, Ganymede, and Callisto). The horizontal dashed line in the bottom panel represents the astrometric observation of the frequency independent k2/Q2 from Lainey et al. (2009).

3 Results

In this paper, we focus on the dominant tidal component ${\text{Ψ}}_2^2$ and a high-degree tesseral component ${\text{Ψ}}_4^2$, for which non-negligible dynamical corrections have been detected as we have discussed in Sect. 1. Our calculations are limited to the frequency range of −2 ≤ ω/Ω ≤ −1, which is relevant to the tidal frequencies of the Galilean moons. The negative tidal frequency means that the tidal forcing is retrograde in the corotating frame with the planet based on our convention. For the real part of Love numbers, we show the fractional correction ∆k_lm. In order to make comparisons with the Juno observation, the calculated ∆k₄₂ was multiplied by 0.07 to compensate for the centrifugal effect which is neglected in our calculations. Because of the negative tidal frequency, the imaginary part of Love numbers is also negative in our calculations and is related to the tidal quality factor by $k_{lm}/Q_l=-\mathrm{Im}\left[K_l^m\right]$ according to Eq. (1).

3.1 Full polytrope model

Before presenting results for the interior models in Fig. 1, we first show the tidal response of a full isentropic polytrope, that is to say neutrally buoyant in the whole fluid sphere. This model serves as a reference for other models and has been used to investigate the dynamical tides of Jupiter in recent analytical studies (Idini & Stevenson 2021; Lai 2021). Figure 2 shows both the real and imaginary parts of the Love numbers as a function of the tidal frequency for the full polytrope model. We can see that ∆k₂₂ is negative in the frequency range we considered and it smoothly varies as the tidal frequency except at a burst around ω/Ω = −1.08, which corresponds to a resonance with an inertial mode. Away from resonances, our numerical results are consistent with recent theoretical calculations and produce ∆k₂₂ ≈ −4% at the tidal frequency of Io (Lai 2021; Idini & Stevenson 2021). These studies also revealed that the dynamical correction ∆k₂₂ can be attributed to the Coriolis effect on the f -modes. Apart from the f-modes, the rotating sphere of isentropic fluid also supports smooth inertial modes restored by the Coriolis force in the frequency range of 0 < |ω/Ω| < 2 (Greenspan 1968; Lockitch & Friedman 1999). The burst of ∆k₂₂ at ω/Ω = −1 08 is indeed due to the resonant excitation of the inertial mode as shown in Fig. 3a, but we noticed that the resonance occurs only in a very narrow frequency range.

However, this inertial mode has more significant contributions to ∆k₄₂. The angular structure of an inertial mode cannot be described by single spherical harmonics in general (Lockitch & Friedman 1999), but the density perturbations (and thus the gravitational perturbations) are dominated by the spherical harmonics $Y_4^2\left(\theta ,\varphi \right)$ for the resonant inertial mode at ω/Ω = −1.0836 as we can see from Fig. 3a. This suggests a likely strong coupling between the tidal potential component ${\text{Ψ}}_4^2$ and the inertial mode in Fig. 3a, that is to say large tidal overlap as described in Wu (2005b), leading to significant dynamical corrections to k₄₂. The dynamical correction can reach ∆k₄₂ ≈ −10% (after the centrifugal correction) near the resonance at ω/Ω = −1.0836. However, the tidal frequencies of the Galilean satellites are too far away from this resonance.

The curve of ∆k₄₂ also shows a spike around ω/Ω = −1.51, corresponding to a narrow resonance with a high-degree inertial mode (ρ′ is dominated by $Y_6^2\left(\theta ,\varphi \right)$ as shown in Fig. 3b). Interestingly, the tidal frequency of Io is close to this resonance, but the dynamical correction caused by this resonant mode is insufficient to account for the observed ∆k₄₂ ≈ −11%. The frequencies of inertial modes in Fig. 3 are slightly shifted compared to the results of Lockitch & Friedman (1999) for a polytrope of index 1 (see their table 6 and note the different conventions for the sign of frequencies) because they assumed ϵ → 0 whereas we set ϵ = 0.288.

The imaginary parts of the Love numbers in Fig. 2 show that resonances with inertial modes significantly enhance the tidal dissipation. The enhanced dissipation due to resonant inertial modes in a neutrally buoyant sphere has been demonstrated by Wu (2005b), but using different density profiles. When the tidal frequency is away from resonances, the dissipation rate for the full isentropic polytrope is too small to account for the observed tidal quality factor Q (Lainey et al. 2009).

Fig. 3 Density perturbations (left half) and radial velocity perturbations (right half) in the meridional plane to the tidal component ${\text{Ψ}}_4^2$ at two resonant frequencies in Fig. 2. Amplitudes were normalized by the maximum absolute values.

3.2 Compact rigid core model

We now consider tidal responses for the interior model with a compact rigid core. Basically, the inner region (r ≤ 0.25R) of a whole fluid polytrope becomes solid for this model. Figure 4 shows the frequency dependence of the Love numbers for the compact rigid core model. We can see that the real parts are largely similar to those of a full polytrope, but the imaginary parts are rather different from those of a full polytrope, showing enhanced tidal dissipation by introducing the rigid core. The rigid core model also supports inertial waves in the fluid envelope, but these waves have some peculiar behaviors due to the singularity in a spherical shell (Stewartson & Rickard 1969). Smooth inertial modes generally do not exist in a spherical shell even with uniform density (Rieutord et al. 2001), and localized wave beams spawned from the critical latitudes propagate in the bulk along the characteristics of the inertial wave equations (e.g., Ogilvie 2009). However, Lin & Ogilvie (2021) recently revealed that resonant tidal responses in a spherical shell correspond to eigen modes with large-scale flows hidden beneath localized wave beams using a uniform density model. Furthermore, it was shown that the hidden large-scale structures basically resemble inertial modes in a full sphere. This is in line with our results for the nonuniform density model in this study. The real parts, k₂₂ and k₄₂, are relevant to only large-scale density perturbations, which are similar to inertial modes in a full sphere as one can see from Fig. 5. Therefore, the curves of ∆k_1m for the rigid core model resemble those of a full polytrope, but we note slight shifts of the resonant frequencies due to the presence of a rigid core. As for the full polytrope, the compact rigid core model can produce ∆k₂₂ = −4% as observed, but it cannot produce sufficient dynamical correction in the high-degree Love number k42 near the tidal frequency of Io to account for the observed ∆k₄₂ = −11%.

On the other hand, the imaginary parts are largely modified by the presence of a small rigid core. We can see that the tidal dissipation is significantly enhanced by the localized wave beams spawned from the critical latitudes both in and out of resonances. The velocity perturbations in Fig. 5b indeed exhibit localized waves propagating in the bulk, which can generate significant viscous dissipation but they do not produce significant density and gravitational perturbations. In Fig. 4, we also see that several peaks in the tidal dissipation (bottom panel) do not lead to obvious fluctuations in ∆k_lm (top panel), corresponding to resonances with higher degree modes that contribute little to the low degree (i.e., l = 2 and l = 4) gravitational perturbations.

In summary for the compact rigid core model, the tidal dissipation is significantly enhanced with respect to the full polytrope case. This is in line with the early work of Ogilvie & Lin (2004), who showed the enhanced tidal dissipation due to inertial waves in the convective envelope of rotating stars and planets. The averaged dissipation in the tidal frequency range of Galilean moons gives rise to a comparable tidal quality factor as observed by Lainey et al. (2009). However, the fractional correction to the real part of Love number ∆k₄₂ is insufficient to explain the observation.

Fig. 4 As for Fig. 2, but for the interior model with a compact rigid core. The fractional correction ∆k42 (orange curve in the top panel) was multiplied by 0.07.

Fig. 5 Density perturbations (left half) and gravitational perturbations (right half) in the meridional plane to the tidal component ${\text{Ψ}}_4^2$ for the interior model with a compact rigid core at (a) ω/Ω = −1.5092 (resonance) and (b) ω/Ω = 1.53 (nonresonance) with Ek = 10−7. Amplitudes were normalized by the maximum absolute values.

Fig. 6 As for Fig. 2, but for the interior model with an extended dilute core. Green lines represent results without including the Coriolis force. The fractional correction ∆k42 (orange and green curves in the top panel) was multiplied by 0.07.

3.3 Dilute core model

An extended dilute core rather than a compact core in Jupiter has been suggested recently based on Juno gravitational measurements (Wahl et al. 2017; Militzer et al. 2022). In this subsection, we consider tidal responses for the interior model with an extended dilute core as shown in Fig. 1b. The dilute core was treated as a stably stratified layer which supports gravity waves restored by the buoyancy. If the Coriolis force is fully taken into account, dynamical tides in the dilute core region would be in the form of gravito-inertial waves (Dintrans et al. 1999; Xu & Lai 2017). Idini & Stevenson (2022b) recently calculated the tidal response of Jupiter with an extended dilute core, but they did not fully consider the Coriolis effect, which turns out to be important as we subsequently show.

Figure 6 shows the frequency dependence of the Love numbers for the dilute core model. For the tidal component Ψ²₂ (blue curves), the dynamical correction ∆k₂₂ is generally similar to that of the full polytrope, except for the absence of obvious spikes in the dilute core model. However, the imaginary part exhibits several peaks and troughs, suggesting possible resonances with high-degree mixed modes that enhance the tidal dissipation but do not significantly contribute to the l = 2 gravitational perturbations. The overall tidal dissipation is also enhanced with respect to the full polytrope due to the excitation of gravito-inertial waves in the dilute core and inertial waves in the convective envelope. The frequency-averaged tidal dissipation tends to be compatible with the observed tidal quality factor as we can see from Fig. 6.

For the tidal component ${\text{Ψ}}_4^2$, Fig. 6 also shows results without including the Coriolis force (green curves) for comparison. We note that the fractional correction ∆k₄₂ is always positive when the Coriolis force is neglected, probably because the pure gravity modes enhance the in-phase gravitational perturbations and thus produce positive dynamical corrections. Nevertheless, we observed distinct resonant responses at certain tidal frequencies from both real and imaginary parts of the Love number for the non-Coriolis case. For instance, the resonance at ω/Ω = −1.5193, which is close to the tidal frequency of Io, corresponds to the first gravity mode of l = 4 and m = 2 as shown in Fig. 7a. Indeed, Idini & Stevenson (2022b) propose the resonant locking between this gravity mode¹ (referred to as ${}_4^2{g}_1$) and the Jupiter-Io orbital evolution to explain the observed ∆k₄₂ for Jupiter. In Idini & Stevenson (2022b), the Coriolis force was neglected for the calculation of gravity modes, but approximated rotational corrections were made to obtain the Love number. However, fully taking the Coriolis force into account significantly alters the tidal responses as we can see from Fig. 6. The dynamical correction ∆k₄₂ exhibits several large fluctuations especially in the frequency range of −1.5 < ω/Ω < −1. This is due to the mixing of gravity modes and inertial modes in the dilute core, leading to more chances for resonances. The most significant dynamical corrections are produced near the tidal frequency ω/Ω = −1.2, which is close to the frequency of the purely inertial mode as shown in Fig. 3a. Of course, the inertial mode is mixed with gravity modes in the dilute core for this model. The resonance around ω/Ω = −1.2 can produce more than −10% dynamical corrections in k₄₂ (after the centrifugal correction), but it is too far away from the tidal frequency of Io. The resonance close to the tidal frequency of Io (also close to the frequency of pure gravity mode ${}_4^2{g}_1$) occurs at ω/Ω = −1.4448 when the Coriolis force is considered. Figure 7b shows the spatial structure of this resonant response. The Coriolis effect not only leads to a non-negligible shift in the mode frequency, but also largely modifies the mode structure. The perturbations are in the from of gravitoinertial waves in the dilute core and become pure inertial waves in the neutrally buoyant envelope. Non-negligible dynamical corrections are induced by this resonance at ω/Ω = −1.4448, but the corrections are insufficient (after the centrifugal correction) to account for the observed ∆k₄₂ = −11%. As the resonance is very narrow, we used 200 equally spaced frequency points in the tidal frequency interval of [−1.45, −1.43]. The peak amplitude of ∆k₄₂ in this frequency interval is comparable to the amplitude in the calculations using only 20 frequency points, suggesting that the frequency sampling points are sufficient to capture the resonant peak.

Comparing the orange and green curves in the bottom panel of Fig. 6, we can see that the tidal dissipation is increased by about two orders of magnitude when the Coriolis force is included. This suggests that the excitation of pure gravity waves is a less efficient tidal dissipation mechanism (unless resonances take place) based on our linear calculations, though the nonlinear interaction or wave breaking of gravity waves may lead to efficient tidal dissipation (e.g., Barker 2011; Weinberg et al. 2012).

Fig. 7 Density perturbations (left half) and radial velocity perturbations (right half) in the meridional plane to the tidal component ${\text{Ψ}}_4^2$ for the interior model with an extended dilute core. (a) Without including the Coriolis force at ω/Ω = −1.5193 (resonance); (b) including the Coriolis force at ω/Ω = −1.4448 (resonance). Amplitudes were normalized by the maximum absolute values.

3.4 Outer stable layer model

We finally consider the effect of an outer stable layer, which may exist in Jupiter resulting from H-He immiscibility (Debras & Chabrier 2019). Figure 8 shows the Love numbers as a function of the tidal frequency for the interior model c in Fig. 1, which includes a compact rigid core and a top stable layer between 0.8R and 0.9R. For the tidal responses to ${\text{Ψ}}_2^2$, the dynamical correction ∆k₂₂ is similar to the case without the stable layer, but the presence of the thin stable layer eliminates the spike due to the resonant inertial mode at the tidal frequency around ω/Ω = −1.08. The overall tidal dissipation due to ${\text{Ψ}}_2^2$ is comparable to the counterpart without the top stable layer (blue curve in the bottom panel of Fig. 4), but the fluctuation amplitudes (the differences between peaks and troughs) are smaller.

For the tidal responses to ${\text{Ψ}}_4^2$, we also show results for Ek = 10⁻⁷ (green curves) to illustrate the effect of fluid viscosity in Fig. 8. One can see that the viscosity has little influence on the real part of the Love number. The tidal dissipation weakly depends on viscosity at peaks and troughs, but the overall dissipation tends to be insensitive to viscosity. Indeed, Ogilvie (2013) has shown that the frequency-averaged dissipation is independent of viscosity.

The dynamical correction ∆k₄₂ is also similar to the case without the stable layer. We can see large variations in ∆k₄₂ at the tidal frequency around ω/Ω = −1.165, which corresponds to a resonant mode as shown in Fig. 9. This mode is complicated because it involves three different layers for the interior model considered here. The fluid body is primarily neutrally buoyant and supports inertial waves. However, the fluid domain is separated by the thin stable layer, which suppresses radial fluid motions and creates a “barrier” for the communication between inertial waves in the inner and outer regions (see the radial velocity and vorticity perturbations in Fig. 9). In addition, the thin stable layer supports rotationally modified gravity waves. The density perturbations are mainly restricted in the stable layer and the outer envelope (r > 0.8R). Despite the complicated velocity and density perturbations, the gravitational perturbations are dominated by the l = 4 component with relatively simple radial dependence. In this regard, this complicated mode is relevant to the l = 4 inertial mode without the stable layer, leading to large dynamical corrections around the tidal frequency at ω/Ω ≈ −1.1 as in Fig. 4. However, the dynamical correction ∆k₄₂ is negligible after the centrifugal correction at the tidal frequency of Io.

4 Conclusions

We have developed a numerical method for calculating the tidal responses of a compressible, self-gravitating, rotating, and viscous fluid body. We have fully taken the Coriolis force into account, but neglected the centrifugal distortion, which allowed us to solve the problem in the spherical geometry. We used the pseudo-spectral method based on spherical harmonics in the angular directions and Chebyshev collocation in the radial direction. Different from recent studies on Jupiter’s dynamical tides (Lai 2021; Idini & Stevenson 2022b; Dewberry & Lai 2022), we directly solved the tidally forced problem and explicitly added the fluid viscosity, which allowed us to simultaneously obtain the real and imaginary parts of the tidal Love numbers for a given planetary interior model.

In this study, we have considered three simplified interior models (Fig. 1) of Jupiter based on a polytrope of index 1. We have focused on the tidal components ${\text{Ψ}}_2^2$ and ${\text{Ψ}}_4^2$ in the frequency range of −2 ≤ ω/Ω ≤ −1, which is relevant to the tidal frequencies of Galilean moons. Our numerical results show that the dynamical correction ∆k₂₂ is generally insensitive to the interior models. All of the models we considered can give rise to the observed ∆k₂₂ ≈ −4% at the tidal frequency of Io, which is also in line with previous studies (Idini & Stevenson 2021; Lai 2021). The tidal dissipation is significantly enhanced by the presence of a compact rigid core model or an extended dilute core with respect to the full polytrope, leading to comparable tidal quality factor Q as observed by Lainey et al. (2009).

For the tidal responses to the ${\text{Ψ}}_4^2$ component, all of models we considered are difficult to give rise to ∆k₄₂ ≈ −11% near the tidal frequency of Io. For the interior model with a compact rigid core, significant dynamical corrections are generated at the tidal frequency around ω/Ω ≈ −1.1 due to the resonance with an inertial mode whose gravitational perturbations are dominated by the spherical harmonics of l = 4 and m = 2. However, this resonance is too far away from the tidal frequencies of Galilean moons. For the interior model with an extended dilute core, we demonstrate that the gravity modes in the dilute core can be significantly modified by the Coriolis force, leading to the mixed gravito-inertial modes. Resonances with gravito-inertial modes in the dilute core can produce non-negligible dynamical corrections, but they are insufficient to explain the observed ∆k₄₂ ≈ −11% near the tidal frequency of Io based on our simplified model. We also briefly investigated the effect of a top stable layer on Jupiter’s tides. The thin stable layer acts as a “barrier” and tends to restrict the density and velocity perturbations mainly in the outer envelope. However, our numerical results show that the top stable layer has little influence on the real part of tidal Love numbers.

As we have mentioned, we do not aim to construct a realistic interior model of Jupiter in this study. These simplified models were designed to characterize the tidal responses of some possible scenarios of Jupiter’s interior. Because the dynamical tides highly depend on the tidal frequency, the satellite-dependent tidal Love numbers would provide more constraints on the interior of Jupiter (Idini & Stevenson 2022b). In addition, seismology is the most effective approach to determine the interior structure of planets, though the detection of Jupiter’s oscillations remains a big challenge (Gaulme et al. 2011). Nevertheless, the numerical scheme we developed in this study can also be used for theoretical calculations of oscillation modes of giant planets.

There are some caveats, which should be considered in future. First, we did not consider the centrifugal deformation in order to solve the problem in the spherical geometry. The centrifugal effect plays a significant role in the tidal Love numbers of Jupiter, especially for the high-degree tidal components. Although we have made the centrifugal corrections when the numerical results were qualitatively compared with the observations, both the Coriolis and centrifugal effects should be self-consistently taken into account for quantitative comparisons with the high precision observations in the future. Second, giant planets exhibit differential rotations, which also influence the oscillation modes and thus tidal responses (Dewberry et al. 2021). Finally, Jupiter has the strongest magnetic field among planets in the Solar System and mainly consists of electrically conducting fluid (metallic hydrogen), so the magnetic effects (Lin & Ogilvie 2018; Wei 2022) should also play a part in the tides of Jupiter.

Fig. 8 As for Fig. 2, but for the interior model with a small rigid core and a top stably stratified layer. Green lines represent results at the Ekman number Ek = 10−7. The fractional correction ∆k42 (orange and green curves in the top panel) was multiplied by 0.07.

Fig. 9 Perturbations in the meridional plane to the tidal component ${\text{Ψ}}_4^2$ for the interior model c in Fig. 1. (a) Density (left half) and radial velocity (right half) perturbations; (b) gravitational (left half) and vorticity (right half) perturbations. Amplitudes were normalized by the maximum absolute values. The dashed lines denote r = 0.8R.

Acknowledgements

The author would like to thank an anonymous referee for constructive comments and Dali Kong for fruitful discussions. This study was supported by the B-type Strategic Priority Program of the CAS (XDB41000000), National Natural Science Foundation of China (grant no. 42174215) and the preresearch project on Civil Aerospace Technologies of CNSA (D020308). Numerical calculations were performed on the Taiyi cluster supported by the Center for Computational Science and Engineering of Southern University of Science and Technology.

References

All Figures

Fig. 1 Three nominal models of Jupiter’s interior used in this study. The top panel shows the schematic models and the bottom panel shows the density (normalized by the density at the center) and the Brunt-Väisälä frequency (normalized by the dynamical frequency) as a function of the radius. The blue shadow in the bottom panel indicates solid regions. (a) A compact rigid core model; (b) an extended dilute core model; and (c) a compact rigid core and an outer stable layer model.

In the text

Fig. 2 Complex Love number as a function of the tidal frequency for a full isentropic polytrope of index 1. Top panel: fractional correction ∆klm of the real part of the Love numbers. The fractional correction ∆k42 (orange curve in the top panel) was multiplied by 0.07. Bottom panel: minus imaginary part $-\mathrm{Im}\left[K_l^m\right]$, which is equivalent to klm/Ql. Vertical dashed lines indicate tidal frequencies of four Galilean moons of Jupiter (from right to left: Io, Europa, Ganymede, and Callisto). The horizontal dashed line in the bottom panel represents the astrometric observation of the frequency independent k2/Q2 from Lainey et al. (2009).

In the text

	Fig. 3 Density perturbations (left half) and radial velocity perturbations (right half) in the meridional plane to the tidal component ${\text{Ψ}}_4^2$ at two resonant frequencies in Fig. 2. Amplitudes were normalized by the maximum absolute values.
In the text

	Fig. 4 As for Fig. 2, but for the interior model with a compact rigid core. The fractional correction ∆k42 (orange curve in the top panel) was multiplied by 0.07.
In the text

	Fig. 5 Density perturbations (left half) and gravitational perturbations (right half) in the meridional plane to the tidal component ${\text{Ψ}}_4^2$ for the interior model with a compact rigid core at (a) ω/Ω = −1.5092 (resonance) and (b) ω/Ω = 1.53 (nonresonance) with Ek = 10−7. Amplitudes were normalized by the maximum absolute values.
In the text

	Fig. 6 As for Fig. 2, but for the interior model with an extended dilute core. Green lines represent results without including the Coriolis force. The fractional correction ∆k42 (orange and green curves in the top panel) was multiplied by 0.07.
In the text

	Fig. 7 Density perturbations (left half) and radial velocity perturbations (right half) in the meridional plane to the tidal component ${\text{Ψ}}_4^2$ for the interior model with an extended dilute core. (a) Without including the Coriolis force at ω/Ω = −1.5193 (resonance); (b) including the Coriolis force at ω/Ω = −1.4448 (resonance). Amplitudes were normalized by the maximum absolute values.
In the text

	Fig. 8 As for Fig. 2, but for the interior model with a small rigid core and a top stably stratified layer. Green lines represent results at the Ekman number Ek = 10−7. The fractional correction ∆k42 (orange and green curves in the top panel) was multiplied by 0.07.
In the text

	Fig. 9 Perturbations in the meridional plane to the tidal component ${\text{Ψ}}_4^2$ for the interior model c in Fig. 1. (a) Density (left half) and radial velocity (right half) perturbations; (b) gravitational (left half) and vorticity (right half) perturbations. Amplitudes were normalized by the maximum absolute values. The dashed lines denote r = 0.8R.
In the text