© The Authors 2025

1 Introduction

Tides are among the most fundamental mechanisms driving long-term celestial dynamics. They originate from the differential gravitational forces exerted by orbiting bodies on one another, leading to mass redistribution and energy dissipation. This process induces the transfer of angular momentum and energy between the spin of celestial bodies and their orbits, resulting in internal heating and progressive evolution of planetary systems over extensive timescales (e.g. Ogilvie 2014). In particular, tides shape the system architecture by establishing stable equilibrium states. The main stabilising effects include (i) circularisation, which tends to make orbits circular; (ii) synchronisation, whereby orbital and rotational periods become equal to each other; and (iii) alignment, whereby equatorial and orbital planes converge (e.g. Hut 1980, 1981).

Tidal dissipation also impacts the climate and surface conditions of rocky planets by influencing their capacity to maintain liquid water on their surfaces (see e.g. Bolmont 2018, and references therein). The associated tidal heating can transform the thermal state of such bodies, potentially melting interiors and triggering volcanism, as was seen on Io (Peale et al. 1979; Peale 2003). Similarly, Tyler (2008) showed that tides can produce substantial heat fluxes on icy moons in the outer Solar System, and thus sustain long-lived sub-surface oceans of liquid water (see also Tyler 2009, 2011, 2014, 2020). The dissipative processes generating tidal heating and the interplay between the oceanic response and the upper icy crust were subsequently examined by several authors (e.g. Matsuyama 2014; Kamata et al. 2015; Beuthe 2016; Matsuyama et al. 2018; Hay & Matsuyama 2019; Rekier et al. 2019; Rovira-Navarro et al. 2023).

Finally, tidal interactions are crucial to understanding the Earth-Moon system’s 4.5-billion-year climato-biologic evolution, especially through their effects on Earth’s length of day (LOD) and obliquity (e.g. Daher et al. 2021; Tyler 2021; Farhat et al. 2022a). On Earth, ocean tides dominate, accounting for over 90% of energy dissipation (e.g. Egbert & Ray 2000, 2001). As was proposed by Zahnle & Walker (1987), atmospheric thermal tides – namely, tides caused by Solar radiative flux – may have equally played an important role by partly counteracting the decelerating effect of gravitational tides on Earth’s rotation during the Precambrian. However, despite sparse geological data suggesting that Earth’s mid-Proterozoic LOD was stabilised (Bartlett & Stevenson 2016; Mitchell & Kirscher 2023; Wu et al. 2023), tidal locking likely never occurred (Farhat et al. 2024; Laskar et al. 2024), unlike on Venus, which is maintained in asynchronous rotation by thermal tides (Gold & Soter 1969; Ingersoll & Dobrovolskis 1978; Correia & Laskar 2001, 2003; Leconte et al. 2015).

Tides are usually considered as small perturbations around an equilibrium state, allowing for a linearised approach (e.g. Mathis et al. 2013). The linear response of a body to tidal forces is quantified by the potential Love number, named after A.E.H. Love, which relates the perturbed gravitational potential – resulting from variations in self-attraction – to the tide-raising potential (e.g. Ogilvie 2014). This potential Love number is a complex-valued, dimensionless transfer function that varies with tidal frequency and accounts for both tidal deformation and energy dissipation in the frequency domain¹. Commonly denoted by $k_l^m$, where l and m represent the degree and order of the spherical harmonic, respectively, the Love number is often simplified to a single parameter, k₂, or the degree-2 Love number, which depends solely on the tidal frequency and the internal physics of the tidally perturbed body (e.g. Mathis et al. 2013; Correia & Valente 2022; Valente & Correia 2022). Such an approximation implies rotational invariance around the body’s centre of mass, meaning that the background fields and tidal dynamics equations are independent of the body’s orientation relative to the perturber’s orbit. Throughout this study, we refer to this simplification as the ‘isotropic assumption’, noting that isotropy designates uniformity in all orientations. Although alternative terminology can be used, we opt for this one because it evokes a meaningful analogy between tidal perturbations and light rays propagating through an optical medium.

In solid bodies, the isotropic assumption is generally reasonable as a first-order approximation, since their internal structure and visco-elastic properties are primarily determined by radial self-attraction and pressure forces (e.g. Sotin et al. 2007). The resulting bodily tide is essentially a quasi-static adjustment to tidal forces, without inertial effects². Consequently, all tidal constituents are straightforwardly obtained from the equatorial degree-2 Love number describing the semidiurnal tidal bulge in the coplanar-circular configuration, which is just evaluated across several tidal frequencies – one per constituent. As a corollary, the degree-2 Love number is a symmetric function of the tidal frequency, with an even real part and an odd imaginary part, as was discussed by Efroimsky (2012a).

Unlike solid bodies, fluid layers exhibit anisotropic tidal responses – even in spherically symmetric cases – because fluid particle motion is influenced by Coriolis acceleration. By introducing a preferred direction (the spin axis), Coriolis acceleration distorts tidal flows, making the fluid response markedly different from that of solid materials. Notably, anisotropy induces scattering: the spherical harmonics coefficients describing the fluid tidal response are energetically coupled together instead of being independent as in solid bodies. This symmetry-breaking effect is intensified by the so-called dynamical tide – the component of the response associated with wave propagation (Zahn 1975) – since each oscillatory fluid mode can resonate at a specific frequency (e.g. Longuet-Higgins 1968; Savonije & Witte 2002; Ogilvie & Lin 2004; Ogilvie 2013, 2014; Auclair-Desrotour et al. 2018).

A wide variety of wave types can arise, mixing or remaining distinct depending on the fluid layer’s properties (e.g. Fuller et al. 2024). For instance, ocean tides on rocky planets are predominantly governed by long-wavelength surface gravity modes (Longuet-Higgins 1968; Cartwright 1977). These waves, restored by self-attraction, resemble the planetary-scale compressibility modes known as Lamb waves (e.g. Bretherton 1969; Lindzen & Blake 1972) that can be gravitationally or radiatively driven in planetary atmospheres, including Earth’s (e.g. Farhat et al. 2024). In other cases, the dynamical tide may consist of inertial, gravity, Alfvén, or acoustic waves, each restored by different forces – Coriolis acceleration, buoyancy (in stably stratified layers), magnetic tension, or pressure, respectively (e.g. Mathis et al. 2013; Rieutord 2015) – as is observed in stars and the fluid envelopes of giant planets.

Tidal anisotropy is further amplified by complex geometries. On Earth, continental coastlines redirect tidal flows, shifting resonant frequencies and disrupting symmetries (e.g. Green et al. 2017; Blackledge et al. 2020; Lyard et al. 2021; Auclair- Desrotour et al. 2023), while shallow seas add dynamic coupling effects between shelf and ocean tides (e.g. Brian K. Arbic & Garrett 2009; Arbic & Garrett 2010; Wilmes et al. 2023). Owing to these complexities, the role of anisotropic effects in planetary system evolution remains poorly understood. Yet, these effects may significantly alter the characteristic signatures of the dynamical tide, which manifests as staircase-like variations in orbital parameters during resonance episodes (e.g. Auclair- Desrotour et al. 2014; Tyler 2021; Motoyama et al. 2020; Farhat et al. 2022a; Huang et al. 2024; Zhou et al. 2024; Wu et al. 2024).

In the present study, we examine the limitations of the isotropic approximation to characterise their impact on the modelled orbital evolution of binary systems with misaligned spin and orbital angular momenta. Our approach extends previous research on the past history of the Earth-Moon system (Farhat et al. 2022a,b, 2024; Auclair-Desrotour et al. 2023). In Sect. 2, building on insights from Ogilvie (2013) and Boué (2017), we introduce a new formalism to calculate the time-averaged tidal rates of angular momentum and energy transfers in non-isotropic configurations. We further complement this formalism in Sect. 3 with the equations governing the long-term tidal evolution of a planet-perturber system. Finally, in Sect. 4, we apply these methods to an idealised Earth-Moon system, which allows us to probe the validity domain of the isotropic assumption and quantify related errors. Our findings reveal that the isotropic assumption results in notably inaccurate predictions for evolution models in the dynamical tide regime of fluid bodies, with errors scaling proportionally to the amplification factors associated with resonantly excited tidal waves. Notations and acronyms used throughout this paper are listed in the nomenclature in Appendix A.

2 Tidal torque and tidally dissipated power

2.1 Tidal torque as a function of tidal potentials

In this work, we examine the tidal interaction of two bodies orbiting around their centre of mass. For simplicity, we assume the central body to be a planet of mass M_p and radius R_p, and the tidal perturber a satellite of mass M_s. Nevertheless, the derivations presented in this section apply to any combination of central bodies and perturbers. The satellite induces a tidal potential, U_T, which disturbs the planet’s mass distribution. This, in turn, generates a deformation potential, U_D, which influences the system’s orbital evolution. Our initial aim is to establish the relationship between these two gravitational potentials and the rates of angular momentum and energy exchanges between the planet and satellite. These derivations are framed within linear theory, in which tides are treated as infinitesimal perturbations, and will later be used in Sect. 3 for numerical calculations.

We adopt a non-rotating, geocentric frame of reference, 𝓡_f: (O, e_x, e_y, e_z), centred at the planet’s centre of mass, O. The corresponding Cartesian unit vectors, (e_x, e_y, e_z), point towards distant stars. In 𝓡_f, the position of any point M is described using spherical polar coordinates (r, θ, φ), where r is the radial distance, θ the colatitude, and φ the longitude. The position vector of point M is denoted as r ≡ re_r, while the position vector of the satellite’s centre of mass is given by r_s = r_se_r, where e_r is the unit radial vector. Additionally, we define the rotating reference frame associated with the planet, 𝓡: (O, e_X, e_Y, e_Z), where e_Z aligns with the planet’s spin axis and points towards the north pole, while the other two unit vectors, e_X and e_Y, define orthogonal directions in the planet’s equatorial plane. In this rotating frame, the position of a point is described by the spherical coordinates $(\widehat{r},\widehat{\theta },\widehat{\phi })$.

As is illustrated by Fig. 1, the transformation between the two frames, 𝓡_f : (O, e_x, e_y, e_z) → 𝓡: (O, e_X, e_Y, e_Z), is represented by the 3-2-3 Euler rotation matrix R (α,β, γ) = R₃ (α) R₂ (β) R₃ (γ), where R₂ and R₃ are the rotation matrices around the second and the third axes, respectively. The matrix R is parametrised by the Euler angles (α, β, γ) given by (e.g. Varshalovich et al. 1988, Eq. (54) p. 30) $$R(\alpha ,\beta ,\gamma )\equiv \left[\begin{array}{ccc}{C}_{\alpha }{C}_{\beta }{C}_{\gamma }-S_{\alpha }{S}_{\gamma }& -C_{\alpha }{C}_{\beta }{S}_{\gamma }-S_{\alpha }{C}_{\gamma }& C_{\alpha }{S}_{\beta }\\ S_{\alpha }{C}_{\beta }{C}_{\gamma }+C_{\alpha }{S}_{\gamma }& -S_{\alpha }{C}_{\beta }{S}_{\gamma }+C_{\alpha }{C}_{\gamma }& S_{\alpha }{S}_{\beta }\\ -S_{\beta }{C}_{\gamma }& S_{\beta }{S}_{\gamma }& C_{\beta }\end{array}\right],$$(1)

with C_ϕ = cos ϕ and S_ϕ = sin ϕ for ϕ = α, β, or γ. Here, α and β denote the planet’s precession and nutation angles, respectively, while γ represents the intrinsic rotation angle. In Kaula’s theory, the angles α and β evolve over timescales that are much longer than typical tidal oscillation periods (e.g. Efroimsky & Williams 2009). Therefore, they are treated as constants in the formulation of the tidal problem. This is not the case for the intrinsic rotation angle, γ, which describes the planet’s spin, and explicitly varies with time, t, as $$\gamma (t)=\text{Ω}t,$$(2)

where Ω is the planet’s spin angular velocity.

The torque exerted by a force, F, about the origin point, O, on a planet occupying a volume, 𝒱, is expressed in spherical coordinates (r, θ, φ) as $$𝓣={\displaystyle {\int }_{\mathcal{V}}(r\times F)}\rho \text{d}V,$$(3)

where dV is an infinitesimal volume element, ρ is the local density, and × denotes the cross product. Following the convention of Zahn (1966), the tidal force is expressed as $$\mathit{F}(\mathit{r},t)=\nabla U_{\text{T}},$$(4)

where ∇ designates the gradient operator (see Appendix B). The tide-raising potential is given by (e.g. Kaula 1964; Efroimsky & Williams 2009) $$U_{\text{T}}\equiv G{M}_{\text{s}}\left(\frac{1}{\left|\mathit{r}-{\mathit{r}}_{\text{s}}\right|}-\frac{\mathit{r}\cdot {\mathit{r}}_{\text{s}}}{r_{\text{s}}^3}-\frac{1}{r_{\text{s}}}\right),$$(5)

with G = 6.67430 × 10⁻¹¹ m³ kg⁻¹ s⁻² being the universal gravitational constant (Tiesinga et al. 2021).

Since tidal oscillation periods are typically much shorter than the characteristic evolution timescales of the orbital system, we focus exclusively on the time-averaged tidal effects on the planet and satellite’s motions. We denote by 〈 f 〉 the average of any function of time, f, defined as $$\langle f\rangle \equiv \underset{P\to +\infty }{\lim }\frac{1}{P}{\displaystyle {\int }_0^{P}f}(t)\text{d}t.$$(6)

The time-averaged tidal torque exerted about the three directions of the Galilean frame (Rf ) is thus given by $$𝓣=\langle {\displaystyle {\int }_{\mathcal{V}}\left(ℒU_{\text{T}}\right)}\delta \rho \text{d}V\rangle ,$$(7)

where 𝓛 ≡ r × ∇ is the angular momentum operator (Varshalovich et al. 1988, Sect. 2.2), and δρ is the local tidal density variation. We can express δρ in terms of the gravitational potential of the tidally distorted body, U_D, using Poisson’s equation (e.g. Arfken & Weber 2005, Sect. 1.14), $${\nabla }^2{U}_{\text{D}}=-4\pi G\delta \rho ,$$(8)

with ∇² denoting the Laplacian operator, which leads to $$𝓣=-\frac{1}{4\pi G}\langle {\displaystyle {\int }_{\mathcal{V}}\left(ℒU_{\text{T}}\right)}\left({\nabla }^2{U}_{\text{D}}\right)\text{d}V\rangle .$$(9)

It is important to note that both U_T and U_D are expressed in the system of coordinates associated with the non-rotating frame of reference (𝓡_f), specifically (r, θ, φ).

Ogilvie (2013) demonstrated that the integral in Eq. (9) can, in fact, be carried out over any volume including the tidally distorted body. We therefore considered a simply connected domain, 𝒱_∗, which includes the planet but not the perturber, as is illustrated by Fig. 2, and we introduced the vector, A = (A_x, A_y, A_z), defined as $$\mathit{A}=\underset{{\mathcal{V}}_{*}}{{\displaystyle \int }}\left(ℒU_{\text{T}}\right)\left({\nabla }^2{U}_{\text{D}}\right)\text{d}V.$$(10)

Since 𝒱_* excludes the perturber, we have ∇²U_T = 0 in 𝒱_*. Furthermore, the angular momentum and Laplacian operators can be interchanged (e.g. Varshalovich et al. 1988, Sect. 2.2.2.), giving $${\nabla }^2\left(ℒU_{\text{T}}\right)=ℒ\left({\nabla }^2{U}_{\text{T}}\right)=0.$$(11)

As a result, the integral in Eq. (10) is rewritten as $$\mathit{A}={\displaystyle {\int }_{{\mathcal{V}}_{*}}\left[ℒU_{\text{T}}{\nabla }^2{U}_{\text{D}}-U_{\text{D}}{\nabla }^2\left(ℒU_{\text{T}}\right)\right]}\text{d}V.$$(12)

Thus, by virtue of Green’s second identity (e.g. Arfken & Weber 2005, Sect. 1.11), any component, A_v, of A, where v = x,y or z, is expressed as $$A_v={\displaystyle {\oint }_{\partial {\mathcal{V}}_{*}}\left[\left({ℒ}_v{U}_{\text{T}}\right)\nabla U_{\text{D}}-U_{\text{D}}\nabla \left({ℒ}_v{U}_{\text{T}}\right)\right]}\cdot \text{d}S.$$(13)

Here, dS is an outward-pointing infinitesimal surface element vector, and ∂𝒱_∗ is the boundary of 𝒱_∗. This allows the time- averaged torques exerted about the x, y, and z axes of 𝓡_f to be written as $${𝓣}_x=-\frac{1}{4\pi G}\langle {\displaystyle {\oint }_{\partial {\mathcal{V}}_{*}}\left[\left({ℒ}_x{U}_{\text{T}}\right)\nabla U_{\text{D}}-U_{\text{D}}\nabla \left({ℒ}_x{U}_{\text{T}}\right)\right]}\cdot \text{d}S\rangle ,$$(14) $${𝓣}_y=-\frac{1}{4\pi G}\langle {\displaystyle {\oint }_{\partial {\mathcal{V}}_{*}}\left[\left({ℒ}_y{U}_{\text{T}}\right)\nabla U_{\text{D}}-U_{\text{D}}\nabla \left({ℒ}_y{U}_{\text{T}}\right)\right]}\cdot \text{d}S\rangle ,$$(15) $${𝓣}_{ɀ}=-\frac{1}{4\pi G}\langle {\displaystyle {\oint }_{\partial {\mathcal{V}}_{*}}\left[\left({ℒ}_{ɀ}{U}_{\text{T}}\right)\nabla U_{\text{D}}-U_{\text{D}}\nabla \left({ℒ}_{ɀ}{U}_{\text{T}}\right)\right]}\cdot \text{d}S\rangle .$$(16)

Finally, the rate of energy transferred from the perturber’s orbital motion to the planet’s rotation is the time-averaged tidal power, 𝒫_T, which we express as a function of both the forcing and deformation tidal potentials too (e.g. Auclair-Desrotour et al. 2023), $${\mathcal{P}}_{\text{T}}=-\frac{1}{4\pi G}\langle \underset{{\mathcal{V}}_{*}}{{\displaystyle \int }}{U}_{\text{T}}{\nabla }^2\left({\partial }_t{U}_{\text{D}}\right)\text{d}V\rangle .$$(17)

By invoking Green’s second identity once again, we rewrite this expression as $${\mathcal{P}}_{\text{T}}=-\frac{1}{4\pi G}\langle \underset{\partial {\mathcal{V}}_{*}}{{\displaystyle \oint }}\left[U_{\text{T}}\nabla \left({\partial }_t{U}_{\text{D}}\right)-\left({\partial }_t{U}_{\text{D}}\right)\nabla U_{\text{T}}\right]\cdot \text{d}S\rangle .$$(18)

The time-averaged tidally dissipated power, 𝒫_diss, can be readily deduced from 𝒯 and 𝒫_T, as was demonstrated by Ogilvie (2013) and as is discussed in Sect. 3.2. Since Eqs. (14–16) and Eq. (18) provide the rates of momentum and energy exchanges as functions of the tide-raising and deformation gravitational potentials, U_T and U_D, we shall specify these two potentials in the following.

Fig. 1 Euler angles (α, β, γ) corresponding to the 3-2-3 Euler rotation matrix defined by Eq. (1), which describes the change of coordinate systems 𝓡f : O, ex, ey, ez → 𝓡: (O, eX, eY, eZ).

Fig. 2 Diagram illustrating the domain 𝒱∗, which includes the central body while excluding the perturber.

2.2 Tidal gravitational potential of the distorted planet

In linear theory, the tidal forcing and deformation potentials in Eqs. (14–16) and Eq. (18) are both series of periodic functions of time, since U_T is composed of oscillatory components. Consequently, it is appropriate to transition from the temporal to the frequency domain by expanding both potentials in Fourier series. These manipulations yield $$U_{\text{T}}(r,\theta ,\phi ,t)=\Re \left\{{\displaystyle \sum _{k=0}^{+\infty }{\tilde{U}}_{\text{T}}^k}(r,\theta ,\phi ){\text{e}}^{i{\sigma }_{k}t}\right\},$$(19) $$U_{\text{D}}(r,\theta ,\phi ,t)=\Re \left\{{\displaystyle \sum _{k=0}^{+\infty }{\displaystyle \sum _{q,p=-\infty }^{+\infty }{\tilde{U}}_{\text{D}}^{k,q,p}}}(r,\theta ,\phi ){\text{e}}^{i{\sigma }_{k,q,p}t}\right\},$$(20)

where ℜ denotes the real part of a complex number (with ℑ referring to the imaginary part), i represents the imaginary unit (i² = −1), σ_k is the k-th forcing tidal frequency, σ_k,q,p is the frequency of the tidal response defined by the triplet (k, q, p), $${\sigma }_{k,q,p}\equiv {\sigma }_k+(q-p)\text{Ω},$$(21)

and ${\tilde{U}}_{\text{T}}^k$ and ${\tilde{U}}_{\text{D}}^{k,q,p}$ designate the corresponding complex distributions for each potential. The frequencies, σ_k, are expressed as $${\sigma }_k=k{n}_{\text{s}},$$(22)

where n_s is the anomalistic mean motion of the perturber in the Galilean frame, and the integer, k, ranges from 0 to +∞³. It is noteworthy that the summations over q and p in Eqs. (19) and (20) result from the forward and backward rotations between the Galilean and rotating frames of reference, as will be discussed further.

Outside the planet, the coefficients, ${\tilde{U}}_{\text{T}}^k$ and ${\tilde{U}}_{\text{D}}^{k,q,p}$, are expressed as series of spherical harmonics, $Y_l^m$, of degree l and order m, as follows: $${\tilde{U}}_{\text{T}}^k(r,\theta ,\phi )={\displaystyle \sum _{l=2}^{+\infty }{\displaystyle \sum _{m=-l}^l{\tilde{U}}_{\text{T};l}^{k,m}}}{\left(\frac{r}{R_{\text{p}}}\right)}^l{Y}_l^m(\theta ,\phi ),$$(23) $${\tilde{U}}_{\text{D}}^{k,q,p}(r,\theta ,\phi )={\displaystyle \sum _{l=\left|p\right|}^{+\infty }{\displaystyle \sum _{m=-l}^l{\tilde{U}}_{\text{D};l}^{k,q,p,m}}}{\left(\frac{r}{R_{\text{p}}}\right)}^{-(l+1)}{Y}_l^m(\theta ,\phi ),$$(24)

where ${\tilde{U}}_{\text{T};l}^{k,m}$ and ${\tilde{U}}_{\text{p};l}^{k,q,p,m}$ are complex weighting coefficients associated with the triplet (k, l, m) and the quintuplet (k, l, m, p, q), respectively. The spherical harmonics, $Y_l^m$, in Eqs. (23) and (24) are explicitly defined in Appendix D.

The gravitational potential induced by the planet’s tidal response is usually obtained in the rotating frame of reference (𝓡), where the equations governing the planet’s tidal response are straightforwardly formulated. This requires expressing the tidal components given by Eqs. (23) and (24) as functions of the coordinates associated with 𝓡, namely $(\widehat{\theta },\widehat{\phi })$. This step is easily completed by considering the rotation operator of spherical harmonics: the spherical harmonics sharing the same degree in the Galilean and rotated systems of coordinates are linked through the relations (Varshalovich et al. 1988, Sect. 4.1) $$Y_l^m(\widehat{\theta },\widehat{\phi })={\displaystyle \sum _{q=-l}^l{D}_{q,m}^l}(\alpha ,\beta ,\gamma )Y_l^q(\theta ,\phi ),$$(25) $$Y_l^m(\theta ,\phi )={\displaystyle \sum _{q=-l}^l\overline{D_{m,q}^l}}(\alpha ,\beta ,\gamma )Y_l^q(\widehat{\theta },\widehat{\phi }).$$(26)

In the above equations, the notation $D_{q,m}^l$ refers to the Wigner D- functions (Varshalovich et al. 1988, Sect. 4.3, Eq. (1)), which are detailed in Appendix E, and $\overline{D_{m,q}^l}$ denotes the complex conjugate of $D_{m,q}^l$. Besides, we recall that γ is a function of time, as is defined in Eq. (2).

Using Eq. (25) in Eq. (23), we express the tideraising potential in the coordinates of the rotating frame of reference, $$U_{\text{T}}(\widehat{r},\widehat{\theta },\widehat{\phi },t)=\Re \left\{{\displaystyle \sum _{k,q}{\displaystyle \sum _{l=\max (|q|,2)}^{+\infty }{\widehat{U}}_{\text{T},l}^{k,q}}}{\left(\frac{\widehat{r}}{R_{\text{p}}}\right)}^l{Y}_l^q(\widehat{\theta },\widehat{\phi }){\text{e}}^{i{\widehat{\sigma }}_{k,q}t}\right\},$$(27)

where the integer k runs from 0 to +∞ and q from −∞ to +∞. In this expression, the tidal frequencies, ${\widehat{\sigma }}_{k,q}$, are defined as $${\widehat{\sigma }}_{k,q}\equiv q\text{Ω}+{\sigma }_k,$$(28)

and the complex weighting coefficients, ${\widehat{U}}_{\text{T};l}^{q,k}$, as $${\widehat{U}}_{\text{T};l}^{k,q}(\alpha ,\beta )={\displaystyle \sum _{m=-l}^l{\tilde{U}}_{\text{T};l}^{k,m}}{\text{e}}^{im\alpha }{d}_{m,q}^l(\beta ).$$(29)

The component of the forcing associated with the triplet (l, q, k) generates a tidal gravitational potential oscillating with the same frequency, but not necessarily with the same spatial distribution. The spatially dependent factor of this potential is expressed in general as $${\widehat{U}}_{\text{D};l}^{k,q}(\widehat{r},\widehat{\theta },\widehat{\phi })={\displaystyle \sum _{s=0}^{+\infty }{\displaystyle \sum _{p=-s}^s{H}_{l,s}^{k,q,p}}}{\widehat{U}}_{\text{T};l}^{k,q}{\left(\frac{\widehat{r}}{R_{\text{p}}}\right)}^{-(s+1)}{Y}_s^p(\widehat{\theta },\widehat{\phi }),$$(30)

with $H_{l,s}^{k,q,p}$ being the transfer function associated with the harmonic of degree s and order p. It is noteworthy that $H_{l,s}^{k,q,p}$ stands for a generalised Love number accounting for the fact that the spatial dependence of the tidal response generally differs from that of the tide-raising potential because of rotational scattering. For spherically isotropic bodies, this transfer function simplifies to $H_{l,s}^{k,q,p}=H_{l,l}^{k,q,q}{\delta }_{l,s}{\delta }_{q,p}$, where δ_q,p designates the Kronecker delta function, such that δ_q,p = 1 if q = p and δ_q,p = 0 otherwise.

In the coordinates of the rotating frame,, the gravitational potential of the tidally distorted body can be expressed as $$U_{\text{D}}(\widehat{r},\widehat{\theta },\widehat{\phi },t)=\Re \left\{{\displaystyle \sum _{k=0}^{+\infty }{\displaystyle \sum _{q=-\infty }^{+\infty }{\widehat{U}}_{\text{D}}^{k,q}}}(\widehat{r},\widehat{\theta },\widehat{\phi }){\text{e}}^{i{\widehat{\sigma }}_{k,q}t}\right\},$$(31)

where the spatial distributions associated with the frequencies, ${\widehat{\sigma }}_{k,q}$, are given by $${\widehat{U}}_{\text{D}}^{k,q}(\widehat{r},\widehat{\theta },\widehat{\phi })={\displaystyle \sum _{s=0}^{+\infty }{\displaystyle \sum _{p=-s}^s{\widehat{U}}_{\text{D};s}^{k,q,p}}}{\left(\frac{\widehat{r}}{R_{\text{p}}}\right)}^{-(s+1)}{Y}_s^p(\widehat{\theta },\widehat{\phi }).$$(32)

In this expression, the complex weighting coefficients, ${\widehat{U}}_{\text{D};s}^{k,q,p}$, are defined as $${\widehat{U}}_{\text{D};s}^{k,q,p}(\alpha ,\beta )={\displaystyle \sum _{l=\max (|q|,2)}^{+\infty }{H}_{l,s}^{k,q,p}}{\widehat{U}}_{\text{T};l}^{k,q}.$$(33)

By changing the coordinates from $(\widehat{r},\widehat{\theta },\widehat{\phi })$ to (r, θ, φ), we obtain the gravitational potential components ${\tilde{U}}_{\text{D};l}^{k,q,p,m}$, introduced in Eq. (24), first as functions of the ${\widehat{U}}_{\text{D};s}^{k,q,p}$ coefficients, $${\tilde{U}}_{\text{D};l}^{k,q,p,m}(\alpha ,\beta )={\widehat{U}}_{\text{D};l}^{k,q,p}{\text{e}}^{-im\alpha }{d}_{m,p}^l(\beta ),$$(34)

and second, using Eq. (33), in terms of the ${\widehat{U}}_{\text{T};s}^{k,q}$, $${\tilde{U}}_{\text{D};l}^{k,q,p,m}(\alpha ,\beta )={\displaystyle \sum _{s=\max (|q|,2)}^{+\infty }{H}_{s,l}^{k,q,p}}{\text{e}}^{-im\alpha }{d}_{m,p}^l(\beta ){\widehat{U}}_{\text{T};s}^{k,q}.$$(35)

Finally, substituting Eq. (29) into the above equation, we rewrite ${\tilde{U}}_{\text{D};l}^{k,q,p,m}$ as a function of the weighting coefficients of the forcing tidal potential in the inertial frame, 𝓡_f, $${\tilde{U}}_{\text{D};l}^{k,q,p,m}(\alpha ,\beta )={\displaystyle \sum _{s=\max (|q|,2)}^{+\infty }{\displaystyle \sum _{j=-s}^s{H}_{s,l}^{k,q,p}}}{\tilde{U}}_{\text{T};s}^{k,j}{\text{e}}^{i(j-m)\alpha }{d}_{m,p}^l(\beta )d_{j,q}^s(\beta ).$$(36)

This formulation of the deformation potential accounts for the possible coupling of spherical modes in the tidal response. When the central body is assumed to be spherically isotropic, no such coupling occurs. In that case, each component of the response corresponds to the same spherical harmonic as the forcing component generating it. Therefore, the summation over s in Eq. (36) vanishes and ${\tilde{U}}_{\text{D};l}^{k,q,p,m}$ simplifies to $${\tilde{U}}_{\text{D};l}^{k,q,p,m}(\alpha ,\beta )={\displaystyle \sum _{j=-l}^l{H}_{l,l}^{k,q,q}}{\tilde{U}}_{\text{T};l}^{k,j}{\text{e}}^{i(j-m)\alpha }{d}_{m,q}^l(\beta )d_{j,q}^l(\beta ){\delta }_{q,p}.$$(37)

2.3 Components of the tidal torque

As a final step, we express the tidal torque and tidally dissipated power, derived in Eqs. (14–16) and Eq. (18), in terms of the complex coefficients from the multipole expansions of the forcing and deformation tidal potential. These are the coefficients ${\tilde{U}}_{\text{T};l}^{k,m}$ and ${\tilde{U}}_{\text{D};l}^{k,q,p,m}$ introduced in Eqs. (23) and (24), respectively. This step is not straightforward, as the angular momentum operator introduces significant mathematical complexities. However, these difficulties can be elegantly addressed using angular momentum theory, as is demonstrated by Boué (2017). In this approach inspired from quantum mechanics, the problem is simplified by changing to a system of coordinates in which the action of the angular momentum operator on spherical harmonics is easily formulated. The terms involving the angular momentum operator in Eqs. (14–16) and Eq. (18) are then directly obtained by converting back to the Cartesian system, (x, y, z).

Following Boué (2017), we introduced the set of complex unit vectors (e₊, e₀, e₋), where the coordinates of any vector, a, are represented as (a₊, a₀, a₋). These are related to the Cartesian coordinates (a_x, a_y, a_z} in the reference frame, 𝓡_f, by $$a_{+}=-\frac{1}{\sqrt{2}}\left(a_x+i{a}_y\right),a_0=a_{ɀ},a_{-}=\frac{1}{\sqrt{2}}\left(a_x-i{a}_y\right).$$(38)

For v ∈ {−, 0, +} (i.e. v = −1,0,1, respectively), the v component, 𝓛_ν, of the angular momentum operator applied to a spherical harmonic, $Y_l^m$, is expressed as $${ℒ}_v{Y}_l^m=i{L}_{v,l}^m{Y}_l^{m+v},$$(39)

where the coefficients, $L_{v,l}^m$are real and defined as (Varshalovich et al. 1988) $$L_{v,l}^m\equiv \{\begin{array}{ll}m\hfill & \text{ if }v=0\hfill \\ -v\sqrt{\frac{l(l+1)-m(m+v)}{2}}\hfill & \text{ if }v=\pm 1.\hfill \end{array}$$(40)

The forcing and deformation tidal potentials, U_T and U_D, are written in terms of the complex potentials, Ũ_T and Ũ_D, as follows: $$U_{\text{T}}=\Re \left\{{\tilde{U}}_{\text{T}}\right\},U_{\text{D}}=\Re \left\{{\tilde{U}}_{\text{D}}\right\}.$$(41)

Since the real part and angular momentum operator can be interchanged, we have $$\mathbf{ℒ}\Re \left\{{\tilde{U}}_{\text{T}}\right\}=\Re \left\{\mathbf{ℒ}{\tilde{U}}_{\text{T}}\right\}.$$(42)

Using the change of coordinates introduced in Eq. (38), we obtain $${ℒ}_x\Re \left\{{\tilde{U}}_{\text{T}}\right\}=-\frac{1}{\sqrt{2}}\Re \left\{{\displaystyle \sum _{v=\pm 1}v}{ℒ}_v{\tilde{U}}_{\text{T}}\right\},$$(43) $${ℒ}_y\Re \left\{{\tilde{U}}_{\text{T}}\right\}=\frac{1}{\sqrt{2}}\Re \left\{i{\displaystyle \sum _{v=\pm 1}{ℒ}_v}{\tilde{U}}_{\text{T}}\right\},$$(44) $${ℒ}_{ɀ}\Re \left\{{\tilde{U}}_{\text{T}}\right\}=\Re \left\{{ℒ}_0{\tilde{U}}_{\text{T}}\right\}.$$(45)

To express the torques defined in Eqs. (14–16) in terms of the Fourier components of Ũ_T and Ũ_D, one must compute the conju- gates of 𝓛_xŨ_T, 𝓛_yŨ_T, and 𝓛_zŨ_T, as well as the radial component of their gradients and of the gradient of Ũ_D. The quantity, 𝓛_v Ũ_T, is given by $${ℒ}_v{\tilde{U}}_{\text{T}}={\displaystyle \sum _{k=0}^{+\infty }{\displaystyle \sum _{l=2}^{+\infty }{\displaystyle \sum _{m=-l}^{l}i}}}{L}_{v,l}^m{\tilde{U}}_{\text{T};l}^{k,m}{\left(\frac{r}{R_{\text{p}}}\right)}^l{Y}_l^{m+v}(\theta ,\phi ){\text{e}}^{i{\sigma }_{k}t}.$$(46)

This results in $$\overline{{ℒ}_x{\tilde{U}}_{\text{T}}}=\frac{i}{\sqrt{2}}{\displaystyle \sum _{k,v}{\displaystyle \sum _{l=2}^{+\infty }{\displaystyle \sum _{m=-l}^{l}v}}}{L}_{v,l}^m\overline{{\tilde{U}}_{\text{T};l}^{k,m}}{\left(\frac{r}{R_{\text{p}}}\right)}^l\overline{Y_l^{m+v}}(\theta ,\phi ){\text{e}}^{-i{\sigma }_{k}t},$$(47) $$\overline{{ℒ}_y{\tilde{U}}_{\text{T}}}=-\frac{1}{\sqrt{2}}{\displaystyle \sum _{k,v}{\displaystyle \sum _{l=2}^{+\infty }{\displaystyle \sum _{m=-l}^l{L}_{v,l}^m}}}\overline{{\tilde{U}}_{\text{T};l}^{k,m}}{\left(\frac{r}{R_{\text{p}}}\right)}^l\overline{Y_l^{m+v}}(\theta ,\phi ){\text{e}}^{-i{\sigma }_{k}t},$$(48) $$\overline{{ℒ}_z{\tilde{U}}_{\text{T}}}=-i{\displaystyle \sum _k{\displaystyle \sum _{l=2}^{+\infty }{\displaystyle \sum _{m=-l}^{l}m}}}\overline{{\tilde{U}}_{\text{T};l}^{k,m}}{\left(\frac{r}{R_{\text{p}}}\right)}^l\overline{Y_l^m}(\theta ,\phi ){\text{e}}^{-i{\sigma }_{k}t},$$(49)

with k ranging from 0 to +∞ and v = ±1.

Similarly, the radial component of the gradient of Ũ_D is given by $${\partial }_r{\tilde{U}}_{\text{D}}=-{\displaystyle \sum _{k,q,p}{\displaystyle \sum _{l=\left|p\right|}^{+\infty }{\displaystyle \sum _{m=-l}^l\frac{l+1}{r}}}}{\tilde{U}}_{\text{D};l}^{k,q,p,m}{\left(\frac{r}{R_{\text{p}}}\right)}^{-(l+1)}{Y}_l^m(\theta ,\phi ){\text{e}}^{i{\sigma }_{k,q,p}t},$$(50)

where k runs from 0 to +∞, and q and p from −∞ to +∞. However, only the components associated with the frequencies σ_k of the forcing tidal potential contribute to the tidal torque. These correspond to terms where p = q. Thus, the gradient of the tidal potential is rewritten as $$\begin{array}{l}{\partial }_r{\tilde{U}}_{\text{D}}=-{\displaystyle \sum _k{\displaystyle \sum _{l=0}^{+\infty }{\displaystyle \sum _{m=-l}^l\frac{l+1}{r}}}}\left({\displaystyle \sum _{q=-l}^l{\tilde{U}}_{\text{D};l}^{k,q,q,m}}\right){\left(\frac{r}{R_{\text{p}}}\right)}^{-(l+1)}{Y}_l^m(\theta ,\phi ){\text{e}}^{i{\sigma }_{k}t}\\ +{\partial }_r\mathcal{U},\end{array}$$(51)

where 𝒰 represents the component of Ũ_D that does not generate any torque in average.

The time-averaged components of the torque can be straightforwardly deduced from the preceding derivations. In essence, we consider the function f defined as $$f(t)={\displaystyle {\int }_{\mathcal{V}}A}(\mathit{r},t)B(\mathit{r},t)\text{\hspace{0.17em}d}V,$$(52)

where 𝒱 denotes the spatial domain over which the integral is evaluated, dV represents an infinitesimal volume element, and A and B are two functions expressed as $$A(\mathit{r},t)=\Re \{\tilde{A}(r,\theta ,\phi ,t)\},$$(53) $$B(\mathit{r},t)=\Re \{\tilde{B}(r,\theta ,\phi ,t)\}.$$(54)

In the above equations, the complex functions, A and B, are sinusoidal time-oscillating functions of positive frequencies, σ₁ and σ₂, respectively, multiplied by complex spatial distributions, $$\tilde{A}(r,\theta ,\phi ,t)={\tilde{A}}_0(r,\theta ,\phi ){\text{e}}^{i{\sigma }_{1}t},$$(55) $$\tilde{B}(r,\theta ,\phi ,t)={\tilde{B}}_0(r,\theta ,\phi ){\text{e}}^{i{\sigma }_{2}t}.$$(56)

As is demonstrated in Appendix C, the time-averaged value of f, defined by Eq. (6), can be expressed in terms of the complex spatial distributions Ã₀ and ${\tilde{B}}_0$ as $$\langle f\rangle =\{\begin{array}{ll}\Re \left\{\frac{1}{2}{\displaystyle {\int }_{\mathcal{V}}\overline{{\tilde{A}}_0}}{\tilde{B}}_0 \text{d}V\right\}\hfill & \text{ if }{\sigma }_1={\sigma }_2,\hfill \\ 0\hfill & \text{ if }{\sigma }_1\ne {\sigma }_2.\hfill \end{array}$$(57)

By applying this identity to the components of the tidal torque (Eqs. (14–16)) and tidal power (Eq. (18)), and following the steps detailed in Appendix F, we arrive at $${𝓣}_x=-\frac{K}{\sqrt{2}}\mathfrak{J}\left\{{\displaystyle \sum _{k,v}{\displaystyle \sum _{l=2}^{+\infty }{\displaystyle \sum _{m=-l}^{l}v}}}(2l+1)L_{v,l}^m\overline{{\tilde{U}}_{\text{T};l}^{k,m}}\left({\displaystyle \sum _{q=-l}^l{\tilde{U}}_{\text{D};l}^{k,q,q,m+v}}\right)\right\},$$(58) $${𝓣}_y=-\frac{K}{\sqrt{2}}\Re \left\{{\displaystyle \sum _{k,v}{\displaystyle \sum _{l=2}^{+\infty }{\displaystyle \sum _{m=-l}^l(2l+1)}}}{L}_{v,l}^m\overline{{\tilde{U}}_{\text{T},l}^{k,m}}\left({\displaystyle \sum _{q=-l}^l{\tilde{U}}_{\text{D};l}^{k,q,q,m+v}}\right)\right\},$$(59) $${𝓣}_{ɀ}=K\mathfrak{J}\left\{{\displaystyle \sum _k{\displaystyle \sum _{l=2}^{+\infty }{\displaystyle \sum _{m=-l}^l(2l+1)}}}m\overline{{\tilde{U}}_{\text{T},l}^{k,m}}\left({\displaystyle \sum _{q=-l}^l{\tilde{U}}_{\text{D};l}^{k,q,q,m}}\right)\right\},$$(60) $${𝒫}_{\text{T}}=-K\mathfrak{J}\left\{{\displaystyle \sum _k{\displaystyle \sum _{l=2}^{+\infty }{\displaystyle \sum _{m=-l}^l{\sigma }_k}}}(2l+1)\overline{{\tilde{U}}_{\text{T};l}^{k,m}}\left({\displaystyle \sum _{q=-l}^l{\tilde{U}}_{\text{D};l}^{k,q,q,m}}\right)\right\},$$(61)

where k runs from 0 to +∞, ν = ±1, and K is a constant defined as $$K=\frac{R_{\text{p}}}{8\pi G}.$$(62)

The formulae given by Eqs. (58–61) enable the calculation of the average rates of angular momentum and energy transfer in a general setting. Notably, these formulae remain valid even if the system does not conform to the spherical isotropy assumptions that typically prevail in the linear theory of bodily tides (see e.g. Efroimsky & Williams 2009, Sect. 4) or in the equilibrium tide model. Specifically, (i) the energy dissipation rate of one component is not solely dependent on the associated tidal frequency within this framework, and (ii) the deformation of the planet’s shape can differ spatially from the tide-generating potential, though it remains linearly related to it.

3 Dynamics of a planet-satellite system

3.1 Tide-raising gravitational potential

The tidal torque and power, as is expressed in Eqs. (58–60) and Eq. (61), determine the long-term evolution of planetary systems. To demonstrate how they specifically influence the spin and the orbital parameters of celestial bodies, we establish in this section the equations governing the tidal evolution of the Keplerian elements of a planet-perturber system. This is done using the formalism and method outlined in Boué & Efroimsky (2019).

We considered a two-body system in which a point-mass satellite orbits a planet that is not spherically isotropic. The satellite’s orbit is characterised by the standard Keplerian elements, as is illustrated by Fig. 3. We denote the semi-major axis with a_s, the eccentricity with e_s, the inclination with i_s, the true anomaly with υ_s, the longitude of the ascending node with ϑ_s, and the argument of the pericentre with ω_s. The reference frame associated with the satellite’s orbit, ℛ_s: (O, I_s, J_s, K_s), is centred on the planet’s centre of gravity and inclined relative to the inertial frame, ℛ_f. The unit vectors of this orbital frame are defined such that I_s points towards the pericentre of the orbit, K_s is aligned with the orbital angular momentum, and J_s = K_s × I_s. Additionally, the vector, n, points towards the ascending node, thereby defining the line of nodes. The transformation between the inertial frame, ℛ_f: (O, e_x, e_y, e_z), and the orbital frame, ℛ_s: (O, I_s, J_s, K_s), is described by an Euler rotation matrix, characterised by the angles (ɑ_s, β_s, γ_s), $${\alpha }_{\text{s}}\equiv {\vartheta }_{\text{s}}-\frac{\pi }{2},{\beta }_{\text{s}}\equiv i_{\text{s}},{\gamma }_{\text{s}}\equiv {\omega }_{\text{s}}+\frac{\pi }{2}.$$(63)

In this simplified framework, the forcing tidal potential can be expressed as a function of the mean anomaly of the perturber, ℳ_s, rather than the true anomaly. This is achieved using the Hansen coefficients, $X_k^{l,m}$, defined as (e.g. Hughes 1981; Laskar 2005) $${\left(\frac{r_{\text{s}}}{a_{\text{s}}}\right)}^l{\text{e}}^{im{v}_{\text{s}}}={\displaystyle \sum _{k=-\infty }^{\infty }{X}_k^{l,m}}\left(e_{\text{s}}\right){\text{e}}^{ik{ℳ}_{\text{s}}}.$$(64)

Here, the mean anomaly, ℳ_s, is defined as ℳ_s ≡ n_s t + ℳ_s;0, where ℳ_s;0 is a constant phase angle. After performing the manipulations detailed in Appendix G, we obtain $$\begin{array}{l}{\tilde{U}}_{\text{T};l}^{k,m}=\left(2-{\delta }_{k,0}\right)\frac{4\pi }{2l+1}\frac{G{M}_{\text{s}}}{R_{\text{p}}}{\left(\frac{R_{\text{p}}}{a_{\text{s}}}\right)}^{l+1}\\ \times {\displaystyle \sum _{q=-l}^l{D}_{m,q}^l}\left({\alpha }_{\text{s}},{\beta }_{\text{s}},{\gamma }_{\text{s}}\right)Y_l^q\left(\frac{\pi }{2},0\right)X_k^{-(l+1),-q}\left(e_{\text{s}}\right).\end{array}$$(65)

In this equation, we see the Wigner D-functions previously introduced in Eqs. (25) and (26). Notably, the factor $Y_l^q\left(\frac{\pi }{2},0\right)$ in Eq. (65) represents a specific value of the spherical harmonics (refer to Appendix D for details).

Fig. 3 Frame of reference ℛs: (O, Is, Js, Ks) and Keplerian elements describing the orbit of the perturber.

3.2 Dynamical equations

The planet’s spin angular velocity is typically much greater than the variation rates of the precession and nutation angles introduced in Eq. (1) as α and β, respectively. Mathematically, $\left|\dot{\alpha }\right|\ll \text{Ω}$ and $\left|\dot{\beta }\right|\ll \text{Ω}$. As a result, the gyroscopic approximation can be applied (see e.g. Boué & Efroimsky 2019, Sect. 2.8), which involves disregarding the contributions of the terms related to $\dot{\alpha }$ and $\dot{\beta }$ in the calculation of the planet’s angular momentum, L_p. Under this approximation, L_p is simply expressed as $${\mathit{L}}_{\text{p}}=C\text{Ω},$$(66)

where C designates the principal moment of inertia of the planet, and Ω ∈ Ωe_Z its spin vector. The evolution rate of the planet’s angular momentum is the tidal torque established in Eqs. (58–60), $$\stackrel{.}{{\mathit{L}}_{\text{p}}}=𝓣.$$(67)

Expanding $\dot{\text{Ω}}$ in terms of Ω, α and β, and using the basis unit vectors introduced in Sect. 2.1 and Fig. 1, we obtain $$\stackrel{.}{\text{Ω}}=\dot{\text{Ω}}{e}_Z+\text{Ω}\dot{\alpha }\left(e_{ɀ}\times e_Z\right)+\text{Ω}\dot{\beta }\left(e_y^{\prime }\times e_Z\right),$$(68)

which, owing to the orthogonality of the vectors, e_Z, e_z × e_Z, and e′_y × e_Z, yields the variation rates of Ω, α, and β, $$\dot{\text{Ω}}=\frac{𝓣\cdot e_Z}{C},\dot{\alpha }=\frac{𝓣\cdot \left(e_{ɀ}\times e_Z\right)}{C\text{Ω}{\sin }^2\beta },\dot{\beta }=\frac{𝓣\cdot \left(e_y^{\prime }\times e_Z\right)}{C\text{Ω}}.$$(69)

Analogously with L_p, the angular momentum of the satellite’s orbit, L_s, is defined as $${\mathit{L}}_{\text{s}}\equiv L_{\text{s}}{\mathit{K}}_{\text{s}},$$(70)

and its variation rate is expressed as $${\dot{\mathit{L}}}_{\text{s}}=-𝓣.$$(71)

Thus, expanding ${\dot{L}}_{\text{s}}$ in terms of the Euler angles describing the inclination of the orbit, α_s and β_s (see Eq. (63)), we end up with $${\stackrel{.}{\mathit{L}}}_{\text{s}}=\stackrel{.}{L_{\text{s}}}{\mathit{K}}_{\text{s}}+L_{\text{s}}{\dot{\alpha }}_{\text{s}}\left(e_{ɀ}\times {\mathit{K}}_{\text{s}}\right)+L_{\text{s}}{\dot{\beta }}_{\text{s}}\left(\mathit{n}\times {\mathit{K}}_{\text{s}}\right),$$(72)

and the variation rates of L_s, α_s, and β_s, $$\stackrel{.}{L_{\text{s}}}=-𝓣\cdot {\mathit{K}}_{\text{s}},\stackrel{.}{{\alpha }_{\text{s}}}=-\frac{𝓣\cdot \left(e_{ɀ}\times {\mathit{K}}_{\text{s}}\right)}{L_{\text{s}}{\sin }^2{\beta }_{\text{s}}},\stackrel{.}{{\beta }_{\text{s}}}=-\frac{𝓣\cdot \left(\mathit{n}\times {\mathit{K}}_{\text{s}}\right)}{L_{\text{s}}}.$$(73)

We remark that all the equations established until now are general and do not depend on the choice made for the fixed frame of reference (ℛ_f ).

It is actually appropriate to define ℛ_f so that e_z is aligned with the total angular momentum of the two-body system, $${\mathit{L}}_{\text{tot}}\equiv {\mathit{L}}_{\text{p}}+{\mathit{L}}_{\text{s}}=L_{\text{tot}}{e}_{ɀ},$$(74)

with the vectors e_x and e_y defining the so-called invariable plane; namely, the plane perpendicular to the total angular momentum (e.g. Tremaine et al. 2009; Rubincam 2016; Boué et al. 2016; Boué & Efroimsky 2019). In this configuration, α − α_s = π, which allows the evolution equations of the set of variables {Ω, L_s, β, β_s} to be deduced from Eqs. (69) and (73), $$\dot{\text{Ω}}=C^{-1}\left[\sin \beta \left(\cos \alpha {𝓣}_x+\sin \alpha {𝓣}_y\right)+\cos \beta {𝓣}_{ɀ}\right]$$(75) $$\stackrel{.}{L_{\text{s}}}=-\sin {\beta }_{\text{s}}\left(\cos {\alpha }_{\text{s}}{𝓣}_x+\sin {\alpha }_{\text{s}}{𝓣}_y\right)-\cos {\beta }_{\text{s}}{𝓣}_{ɀ},$$(76) $$\dot{\beta }={(C\text{Ω})}^{-1}\left[\cos \beta \left(\cos \alpha {𝓣}_x+\sin \alpha {𝓣}_y\right)-\sin \beta {𝓣}_{ɀ}\right],$$(77) $${\dot{\beta }}_{\text{s}}=L_{\text{s}}^{-1}\left[-\cos {\beta }_{\text{s}}\left(\cos {\alpha }_{\text{s}}{𝓣}_x+\sin {\alpha }_{\text{s}}{𝓣}_y\right)+\sin {\beta }_{\text{s}}{𝓣}_{ɀ}\right].$$(78)

The perturber’s semi-major axis and eccentricity can be obtained from energy conservation principles. The total energy of the system, E_tot, is formulated as $$E_{\text{tot}}=\frac{1}{2}C{\text{Ω}}^2+E_{\text{s}},$$(79)

where E_s designates the mechanical energy of the perturber’s motion; and its variation rate as $E_{\text{tot}}^{\cdot }=-{𝒫}_{\text{diss}}$. Since the power transferred from the planet’s spin to the planet’s orbital motion is expressed as $${\dot{E}}_{\text{s}}=-{\mathcal{P}}_T,$$(80)

taking the time derivative of Eq. (79) yields the relationship $$\begin{array}{l}{\mathcal{P}}_{\text{diss}}={\mathcal{P}}_{\text{T}}-\text{Ω}\cdot 𝓣,\\ ={\mathcal{P}}_{\text{T}}-C\text{Ω}\dot{\text{Ω}}.\end{array}$$(81)

The semi-major axis of the perturber is readily deduced from E_s (e.g. MacDonald 1964), $$a_{\text{s}}=-\frac{G{M}_{\text{s}}{M}_{\text{p}}}{2E_{\text{s}}}.$$(82)

Analogously, the angular momentum of the perturber is given by (e.g. MacDonald 1964) $$L_{\text{s}}=\frac{M_{\text{p}}{M}_{\text{s}}}{M_{\text{p}}+M_{\text{s}}}\sqrt{G\left(M_{\text{p}}+M_{\text{s}}\right)a_{\text{s}}{\left(1-e_{\text{s}}\right)}^2},$$(83)

which allows the eccentricity to be expressed as a function of the other quantities, $$e_{\text{s}}=1-\frac{\left(M_{\text{p}}+M_{\text{s}}\right)L_{\text{s}}^2}{{\left(M_{\text{p}}{M}_{\text{s}}\right)}^{2}G{a}_{\text{s}}}.$$(84)

We remark that the tangential component of L_tot is zero by construction, given that L_tot is aligned with e_z in the adopted frame of reference. Consequently, the tangential components of L_p and L_s exactly compensate each other, which is formulated as $$C\text{Ω}\sin \beta =L_{\text{s}}\sin {\beta }_{\text{s}}.$$(85)

The equatorial plane of the planet and the orbital plane of the perturber are thus inclined in two opposite directions, and the angle, β_s, can be deduced from β, Ω, and L_s, $${\beta }_{\text{s}}=\arcsin \left(\frac{C\text{Ω}}{L_{\text{s}}}\sin \beta \right).$$(86)

Similarly, L_tot appears to be a first integral of the problem, as no external torque acts on the system. The constancy of L_tot over time $\left(L_{\text{tot}}^{\cdot }=0\right)$ imposes an additional constraint on Ω, L_s, β, and β_s, expressed by the equation $$C\text{Ω}\cos \beta +L_{\text{s}}\cos {\beta }_{\text{s}}=L_{\text{tot}}.$$(87)

This constraint can be used to verify the consistency of the system’s dynamical evolution over time, as the left-hand side must remain unchanged despite variations in Ω, L_s, β, and β_s. The system’s evolution is then determined by integrating Eqs. (75–77) and Eq. (80) for {Ω, β, L_s, E_s}. The quantities a_s, e_s, and β_s are subsequently derived from the former using the relations provided in Eqs. (82), (84), and (86), respectively.

4 Application to the Earth–Moon system

4.1 Physical set-up of tidal solutions

In the previous sections, we derived the general expressions of tidal torque and power as functions of the forcing and deformation tidal potentials, incorporating the non-isotropic effects that are usually ignored. This formalism now enables us to investigate, as a next step, how deviations from the isotropic assumption may alter the estimates of these quantities and potentially lead to inaccurate predictions regarding the evolution of planetary systems.

For a didactical purpose, we applied the described methods to an idealised Earth-Moon system. The Earth was modelled as a two-layer, spherically symmetric body with a massive solid part and a thin, uniformly deep liquid water ocean, with a depth of H << R_p. The formalism used to compute the linear tidal response of such a planet has been extensively detailed in previous works (Auclair-Desrotour et al. 2018, 2019; Farhat et al. 2022b; Auclair-Desrotour et al. 2023). Therefore, we shall only briefly outline the key aspects of the approach followed here and direct readers to these studies for further details.

As is outlined in Auclair-Desrotour et al. (2023), the tidal response of the Earth’s solid part is modelled using Andrade rheology, with parameter values prescribed by Bolmont et al. (2020) for the actual Earth. Andrade rheology captures both viscoelastic and inelastic deformations of the solid regions. Although initially developed from lab experiments on metals (Andrade 1910), later studies demonstrated its applicability to silicate rocks and ices (see e.g. Efroimsky 2012b, and references therein). In this framework, the tidal response of the solid Earth is governed by the mantle’s effective shear modulus, µ, the Maxwell time, τ_M, which parametrises viscous friction, and two additional parameters linked to inelastic elongation: the Andrade time, τ_A, and a dimensionless parameter, α, which characterises the unrecoverable creep (e.g. Castelnau et al. 2008; Castillo-Rogez et al. 2011; Efroimsky 2012b). The latter parameter scales the frequency-dependent energy dissipation in the high-frequency regime.

The Andrade model is believed to accurately represent the tidal response of rocky bodies such as terrestrial planets and icy satellites, particularly in the high-frequency range where inelasticity dominates (e.g. Efroimsky & Lainey 2007). However, alternative models can also be applied (see e.g. Henning et al. 2009). The detailed computation of the solid Love numbers is presented in Auclair-Desrotour et al. (2023), Appendix H. This approach is based on Kelvin’s closed-form solution for an incompressible, homogeneous sphere (e.g. Munk & MacDonald 1960, Chapter 5). Here, we have employed this simplified model as a general method to describe the contribution of solid tides to energy dissipation in Earth-like scenarios, in which oceanic tides are predominant. Nevertheless, the calculation of the solid Love numbers can be refined later using more sophisticated methods based on elasto-gravitational theory (e.g. Takeuchi & Saito 1972; Tobie et al. 2005), along with realistic radial profiles of material properties (e.g. Sotin et al. 2007).

The oceanic tidal theory used here follows the approach detailed in Auclair-Desrotour et al. (2023). The ocean’s tidal response was obtained by solving the linear Laplace tidal equations (hereafter LTEs) over a spherical surface. Originally formulated by Laplace in his ‘Traité de Mécanique Céleste’ (Laplace 1798, Book IV), these equations were given a modern formulation in the foundational work by Longuet-Higgins (1968). We considered the LTEs in the shallow water approximation, in which tidal variables are depth-averaged (Vallis 2006). In this approach, the tidal solution reduces to barotropic tidal flows – those independent of vertical structure – which are dominant on Earth (Gerkema 2019). It represents planetary-scale surface gravity waves, restored by variations in self-attraction due to local ocean surface elevation (Vallis 2006). These waves are directly forced by the tide-raising potential, with their typical phase speed given by $\sqrt{𝑔H}$ (e.g. Cartwright 1977).

Consequently, vertical oceanic structure effects on tidal dynamics are neglected, meaning the contribution of internal gravity waves (i.e. waves restored by the Archimedean force; e.g. Gerkema et al. 2008) is formally excluded. However the impact of these waves is implicitly accounted for in the effective damping caused by dissipative processes. Using high-resolution TOPEX/Poseidon altimetric data, Egbert & Ray (2000) and Egbert & Ray (2001) demonstrated that energy dissipation on Earth primarily arises from bottom drag in shallow seas (~70%) and barotropic-to-baroclinic tidal conversion in deep oceans (25–35%). This conversion results from internal gravity waves generated by the interaction between tidal flows and ocean floor topography (Bell Jr 1975; Jayne & St. Laurent 2001; Garrett & Kunze 2007).

The LTEs consist of the horizontal momentum and continuity equations of fluid dynamics linearised for tidal perturbations, which are assumed to be small disturbances relative to the background fields (see e.g. Hay & Matsuyama 2019, for an exhaustive formulation of the LTEs including non-linear terms). These equations are expressed, respectively, as (e.g. Hendershott 1972; Matsuyama 2014; Matsuyama et al. 2018; Auclair-Desrotour et al. 2023) $${\partial }_{t}V+{\sigma }_{\text{R}}V+\mathit{f}\times V+𝑔\nabla \left({\Gamma }_{\text{D}}\zeta -{\Gamma }_{\text{T}}{\zeta }_{\text{eq}}\right)=0,$$(88) $${\partial }_t\zeta +\mathit{\nabla }\cdot (H\mathit{V})=0,$$(89)

where $𝑔\equiv G{M}_{\text{p}}/R_{\text{p}}^2$ is the planet’s surface gravity, σ_R is the Rayleigh drag frequency accounting for the efficiency of dissipative processes, and $f\equiv 2\text{Ω}\cos \widehat{\theta }{e}_r$ denotes the Coriolis parameter, with e_r being the outward-pointing radial unit vector. Here, V represents the horizontal velocity field, ζ the vertical displacement of the ocean’s surface relative to the ocean floor, and ζ_eq = U_T/𝑔 the equilibrium displacement caused by the tidal gravitational potential in the zero-frequency limit.

For simplification, mean flows are ignored in Eq. (88), treating the ocean as a static fluid layer. This assumption is based on the separation of spatial scales and timescales between tidal and mean flows, with the latter characterised by smaller spatial structures and longer timescales, reducing the likelihood of significant coupling. However, this assumption does not always hold for thick fluid envelopes, where large-scale circulation patterns, such as differential rotation or latitudinal shear, can affect tidal flows by inducing Doppler-shift effects (e.g. Boyd 1978; Ortland 2005a,b; Baruteau & Rieutord 2013; Guenel et al. 2016). While permanent oceanic currents are generally thought to have minimal impact on tidal dynamics, they are influenced by tidal energy dissipation and internal waves, which interact with oceanic mesoscale eddies (see e.g. Arbic 2022, and references therein).

The Rayleigh drag frequency, σ_R, plays a crucial role as it governs the frequency-dependent rates of energy and momentum exchanges within the planet-perturber system. As σ_R decreases, both energy and momentum transfer rates become more sensitive to tidal frequency (e.g. Auclair-Desrotour et al. 2018). This behaviour is linked to the phenomenon of the dynamical tide, which refers to the component formed by forced tidal waves propagating through a fluid layer (Zahn 1975). In this context, the dynamical tide results from long-wavelength surface gravity waves, which can be resonantly excited by the tide-raising potential, leading to significant variations in the tidally dissipated energy (e.g. Webb 1980; Tyler 2014, 2021; Kamata et al. 2015; Auclair-Desrotour et al. 2018, 2019, 2023; Matsuyama et al. 2018; Hay & Matsuyama 2019; Farhat et al. 2022b; Rovira-Navarro et al. 2023). As σ_R diminishes, the strength of the dynamical tide increases, with the height of resonant peaks scaling as $\propto {\sigma }_{\text{R}}^{-1}$ in linear theory (see e.g. Auclair-Desrotour et al. 2019, Eq. (58)).

This relationship generally holds for fluid bodies. The frequency dependence of the tidal response driven by forced wave propagation in planetary fluid layers and stars becomes more pronounced as the damping coefficient for frictional forces (or other dissipative mechanisms) decreases. This can result in energy dissipation rates that span several orders of magnitudes in stars and giants planets (e.g. Savonije & Witte 2002; Ogilvie & Lin 2004, 2007; Ogilvie 2009, 2013; Auclair Desrotour et al. 2015; Mathis et al. 2016; Astoul & Barker 2023), or in the sub-surface oceans of icy moons (Tyler 2011, 2014; Chen et al. 2014; Matsuyama 2014; Matsuyama et al. 2018; Beuthe 2016). Rayleigh drag frequency is typically unconstrained and can vary widely. For instance, Matsuyama et al. (2018) report values as low as σ_R ~ 10⁻¹¹ s⁻¹ for icy moons in the Solar System, with σ_R ~ 10⁻⁹ s⁻¹ for obliquity tides on Europa. In the present study, we adopt σ_R = 10⁻⁵ s⁻¹, in line with estimates reckoned from high-precision measurements of Earth’s tidal energy dissipation and advanced tidal models (e.g. Webb 1980; Egbert & Ray 2001, 2003; Farhat et al. 2022b). Given the strong sensitivity of tidal dissipation to σ_R, this value represents a conservative estimate for examining the impact of anisotropy on tidal dissipation. Any smaller value would amplify the frequency dependence of the oceanic tidal response.

The operators, Γ_D and Γ_T, in Eq. (88) account for ocean loading, self-attraction, and the deformation of the planet’s solid regions, which slightly modify the oceanic tidal response (e.g. Hendershott 1972; Ray 1998; Egbert et al. 2004). These operators are expressed in the frequency domain as functions of the solid tidal and load Love numbers, as is detailed in Auclair-Desrotour et al. (2023). Also, it is important to note that the three-dimensional gradient operator, ∇, introduced in Eq. (7), is applied over the spherical surface in Eq. (88), along with the divergent operator, ∇· in Eq. (89). These operators are explicitly formulated in spherical coordinates in Appendix B. Equations (88) and (89) are solved for the velocity field, V, and ocean surface displacement, ζ, in the frequency domain for every component of the tide-raising gravitational potential expressed in the rotating frame of reference, as is provided by Eq. (27). This potential is computed using Eq. (65), with the degree-2 truncation (l ≤ 2) corresponding to the classical quadrupolar approximation (e.g. Mathis et al. 2013, Sect. 7.2.2.5).

Following the approach of Auclair-Desrotour et al. (2023), we numerically solved the LTEs by employing a spectral method that expands the tidal quantities in spherical harmonics series. Specifically, using the notations introduced in Sect. 2.2, the variations in ocean depth are expressed as $$\zeta (\widehat{\theta },\widehat{\phi },t)=\mathit{\Re }\left\{{\displaystyle \sum _{k=0}^{+\infty }{\displaystyle \sum _{q=-\infty }^{+\infty }{\widehat{\zeta }}^{k,q}}}(\widehat{\theta },\widehat{\phi }){\text{e}}^{i{\widehat{\sigma }}_{k,q}t}\right\}.$$(90)

The spatial functions ${\widehat{\zeta }}^{k,q}(\widehat{\theta },\widehat{\phi })$ are further expanded as $${\widehat{\zeta }}^{k,q}(\widehat{\theta },\widehat{\phi })={\displaystyle \sum _{s=0}^N{\displaystyle \sum _{p=-s}^s{\widehat{\zeta }}_s^{k,q,p}}}{Y}_s^p(\widehat{\theta },\widehat{\phi }),$$(91)

where ${\widehat{\zeta }}_s^{k,q,p}$ are the complex coefficients of the multipole expansion and N denotes the truncation degree of the series. As discussed in Appendix H from a convergence analysis, the truncation degree in Eq. (91) is set to N = 30. This value empirically appears to be sufficiently large to account for the complex mapping between the ocean’s forced modes and spherical harmonics.

Finally, from the solution, we derived the Fourier coefficients of the deformation tidal potential needed to evaluate the tidal torque and power, as is formulated in Eqs. (58–61). These coefficients, introduced in Eq. (32), are expressed in terms of the components of the tide-raising potential and the variation in ocean depth, as is shown in (e.g. Auclair-Desrotour et al. 2023, Eq. (72))⁴, $${\widehat{U}}_{\text{D};s}^{k,q,p}=k_s^{{\widehat{\sigma }}_{k,q}}{\widehat{U}}_{\text{T};s}^{k,q}+\left(1+k_{\text{L};s}^{{\widehat{\sigma }}_{k,q}}\right)\frac{3}{2s+1}\frac{{\rho }_{\text{w}}}{\overline{\rho }}\mathit{𝑔}{\widehat{\zeta }}_s^{k,q,p},$$(92)

where $\overline{\rho }\equiv 3M_{\text{p}}/\left(4\pi R_{\text{p}}^3\right)$ denotes the planet’s mean density, ρ_w ≈ 1035 kg m⁻³ represents the average density of seawater (e.g. Ray et al. 2001), and $k_s^{{\widehat{\sigma }}_{k.q}}$ and $k_{\text{L};s}^{{\widehat{\partial }}_{k,q}}$ are the complex solid tidal and load Love numbers that describe the visco-elastic deformation of solid regions at the frequency ${\widehat{\sigma }}_{k.q}$. The first term on the righthand side of Eq. (92) accounts for the solid elongation caused by the forcing tidal potential. The second term captures the contribution of oceanic tides to self-attraction variations, consisting of two components: the gravitational potential created by surface elevation changes and the deformation of solid regions due to ocean loading.

4.2 Sensitivity of tidal dissipation to planet’s spin rotation and obliquity

In our calculations, we assume that the perturber causing the tidal forcing is a Lunar-mass satellite in a circular orbit around the planet (e_s = 0), with the same semi-major axis as the current Earth-Moon system, and a fixed orbital period, P_s. The planet’s rotation period, P_rot, is varied from 0.1 to 100 days, and its obliquity, β, from 0° to 90°. We investigated two configurations. In the first configuration (‘DRY’), the ocean was neglected, and only the solid tide was considered. In the second configuration (‘GLOBAL OCEAN’), the full tidal response, including the ocean, was evaluated. For each configuration, the calculations were performed in two scenarios: (i) using the full formalism introduced in Sect. 2, where the planet’s tidal response is computed comprehensively and self-consistently (‘FULL’), and (ii) under the classical isotropic approximation, where all tidal components are derived from the equatorial degree-2 Love number defined in the circular-coplanar configuration for the semidiurnal tide (’APPROX’). We computed the tidal solutions using the numerical tools of the TRIP software (Gastineau & Laskar 2011)⁵.

In the APPROX case, the equatorial degree-2 Love number was calculated for every tidal frequency of the forcing potential in the planet’s rotating frame. Each tidal frequency, σ, was treated as a semidiurnal tidal frequency, σ = 2 (Ω − n_s), where the fictitious orbital frequency is defined as n_s = Ω − σ/2. The corresponding tidal response was then evaluated using the spatial distribution for the semidiurnal quadrupolar tidal potential, represented by the spherical harmonic of degree 2 and order 2 (denoted by $Y_2^2$). From the resulting solution, the degree-2 Love number was derived and used to compute the deformation potential for the relevant tidal component. This potential is simply the product of the degree-2 Love number and the tide-raising potential, following the theory of bodily tides. It is important to note that the Love number depends solely on σ in the dry configuration due to the spherical isotropy of the solid part. However, in the global ocean configuration, it also depends on the planet’s spin period because of the influence of Coriolis forces, as was described earlier.