© The Authors 2023

1 Introduction

‘Can one hear the shape of a drum?’ This question, raised by Kac (1966), emphasises the link connecting the geometry of a vibrating system to the frequencies at which it can vibrate. Mathematically, it defines a broad class of problems referred to as spectral geometry, which aims to establish relationships between the eigenmodes of Riemannian manifolds and their geometric features. In the aforementioned study, Kac investigates whether it is possible or not to infer some information about the shape of a vibrating drumhead from the sound it makes, which amounts to recovering the drumhead’s geometry from its acoustic signature in an inverse problem approach. Unfortunately, the answer to this question was shown to be negative in the general case, as two different shapes are able to produce the very same sound (e.g. Gordon et al. 1992). Nevertheless, it remains possible to infer some specific geometric features of the system to a certain extent.

Interestingly, this statement is also applicable to tidally forced ocean planets, which may be regarded as giant celestial vibrating drumheads. In the general case, one calls ‘ocean planets’ rocky bodies that are partly or totally covered by a surface liquid layer¹ (e.g. Léger et al. 2004). With more than 70% of its surface covered by oceans (e.g. Eakins & Sharman 2010), the Earth thus appears as a typical ocean planet. It is actually the only planet of this kind that can be found presently in the Solar System, although observational evidences indicate that Mars possibly had large paleo-oceans as well before it dried (Carr & Head 2003; Dohm et al. 2009; Scheller et al. 2021).

This uniqueness of the Earth, however, diminishes if one considers extrasolar systems. The remarkable exoplanetary diversity established thus far suggests that ocean planets are rather ordinary, if not widespread. Particularly, a substantial number of rocky planets were found to orbit in the habitable zone of their host stars, which is the region where temperature conditions can sustain liquid water on the planets’ surfaces (Kasting et al. 1993). This is the case of many rocky planets orbiting red and brown dwarfs (e.g. Payne & Lodato 2007; Raymond et al. 2007; Kopparapu et al. 2017) such as the super-Earths GJ 1214 b (e.g. Charbonneau et al. 2009), LHS 1140 b (e.g. Lillo-Box et al. 2020), and TOI-1452 b (Cadieux et al. 2022), or the Earth-sized planets TRAPPIST-1 e, f, and g (e.g. Grimm et al. 2018), with the latter being suspected to harbour ocean-scale volumes of liquid water (e.g. Bolmont et al. 2017; Bourrier et al. 2017; Turbet et al. 2018). Improving our understanding of the surface conditions, climate, and fate of these planets requires to constrain their long-term orbital evolution better (e.g. Pierrehumbert 2010).

Tides are the main mechanism that controls the evolution of planetary systems over timescales ranging from millions to billions of years (e.g. Hut 1981; Correia & Laskar 2001; Levrard et al. 2009). Apart from the atmospheric thermal tides that are induced by insolation variations (e.g. Lindzen & Chapman 1969; Leconte et al. 2015; Auclair-Desrotour et al. 2019b), they result from the mutual gravitational interactions between celestial bodies, which undergo the differential attraction of their neighbours. Combined with dissipative processes, this tidal forcing generates a delayed mass redistribution, thereby leading to exchanges of angular momentum between the orbit of the tidal perturber and the spin of the tidally forced body (e.g. Correia et al. 2014). Additionally, tides are accompanied with energy dissipation where the mechanical energy lost by the orbital system is converted into heat. While it is negligible on Earth, this tidal heating is sometimes able to partly melt the body and to generate surface volcanism, as observed on the Jovian moon Io (e.g. Peale et al. 1979).

Oceanic tides take the form of frequency-resonant surface gravity modes distorted by the planet’s rotation (e.g. Auclair-Desrotour et al. 2018), which strongly differs from the smooth visco-elastic elongation of solid bodies (e.g. Correia et al. 2014). As a consequence, tidal dissipation in oceans may be increased by several orders of magnitude while crossing a resonance associated with a predominant mode (e.g. Arbic & Garrett 2010). This resonant behaviour notably explains why, for present day Earth, the oceanic contribution (~2.5 TW) to the semidiurnal component (M₂) of dissipated tidal energy is an order of magnitude greater than the solid counterpart (Lambeck 1977; Provost & Lyard 1997; Egbert & Ray 2001; Tyler 2021). Analogously, oceanic tidal flows presumably warm up the icy moons of the Solar System outer planets that are suspected to harbour subsurface oceans, such as Europa, Callisto, or Titan (e.g. Tyler 2008, 2014; Kamata et al. 2015; Beuthe 2016; Matsuyama et al. 2018).

Starting with MacDonald (1964), several authors investigated the role played by oceanic tides in the evolution of the Earth-Moon system, using either semi-analytical approaches with simplified geometries (e.g. Webb 1980, 1982; Bills & Ray 1999; Farhat et al. 2022a) - including global ocean models (e.g. Auclair-Desrotour et al. 2018; Motoyama et al. 2020; Tyler 2021) –, or numerical methods based upon realistic land-ocean distributions (e.g. Le Provost et al. 1994, 1998; Arbic et al. 2010; Kodaira et al. 2016; Green et al. 2017; Blackledge et al. 2020; Daher et al. 2021). For a didactical review of the pioneering studies of the field, the reader is referred to Lambeck (1977) and Bills & Ray (1999). This long series of works particularly emphasises the strong interplay between the resonant oceanic tidal flows and continents, which shapes the four billion year-history of the Earth’s length of the day (LOD) and Earth-Moon distance (see e.g. Daher et al. 2021; Farhat et al. 2022a). Moreover, Green et al. (2017) show that the topography significantly increases tidal dissipation, and Blackledge et al. (2020) underline the sensitivity of the latter to coastlines’ fractality. However, the way the large-scale geometry of the ocean basin alters the planet’s tidal response cannot be easily disentangled from other effects in realistic models due to its complexity.

The present work is an attempt to address this issue by generalising the semi-analytical tidal theory of hemispherical oceans (Longuet-Higgins & Pond 1970; Webb 1980, 1982; Farhat et al. 2022a) to ocean basins of arbitrary sizes. In this approach, the oceanic depth is assumed to be uniform and thin compared to the planet radius. Adopting the so-called ‘shallow water’ approximation (e.g. Vallis 2006), we ignore both the ocean stratification and the associated baroclinic component of tidal flows - that is the contribution of internal gravity waves (e.g. Gerkema & Zimmerman 2008) – so that the described tide is purely barotropic. Also, the geometry of the continent is simplified to a spherical cap of specified angular radius in order to avoid mathematical complications. Finally, bottom friction is modelled by a standard Rayleigh drag, and the coupling effect between the ocean and the deformable solid surface is taken into account. Analogously with the examples of the global and hemispherical oceans, this formalism allows the planet’s tidal response and the resulting dissipated energy to be formulated in terms of explicit eigenmodes, each of these eigenmodes being resonant for a specific tidal frequency. Additionally, the ocean dynamics is controlled by a small number of dimensionless parameters, which provides an appropriate framework for probing the parameter space and characterising the continentality effect.

In Sect. 2, we detail the theory by successively introducing the Laplace’s Tidal Equations (LTEs) that govern the oceanic tidal response, the eigenmodes of the ocean basin, and the expressions of the tidal Love numbers, torque, and power. In Sect. 3, we examine the frequency-behaviour of an Earth-sized planet with an equatorial hemispherical ocean. This reference case is used in Sect. 4 to characterise the sensitivity of the planet’s tidal response to the position of the continent on the globe, its size, the oceanic depth, the dissipation timescale, the elasticity of the solid part, and its anelasticity timescale. We elaborate on the continentality effect in Sect. 5 by introducing a metric that quantifies the dependence of the tidal torque on the position and size of the supercontinent. Finally, in Sect. 6, the conclusions of the study are summarised. We stress here that Sect. 2 details technical aspects of the developed model. Thus we invite the reader that might be only interested in the application of the theory to skip this section and to jump directly to Sect. 3. It is also noteworthy that all the notations introduced in the main text can be retrieved in the nomenclature made in Appendix A.

2 Oceanic tidal response

We establish here the equations governing the tidal dynamics of the ocean basin harboured by a rocky planet of radius R_p. These equations are based on the commonly used shallow water approach, where the ocean is considered as a thin liquid layer of uniform density and depth H ≪ R_p (e.g. Vallis 2006). In this approach, the fluid is supposed to be incompressible. The tidal response thus described is said to be barotropic because it does not depend on the vertical structure of the ocean (Vallis 2006). The continental geometry is simplified to a spherical cap of angular radius θ_c, as shown by Fig. 1. Accordingly, the angular radius of the ocean basin is defined as θ₀ ≡ π − θ_c.

2.1 Geometry of the ocean basin

To write down the equations that describe the oceanic tidal waves, we shall preliminarily introduce two frames of references and their associated systems of coordinates. First, we denote by ℛ: (O, e_X, e_Y, e_Z) the frame of reference rotating with the planet and having its centre of gravity, O, as origin. The Cartesian basis of unit vectors (e_X, e_Y, e_Z) is such that e_X and e_Y correspond to two orthogonal directions in the planet’s equatorial plane, and e_Z to the direction of the spin vector. The unit vector e_Z thus designates the position of the north pole on the unit sphere. The geocentric frame of reference ℛ is associated with the usual spherical coordinates ($\widehat{r}$, $\widehat{\theta }$, $\widehat{\phi }$), where $\widehat{r}$, $\widehat{\theta }$, and $\widehat{\phi }$ designate the radial, colatitudinal and longitudinal coordinates, respectively. The position of the circular supercontinent on the globe is defined by the coordinates of its centre, (${\widehat{\theta }}_{\text{S}}$, ${\widehat{\phi }}_{\text{S}}$). As the centre of the oceanic basin corresponds to the antipodal point, its coordinates on the globe are given by $\left({\widehat{\theta }}_{\text{oc}},{\widehat{\phi }}_{\text{oc}}\right)=\left(\pi -{\widehat{\theta }}_{\text{S}},\pi +{\widehat{\phi }}_{\text{S}}\right)$.

The oceanic frame of reference, ℛ_oc: (O, e_x, e_y, e_z), is such that e_z points towards the centre of the ocean basin, while e_x and e_y define two orthogonal axes in the plane normal to e_z. This frame of reference is associated with the spherical coordinates (r, θ, φ), where $r=\widehat{r}$ is the radial coordinate, θ the colatitude (θ = 0 at the centre of the ocean basin), and φ the longitude. The change of basis vectors (e_X, e_Y, e_Z) → (e_x, e_y, e_z) is expressed as a function of the coordinates of the continental centre as $$\left[\begin{array}{c}{e}_x\\ e_y\\ e_z\end{array}\right]=\left[\begin{array}{lll}\cos {\widehat{\phi }}_{\text{S}}\cos {\widehat{\theta }}_{\text{S}}\hfill & \sin {\widehat{\phi }}_{\text{S}}\cos {\widehat{\theta }}_{\text{S}}\hfill & -\sin {\widehat{\theta }}_{\text{S}}\hfill \\ sin{\widehat{\phi }}_{\text{S}}\hfill & -\cos {\widehat{\phi }}_{\text{S}}\hfill & 0\hfill \\ -\cos {\widehat{\phi }}_{\text{S}}\sin {\widehat{\theta }}_{\text{S}}\hfill & -\sin {\widehat{\phi }}_{\text{S}}\sin {\widehat{\theta }}_{\text{S}}\hfill & -\cos {\widehat{\theta }}_{\text{S}}\hfill \end{array}\right]\left[\begin{array}{c}{e}_X\\ e_Y\\ e_Z\end{array}\right].$$(1)

We note that the longitude of the continental centre does not affect the tidal response since the tidal forcing is periodic in longitude. Therefore, we set this longitude to ${\widehat{\phi }}_{\text{S}}=0$ in the following, similar to the example shown by Fig. 1, which yields $$\left[\begin{array}{c}{e}_x\\ e_y\\ e_z\end{array}\right]=\left[\begin{array}{lll}\cos {\widehat{\theta }}_{\text{S}}\hfill & 0\hfill & -\sin {\widehat{\theta }}_{\text{S}}\hfill \\ 0\hfill & -1\hfill & 0\hfill \\ -\sin {\widehat{\theta }}_{\text{S}}\hfill & 0\hfill & -\cos {\widehat{\theta }}_{\text{S}}\hfill \end{array}\right]\left[\begin{array}{c}{e}_X\\ e_Y\\ e_Z\end{array}\right].$$(2)

Finally, we introduce the set of unit vectors (e_r, e_θ, e_ϕ) associated with the spherical coordinates (r, θ, ϕ), respectively. The change of basis vectors (e_x, e_y, e_z) → (e_r, e_θ, e_ϕ) is given by the standard relation, $$\left[\begin{array}{c}{e}_r\\ e_{\theta }\\ e_{\phi }\end{array}\right]=\left[\begin{array}{lll}\sin \theta \cos \phi \hfill & \sin \theta \sin \phi \hfill & \cos \theta \hfill \\ \cos \theta \cos \phi \hfill & \cos \theta \cos \phi \hfill & -\sin \theta \hfill \\ -\sin \phi \hfill & \cos \phi \hfill & 0\hfill \end{array}\right]\left[\begin{array}{c}{e}_x\\ e_y\\ e_z\end{array}\right].$$(3)

Fig. 1 Geometry of the studied system and associated parameters. The unit vector $e{\prime }_z$, which indicates the position of the continental centre, is defined as $e{\prime }_z\equiv -e_z$, with ez pointing towards the centre of the ocean basin. In the diagram, the longitude of the continental centre is set to ${\widehat{\phi }}_{\text{S}}=0$.

2.2 Laplace’s tidal equations

The planet is subject to the effects of the tidal gravitational potential generated by the perturber (the Moon or the Sun in the Earth case), which is expressed, in the frame of reference ℛ, as $$U_{\text{T}}\equiv \frac{G{M}_{\text{s}}}{\left|R_{\text{p}}{e}_r-r_{\text{s}}\right|}-\frac{G{M}_{\text{s}}}{r_{\text{s}}^2}{R}_{\text{p}}\cos \widehat{\theta }-\frac{G{M}_{\text{s}}}{r_{\text{s}}},$$(4)

where the constant component and the component responsible for the Keplerian dynamics of the two-body system are removed (second and third terms in the right-hand member of Eq. (4)). In this equation, G is the universal gravitational constant, M_s the mass of the perturber, r_s its position vector, and r_s ≡ |r_s| the planet-perturber distance. We note that the symbol ≡ is used throughout the text to distinguish between definitions and equalities. The tidal force per unit mass exerted by the perturber on the planet is given by F ≡ ∇U_T. Following the formalism introduced in earlier works (e.g. Tyler 2011; Matsuyama 2014; Auclair-Desrotour et al. 2018, 2019a; Motoyama et al. 2020; Farhat et al. 2022a), we write down the momentum and mass conservation equations, respectively, as $${\partial }_{t}V+{\sigma }_{\text{R}}V+f\times V+g\nabla \left({\text{Γ}}_{\text{D}}\zeta -{\text{Γ}}_{\text{T}}{\zeta }_{\text{eq}}\right)=0,$$(5) $${\partial }_t\zeta +\nabla ·\left(HV\right)=0,$$(6)

with t designating the time, g the surface gravity at rest, σ_R the Rayleigh drag frequency used to describe the action of dissipa-tive mechanisms, V the horizontal velocity vector – defined from the horizontal displacement vector ξ as V ≡ ∂_tξ –, f the Coriolis parameter, ζ the vertical displacement of the ocean’s surface with respect to the oceanic floor, and ζ_eq ≡ U_T/g the equilibrium displacement corresponding to the equipotential surface induced by the tidal gravitational potential.

In the momentum equation, the notations Γ_D and Γ_T refer to non-trivial linear operators accounting for the effects of ocean loading, self-attraction, and deformation of the solid regions of the planet (e.g. Hendershott 1972). These operators are called the ‘solid deformation operators’ in the following. They encompass the complex coupling between the oceanic shell and the solid interior described by Poisson’s equation, and the momentum and rheological equations governing the tidal dynamics of the solid part. The two operators simplify to Γ_D = Γ_T = 1 if one neglects both the tidal deformation of the solid part (infinite-rigidity approximation) and the variation of self-attraction induced by the oceanic tidal response (Cowling approximation; see Cowling 1941; Unno et al. 1989). In the general case, they are determined by the rheological behaviour of the solid part in its response to gravitational and surface forcings.

Equations (5) and (6) are known as the Laplace’s tidal equations (hereafter, LTEs) in reference to Laplace’s masterpiece (Laplace 1798). We note that ζ_eq results from the coupled oceanic tidal response and tidal deformation of the solid part, which includes both gravitational and loading interactions, as discussed further. The Coriolis parameter is expressed as a function of the colatitude of the current point in ℛ, $$f\equiv 2\text{Ω\hspace{0.05em}cos}\widehat{\theta }{e}_r,$$(7)

where Ω > 0 is the planet’s spin angular velocity. The horizontal gradient operator and divergence of the horizontal velocity are expressed in ℛ_oc, and the associated basis vectors, (e_θ, e_ϕ), as $$\nabla \equiv R_{\text{p}}^{-1}\left[e_{\theta }{\partial }_{\theta }+e_{\phi }{\left(\sin \theta \right)}^{-1}{\partial }_{\phi }\right],$$(8) $$\nabla ·V={\left(R_{\text{p}}\sin \theta \right)}^{-1}\left[{\partial }_{\theta }\left(\sin \theta V_{\theta }\right)+{\partial }_{\phi }{V}_{\phi }\right].$$(9)

The Rayleigh drag term σ_R V in Eq. (5) accounts for the cumulated effects of dissipative mechanisms on tidal flows. The frequency σ_R is the inverse of the effective dissipation timescale associated with these dissipative mechanisms. It takes values around ~ 10⁻⁵ s⁻¹ for the present day Earth (e.g. Webb 1980; Farhat et al. 2022a)².

2.3 Nondimensional tidal equations

The nondimensional momentum and continuity equations are obtained by choosing as reference time and velocity scales the inertial time, t₀, and the typical velocity of long-wavelength surface gravity waves, V₀, which are defined, respectively, as $$\begin{array}{cc}{t}_0\equiv {\left(2\text{Ω}\right)}^{-1},& V_0\equiv \sqrt{gH}\end{array}.$$(10)

Introducing the normalised time $\tilde{t}$, Coriolis parameter $\tilde{f}$, horizontal gradient $\tilde{\nabla }$, the complex horizontal displacement vector $\tilde{\xi }$, and vertical displacements $\tilde{\zeta }$ and ${\tilde{\zeta }}_{\text{eq}}$, such that $$\begin{array}{lll}t=t_0\tilde{t},\hfill & f=t_0^{-1}\tilde{f}\hfill & \nabla =R_{\text{p}}^{-1}\tilde{\nabla },\hfill \\ \xi =\Re \left(H_{\text{p}}\tilde{\xi }\right),\hfill & \zeta =\Re \left(H\tilde{\zeta }\right),\hfill & {\zeta }_{\text{eq}}=\Re \left(H{\tilde{\zeta }}_{\text{eq}}\right),\hfill \end{array}$$(11)

with ℜ referring to the real part of a complex number, we end up with the nondimensional complex LTEs given by $$\left[{\partial }_{\tilde{t}\tilde{t}}+\left({\tilde{\sigma }}_{\text{R}}+\tilde{f}\text{}\times \right){\partial }_{\tilde{t}}\right]\tilde{\xi }+{\tilde{\sigma }}_{\text{G}}^2\tilde{\nabla }\tilde{\Xi }=0,$$(12) $$\tilde{\zeta }+\tilde{\nabla }·\tilde{\xi }=0,$$(13)

with $\widetilde{\text{Ξ}}\equiv {\text{Γ}}_{\text{D}}\tilde{\zeta }-{\text{Γ}}_{\text{T}}{\tilde{\zeta }}_{\text{eq}}$. In the above equations, tidal dynamics are controlled by two dimensionless parameters, $$\begin{array}{cc}{\tilde{\sigma }}_{\text{G}}\equiv \frac{V_0}{2\text{Ω}{R}_{\text{p}}}=\frac{\sqrt{gH}}{2\text{Ω}{R}_{\text{p}}},& {\tilde{\sigma }}_{\text{R}}\equiv \frac{{\sigma }_{\text{R}}}{2\text{Ω}}.\end{array}$$(14)

The first parameter, ${\tilde{\sigma }}_{\text{G}}$, can be considered as a normalised Rossby deformation length since it compares the typical propagation velocity of surface gravity waves (V₀) with the Earth’s rotation velocity. If ${\tilde{\sigma }}_{\text{G}}\ll 1$ (fast rotator regime), the inertial forces resulting from Coriolis acceleration predominate with respect to the restoring forces of surface gravity waves (pressure forces and gravity). Conversely, if ${\tilde{\sigma }}_{\text{G}}\gg 1$ (slow rotator regime), Coriolis terms are not strong enough to significantly deviate free surface gravity waves. As it describes the ratio of drag forces to Coriolis forces, the second dimensionless parameter, ${\tilde{\sigma }}_{\text{R}}$, may be regarded as an Ekman number. If ${\tilde{\sigma }}_{\text{R}}\ll 1$, the drag does not alter the tidal response much. Conversely, ${\tilde{\sigma }}_{\text{R}}\gg 1$ characterises a frictional (or viscous) regime where inertial effects are annihilated by the strong damping associated with the drag.

In addition to ${\tilde{\sigma }}_{\text{G}}$ and ${\tilde{\sigma }}_{\text{R}}$, the time-derivative operators in Eq. (12) introduce a third dimensionless parameter describing the ratio of tidal forces to Coriolis forces, $$\tilde{\sigma }\equiv \frac{\sigma }{2\text{Ω}}.$$(15)

The notation σ in the above equation designates the typical frequency of tidal flows, which will serve as the tidal frequency of the considered tidal force in the Fourier expansion of tidal quantities detailed further. If $\tilde{\sigma }\ll 1$ (sub-inertial regime), the acceleration term of the momentum equation can be neglected with respect to Coriolis terms. Conversely, if $\tilde{\sigma }\gg 1$ (super-inertial regime), the flow is strongly driven by the tidal forcing and it is thus hardly distorted by the planet’s rotation. We note that $\tilde{\sigma }$ is the inverse of the so-called spin parameter that is commonly used to characterise the oscillations of rotating fluids in planetary and stellar hydrodynamics (e.g. Lee & Saio 1997). The dimensionless parameters that control the planet’s tidal response are summarised in Table 1.

2.4 Helmholtz decomposition

Proudman (1920) demonstrated that Helmholtz’s theorem (e.g. Arfken & Weber 2005) can be used to decompose the horizontal displacement vector field into curl-free and divergence-free vector fields, $$\tilde{\xi }=\tilde{\nabla }\text{Φ\hspace{0.05em}\hspace{0.17em}}+\tilde{\nabla }\text{Ψ}\times e_r.$$(16)

The curl-free $\left(\tilde{\nabla }\times \left(\tilde{\nabla }\text{Φ}\right)=0\right)$ and divergence-free $\left(\tilde{\nabla }\cdot \left(\tilde{\nabla }\text{Ψ}\times e_r\right)=0\right)$ components of $\tilde{\xi }$ are defined from the divergent displacement potential Φ and the rotational displacement stream-function Ψ, respectively (Webb 1980; Tyler 2011). We emphasise that the Helmholtz decomposition is not unique for bounded domains such as the considered ocean basin. This results from the fact that additional physical constraints on the boundary condition are necessary to define Φ and Ψ (e.g. Fox-Kemper et al. 2003). However choosing boundary conditions for the two functions is not straightforward given that the two components of Eq. (16) cannot be disentangled in the total flux. For instance, the impermeability condition at the coastline is formulated as $\tilde{\xi }\cdot n=0$, where n designates the outward pointing unit vector defining the normal to the coast. This theoretically requires to find another boundary condition and to solve it for $\tilde{\nabla }\text{Φ}$ and $\tilde{\nabla }\text{Ψ}$ together with the first condition, which may lead to significant mathematical complications.

To circumvent this difficulty, it is convenient to adopt a commonly used trick (e.g. Webb 1980, 1982; Gent & McWilliams 1983; Watterson 2001; Han & Huang 2020; Farhat et al. 2022a), which consists in applying the impermeability condition to both components of Eq. (16), $$\begin{array}{cc}n·\tilde{\nabla }\text{Φ=0,}& n·\left(\tilde{\nabla }\text{Ψ}\times e_r\right)=0.\end{array}$$(17)

The first condition of Eq. (17) means that the gradient of Φ is zero in the direction normal to the coastline, which is equivalent to Neumann condition (specified value of the derivative of the solution applied at the boundary of the domain; e.g. Morse & Feshbach 1953, Sect. 6.1). The second condition may be reformulated as $\left(e_r\times n\right)\cdot \tilde{\nabla }\Psi =0$, implying that the streamfunction is a constant along the coastline. This corresponds to a Dirich-let condition (specified value of the solution itself; Morse & Feshbach 1953, Sect. 6.1). Since both Φ and Ψ are defined to a constant, the value of ψ at the coastline is set to zero, which simplifies the second condition of Eq. (17) to Ψ = 0. Interestingly, this condition induces the orthogonality of the curl-free and divergence-free components of the horizontal displacement (e.g. Farhat et al. 2022a, Appendix E), formulated as $${\displaystyle {\int }_{𝒪}\overline{\tilde{\nabla }\text{Φ}}·\left(\tilde{\nabla }\text{Ψ}\times e_r\right)\text{d}S=0,}$$(18)

with dS = sin θdθdϕ being an infinitesimal surface element of the unit sphere, 𝒪 the domain occupied by the ocean basin, and $\overline{\tilde{\nabla }\text{Φ}}$ the complex conjugate of $\tilde{\nabla }\text{Φ}$. This property leads to appreciable simplifications in the LTEs, as discussed further.

The divergent potential function and the streamfunction are defined on a compact connected domain, the ocean basin (𝒪), and they satisfy either Neumann or Dirichlet conditions at its boundary, ∂𝒪. As a consequence, the two functions can be expanded in terms of the complete sets of orthogonal eigenfunc-tions {Φ_k}_1≤k≤∞ or {Φ_k}_1≤k≤∞ that are the solutions of the wave equations given by $$\begin{array}{cc}\left({\tilde{\nabla }}^2+\lambda \right)\text{Φ=0\hspace{0.17em}on\hspace{0.17em}}𝒪,& n·\end{array}\tilde{\nabla }\text{Φ\hspace{0.17em}=\hspace{0.17em}}0\text{\hspace{0.17em}at\hspace{0.17em}}\partial 𝒪,$$(19) $$\begin{array}{cc}\left({\tilde{\nabla }}^2+v\right)\text{Ψ\hspace{0.17em}=\hspace{0.17em}}0\text{\hspace{0.17em}on\hspace{0.17em}𝒪,}& \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}Ψ}=0\text{\hspace{0.17em}at\hspace{0.17em}}\partial 𝒪,\end{array}$$(20)

where ${\tilde{\nabla }}^2$ designates the normalised Laplacian operator defined, for any function f, as $${\tilde{\nabla }}^{2}f\equiv {\left(\sin \theta \right)}^{-2}\left[\sin \theta {\partial }_{\theta }\left(\sin \theta {\partial }_{\theta }f\right)+{\partial }_{\phi \phi }f\right].$$(21)

The solutions of the colatitudinal component of Eqs. (19) and (20) are the associated Legendre functions of the first kind (ALFs) of real degrees, which are detailed in Appendix B. The corresponding eigenfunctions Φ_k are associated with the eigenvalues λ_k ∈ ℝ⁺, and the eigenfunctions ψ_k with the eigenvalues v_k ∈ ℝ⁺, the relationship between the real degrees of the ALFs and the eigenvalues being detailed in Appendix C (see Table C.2). These eigenfunctions are obtained by multiplying the ALFs to the solutions of the longitudinal component of Eqs. (19) and (20), namely exp (±_imϕ), with i being the imaginary number. The outcome is a set of functions known as spherical cap harmonics (Haines 1985; Hwang & Chen 1997; Thébault et al. 2004, 2006), and denoted by ‘SCHs’ in the following (see Eq. (C.19)). Thus, the set {Φ_k}_1≤k≤∞ is refered to as the ‘SCHNs’ (Neumann condition), and the set {Φ_k}_1≤k≤∞ as the ‘SCHDs’ (Dirichlet condition). For example, Fig. 2 shows the first basis functions of the two sets for a supercontinent of angular radius θ_c = 10° (i.e. θ₀ = 170°).

One shall notice that the SCHs and their associated eigenvalues both depend on the boundary condition applied at the coastline, meaning that the sets {Φ_k}_1≤k≤∞ and {Φ_k}_1≤k≤∞ are necessarily different by construction. In each set, the basis functions are orthogonal to each other and normalised so that $${\displaystyle {\int }_{𝒪}\overline{{\text{Φ}}_k}{\text{Φ}}_j\text{d}S={\displaystyle {\int }_{𝒪}\overline{{\text{Ψ}}_k}{\text{Ψ}}_j\text{d}S={\delta }_{k,j},}}$$(22)

where δ_k,j designates the Kronecker delta function, such that δ_k,j = 1 for k = j and δ_k,j = 0 otherwise. However, two functions belonging to different sets are not orthogonal in the general case, $${\displaystyle {\int }_{𝒪}\overline{{\text{Φ}}_k}}{\text{Ψ}}_j\text{d}S\ne 0,$$(23)

since the sets of basis functions {Φ_k}_1≤k≤∞ and {Φ_k}_1≤k≤∞ are not the same.

Introducing the time-dependent coefficients p_k and p_−k, we write down the divergent potential function and the streamfunction as $$\text{Φ}\left(\theta ,\phi ,\tilde{t}\right)={\displaystyle \sum _{k=1}^{\infty }{p}_k\left(\tilde{t}\right){\text{Φ}}_k\left(\theta ,\phi \right)},$$(24) $$\text{Ψ}\left(\theta ,\phi ,\tilde{t}\right)={\displaystyle \sum _{k=1}^{\infty }{p}_{-k}\left(\tilde{t}\right){\text{Ψ}}_k\left(\theta ,\phi \right)}.$$(25)

By substituting Eq. (24) into Eq. (13), the continuity equation becomes $$\tilde{\zeta }+{\tilde{\nabla }}^2\text{Φ\hspace{0.17em}=\hspace{0.17em}0,}$$(26)

which implies that $$\tilde{\zeta }={\displaystyle \sum _{k=1}^{\infty }{\lambda }_k{p}_k{\text{Φ}}_k}.$$(27)

However, such a simple expression cannot be obtained for the equilibrium displacement (${\tilde{\zeta }}_{\text{eq}}$) in the general case since it is naturally expanded in series of spherical harmonics (hereafter, SPHs) in the coordinate system associated with the planet’s spin, which do not satisfy Dirichlet or Neumann conditions at coastlines. The form of the expression given by Eq. (27) can be obtained for ${\tilde{\zeta }}_{\text{eq}}$ only in the particular case of the hemispheric ocean (θ₀ = 90°), where the set of eigenfunctions {Φ_k}_1≤k≤∞ actually corresponds to a subset of the SPHs for the hemispherical domain, as shown by Farhat et al. (2022a).

Fig. 2 Eigenfunctions describing the oceanic tidal response for a supercontinent of angular radius θc = 10° (i.e. θ0 = 170°), and associated tidal flows. Top: set {Φk} (SCHNs). Bottom: set {Ψk} (SCHDs). The orange disk designates the continent. The eigenfunctions are computed from the expression given by Eq. (C.19), and their real parts are plotted for ln such that n = 0,…, 3 (vertical axis) and −n ≤ m ≤ n (horizontal axis). Bright or dark colours designate positive or negative values of the eigenfunctions, respectively. Streamlines indicate the tidal flows corresponding to $\tilde{\nabla }{\text{Φ}}_k$ for the set {Φk} and to $\tilde{\nabla }{\text{Ψ}}_k\times e_r$ for the set {Ψk}.

2.5 Temporal equations

The method used to establish the temporal differential equations for the coefficients p_k and p_−k introduced in Eqs. (24) and (25) closely follows that used for the hemispheric ocean configuration (e.g. Webb 1980, 1982; Farhat et al. 2022a). As a first step, we substitute Φ and Ψ by their series expansions in Eq. (16) and in the momentum equation given by Eq. (12). As a second step, we do the dot product of the equation thus obtained by $\overline{\tilde{\nabla }{\text{Φ}}_k}$, and we integrate it over the domain of the ocean basin. It follows $$\begin{array}{l}{\displaystyle \sum _{j=1}^{\infty }\{\left[{\partial }_{\tilde{t}\tilde{t}}+{\tilde{\sigma }}_{\text{R}}{\partial }_{\tilde{t}}\right]p_j{\displaystyle {\int }_{𝒪}\overline{\tilde{\nabla }{\text{Φ}}_k}·}\tilde{\nabla }{\text{Φ}}_j\text{d}S}\hfill \\ +\left({\partial }_{\tilde{t}\tilde{t}}+{\tilde{\sigma }}_{\text{R}}{\partial }_{\tilde{t}}\right)p_{-j}{\displaystyle {\int }_{𝒪}\overline{\tilde{\nabla }{\text{Φ}}_k}·}\left(\tilde{\nabla }{\text{Ψ}}_j\times e_r\right)\text{d}S\hfill \\ +{\partial }_{\tilde{t}}{p}_{-j}{\displaystyle {\int }_{𝒪}\overline{\tilde{\nabla }{\text{Φ}}_k}·}\left[\tilde{f}\times \left(\tilde{\nabla }{\text{Ψ}}_j\times e_r\right)\right]\text{d}S\hfill \\ +{\partial }_{\tilde{t}}{p}_j{\displaystyle {\int }_{𝒪}\overline{\tilde{\nabla }{\text{Φ}}_k}·}\left(\tilde{f}\times \tilde{\nabla }{\text{Φ}}_j\right)\text{d}S\}+{\tilde{\sigma }}_{\text{G}}^2{\displaystyle {\int }_{𝒪}\overline{\tilde{\nabla }{\text{Φ}}_k}·\tilde{\nabla }\widetilde{\text{Ξ}}\text{d}S=0.}\hfill \end{array}$$(28)

The dot product of two eigenfunctions that appears in the first term of Eq. (28) is computed by invoking, successively, Green’s first identity (e.g. Strauss 2007, Chapter 7), $${\displaystyle {\int }_{𝒪}\overline{\tilde{\nabla }{\text{Φ}}_k}·}\tilde{\nabla }{\text{Φ}}_j\text{d}S={\displaystyle {\int }_{\partial 𝒪}{\text{Φ}}_j\left(\overline{\tilde{\nabla }{\text{Φ}}_k}·n\right)\text{d}\ell }-{\displaystyle {\int }_{𝒪}{\text{Φ}}_j\overline{{\tilde{\nabla }}^2{\text{Φ}}_k}\text{d}S,}$$(29)

the impermeability boundary condition on φ_k given by Eq. (17), the fact that the φ_k are eigenfunctions of the Laplacian operator as described by Eq. (19), and the orthogonality property given by Eq. (22). This yields $${\displaystyle {\int }_{𝒪}\overline{\tilde{\nabla }{\text{Φ}}_k}·}\tilde{\nabla }{\text{Φ}}_j\text{d}S={\lambda }_k{\delta }_{k,j},$$(30)

and for the forcing term, $${\displaystyle {\int }_{𝒪}\overline{\tilde{\nabla }{\text{Φ}}_k}·}\tilde{\nabla }\widetilde{\text{Ξ}}\text{d}S={\lambda }_k{\displaystyle {\int }_{𝒪}\overline{{\text{Φ}}_k}\widetilde{\text{Ξ}}}\text{d}S.$$(31)

Since the tidal gravitational potential is expressed in terms of the SPHs associated with the coordinates ($\widehat{\theta }$, $\widehat{\phi }$), the calculation of the above integral requires computing the transition matrix between these SPHs and the oceanic eigenfunctions. This calculation is achieved in two steps. First, one computes the rotated SPHs associated with the change of coordinates $\left(\widehat{\theta },\widehat{\phi }\right)\to \left(\theta ,\phi \right)$, as detailed in Appendices D and E. Second, the transition matrices between the non-rotated SPHs and the SCHs are evaluated following the method described in Appendix F.

The integral of the second term in the left-hand member of Eq. (28) actually vanishes owing to the orthogonality property of the curl-free and divergence-free components of the horizontal displacement (Eq. (18)), $${\displaystyle {\int }_{𝒪}\overline{\tilde{\nabla }{\text{Φ}}_k}·\left(\tilde{\nabla }{\text{Ψ}}_j\times e_r\right)\text{d}S=0.}$$(32)

Finally, the integrals of the third and fourth terms of Eq. (28) are simplified with the help of vectorial identities, yielding $$\begin{array}{l}{\displaystyle {\int }_{𝒪}\overline{\tilde{\nabla }{\text{Φ}}_k}·\left[\tilde{f}\times \left(\tilde{\nabla }{\text{Ψ}}_j\times e_r\right)\right]\text{d}S=}{\displaystyle {\int }_{𝒪}\left(\tilde{f}·e_r\right)}\overline{\tilde{\nabla }{\text{Φ}}_k}·\tilde{\nabla }{\text{Ψ}}_j\text{d}S,\hfill \\ {\displaystyle {\int }_{𝒪}\overline{\tilde{\nabla }{\text{Φ}}_k}·}\left(\tilde{f}\times \tilde{\nabla }{\text{Φ}}_j\right)\text{d}S=-{\displaystyle {\int }_{𝒪}\tilde{f}·\left(\overline{\tilde{\nabla }{\text{Φ}}_k}\times \tilde{\nabla }{\text{Φ}}_j\right)\text{d}S.}\hfill \end{array}$$(33)

The above simplifications lead to the first set of temporal equations for the coefficients p_j and p−j, $$\begin{array}{l}{\displaystyle \sum _{j=1}^{\infty }\{\left[{\partial }_{\tilde{t}\tilde{t}}+{\tilde{\sigma }}_{\text{R}}{\partial }_{\tilde{t}}\right]p_j{\lambda }_k{\delta }_{k,j}+{\partial }_{\tilde{t}}{p}_{-j}}{\displaystyle {\int }_{𝒪}\left(\tilde{f}·e_r\right)\overline{\tilde{\nabla }{\text{Φ}}_k}·\tilde{\nabla }{\text{Ψ}}_j\text{d}S}\hfill \\ -{\partial }_{\tilde{t}}{p}_j{\displaystyle {\int }_{𝒪}\tilde{f}·\left(\overline{\tilde{\nabla }{\text{Φ}}_k}\times \tilde{\nabla }{\text{Φ}}_j\right)\text{d}S}\}+{\tilde{\sigma }}_{\text{G}}^2{\lambda }_k{\displaystyle {\int }_{𝒪}\overline{{\text{Φ}}_k}\widetilde{\text{Ξ}}}\text{d}S=0.\hfill \end{array}$$(34)

The second set of equations is obtained by repeating the same steps with $\overline{\tilde{\nabla }{\text{Ψ}}_k}\times e_r$ instead of $\overline{\tilde{\nabla }{\text{Φ}}_k}$ in the dot product operation, $$\begin{array}{l}{\displaystyle \sum _{j=1}^{\infty }\{\left({\partial }_{\tilde{t}\tilde{t}}+{\tilde{\sigma }}_{\text{R}}{\partial }_{\tilde{t}}\right)p_j}{\displaystyle {\int }_{𝒪}\left(\overline{\tilde{\nabla }{\text{Ψ}}_k}\times e_r\right)·\tilde{\nabla }{\text{Φ}}_j\text{d}S}\hfill \\ +{\partial }_{\tilde{t}}{p}_j{\displaystyle {\int }_{𝒪}\left(\overline{\tilde{\nabla }{\text{Ψ}}_k}\times e_r\right)·\left(\tilde{f}\times \tilde{\nabla }{\text{Φ}}_j\right)\text{d}S}\hfill \\ +{\partial }_{\tilde{t}}{p}_{-j}{\displaystyle {\int }_{𝒪}\left(\overline{\tilde{\nabla }{\text{Ψ}}_k}\times e_r\right)·[\tilde{f}\times \left(\tilde{\nabla }{\text{Ψ}}_j\times e_r\right)]}\text{\hspace{0.17em}d}S\hfill \\ +\left({\partial }_{\tilde{t}\tilde{t}}+{\tilde{\sigma }}_{\text{R}}\right)p_{-j}{\displaystyle {\int }_{𝒪}\left(\overline{\tilde{\nabla }{\text{Ψ}}_k}\times e_r\right)·\left(\tilde{\nabla }{\text{Ψ}}_j\times e_r\right)}\text{d}S\}\hfill \\ -{\tilde{\sigma }}_{\text{G}}^2{\displaystyle {\int }_{𝒪}{e}_r·\left(\overline{\tilde{\nabla }{\text{Ψ}}_k}\times \tilde{\nabla }\widetilde{\text{Ξ}}\right)}\text{d}S=0,\hfill \end{array}$$(35)

which, using the orthogonality properties of the eigenfunctions and the relations $${\displaystyle {\int }_{𝒪}\overline{\tilde{\nabla }{\text{Ψ}}_k}·}\tilde{\nabla }{\text{Ψ}}_j\text{d}S=v_k{\delta }_{k,j},$$(36)

simplifies to $$\begin{array}{l}{\displaystyle \sum _{j=1}^{\infty }\{\left({\partial }_{\tilde{t}\tilde{t}}+{\tilde{\sigma }}_{\text{R}}\right)p_{-j}{v}_k{\delta }_{k,j}-{\partial }_{\tilde{t}}{p}_j{\displaystyle {\int }_{𝒪}\left(\tilde{f}·e_r\right)\left(\overline{\tilde{\nabla }{\text{Ψ}}_k}·\tilde{\nabla }{\text{Φ}}_j\right)\text{d}S}}\hfill \\ -{\partial }_{\tilde{t}}{p}_j{\displaystyle {\int }_{𝒪}\tilde{f}·}\left(\overline{\tilde{\nabla }{\text{Ψ}}_k}\times \tilde{\nabla }{\text{Ψ}}_j\right)\text{d}S\}-{\tilde{\sigma }}_{\text{G}}^2{\displaystyle {\int }_{𝒪}{e}_r·\left(\overline{\tilde{\nabla }{\text{Ψ}}_k}\times \tilde{\nabla }\widetilde{\text{Ξ}}\right)\text{d}S-0.}\hfill \end{array}$$(37)

Thus, following the formalism used in earlier studies (Webb 1980, 1982; Farhat et al. 2022a), Eqs. (34) and (37) form an infinite linear system in the coefficients p_k and p_−k, written as $$\left({\partial }_{\tilde{t}\tilde{t}}+{\tilde{\sigma }}_{\text{R}}{\partial }_{\tilde{t}}\right)p_k+{\lambda }_k^{-1}{\displaystyle \sum _{\underset{j\ne 0}{j=-\infty }}^{+\infty }{\beta }_{k,j}{\partial }_{\tilde{t}}{p}_j+{\tilde{\sigma }}_{\text{G}}^2{\displaystyle {\int }_O\overline{{\text{Φ}}_k}\widetilde{\text{Ξ}}\text{d}S=0}},$$(38) $$\begin{array}{c}\left({\partial }_{\tilde{t}\tilde{t}}+{\tilde{\sigma }}_{\text{R}}{\partial }_{\tilde{t}}\right)p_{-k}+v_k^{-1}{\displaystyle \sum _{\underset{j\ne 0}{j=-\infty }}^{+\infty }{\beta }_{-k,j}{\partial }_{\tilde{t}}{p}_j}\\ -\frac{{\tilde{\sigma }}_{\text{G}}^2}{v_k}{\displaystyle {\int }_O{e}_r\cdot \left(\overline{\tilde{\nabla }{\text{Ψ}}_k}\times \tilde{\nabla }\widetilde{\text{Ξ}}\right)\text{d}S=0.}\end{array}$$(39)

In the above equations, the symbols β_k,j designate the so-called ‘gyroscopic coefficients’ (e.g. Proudman 1920; Longuet-Higgins & Pond 1970; Webb 1980), $${\beta }_{k,j}\equiv -{\displaystyle {\int }_{𝒪}\cos \widehat{\theta }{e}_r\cdot \left(\overline{\nabla {\text{Φ}}_k}\times \nabla {\text{Φ}}_j\right)\text{d}S,}$$(40) $${\beta }_{k,-j}\equiv {\displaystyle {\int }_{𝒪}\cos \widehat{\theta }\cdot \left(\overline{\nabla {\text{Φ}}_k}\cdot \nabla {\text{Ψ}}_j\right)\text{d}S,}$$(41) $${\beta }_{-k,j}\equiv -{\displaystyle {\int }_{𝒪}\cos \widehat{\theta }\cdot \left(\overline{\nabla {\text{Ψ}}_k}\cdot \nabla {\text{Φ}}_j\right)\text{d}S,}$$(42) $${\beta }_{-k,-j}\equiv {\displaystyle {\int }_{𝒪}\cos \widehat{\theta }{e}_r\cdot \left(\overline{\nabla {\text{Ψ}}_k}\times \nabla {\text{Ψ}}_j\right)\text{d}S}.$$(43)

These coefficients account for the coupling effect of the Coriolis terms and the geometry of the ocean basin, which affects how the gradients of the eigenfunctions overlap. Their behaviour and properties are examined in Appendix G. The expression of the forcing term in Eq. (39) is given here for generality. This term is actually zero, as discussed in the next section, meaning that Eq. (39) does not depend on $\widetilde{\text{Ξ}}$.

2.6 Ocean loading and self-attraction variation

At this stage, we still have to express $\widetilde{\text{Ξ}}$ as a function of p_j and p_−j. To do so, we consider the fact that the gravitational tidal force acts on the planet as a periodic perturbation oscillating in time and longitude. Therefore, the tidal gravitational potential given by Eq. (4) can be expanded in Fourier series of time and series of SPHs. In the general case, U_T is expressed as $$U_{\text{T}}{V}_0^2\Re \left\{{\displaystyle \sum _{\sigma }{\displaystyle \sum _{l=2}^{+\infty }{\displaystyle \sum _{m=-l}^l{\tilde{U}}_{\text{T};l}^{m,\sigma }{\widehat{Y}}_l^m\left(\widehat{\theta },\widehat{\phi }\right){\text{e}}^{i\sigma t}}}}\right\},$$(44)

where σ is the tidal frequency, V₀ the reference velocity introduced in Eq. (10), ${\widehat{Y}}_l^m$ the complex SPH of degree l and order m associated with the coordinate system ($\widehat{\theta }$, $\widehat{\phi }$), defined in Eq. (C.9), and ${\tilde{U}}_{\text{T};l}^{m,\sigma }\in ℂ$ the associated frequency-dependent normalised component of the tidal gravitational potential. In the following, these notations are shortened for convenience. Similarly as the eigenfunctions of the ocean basin, Φ_k and Ψ_k, the SPHs are simply denoted by ${\widehat{Y}}_k={\widehat{Y}}_{l_k}^{m_k}$, and the corresponding coefficients by ${\tilde{U}}_{\text{T};k}^{m,\sigma }=U_{\text{T;}{l}_k}^{m_k,\sigma }$, the index k referring to an element of the set of SPHs, ${\left\{{\widehat{Y}}_k\right\}}_{\text{l}\le k\le \infty }$. Moreover, all the indices used from now on (k, j, n, q) are supposed to run from one to infinity when bounds are not specified. It is noteworthy that the Fourier series expansion in Eq. (44) can always be written in term of positive tidal frequencies as long as m runs from −l to l. Therefore, we assume that σ ≥ 0, and that the frequencies are all different from each other (no resonance).

Since the components associated with two different tidal frequencies are not correlated in the linear tidal theory, the resulting tidal responses can be treated separately. We thus consider the contribution of the component associated with a given tidal frequency, ${\tilde{U}}_{\text{T}}^{\sigma }$, and the associated forcing term, ${\widetilde{\text{Ξ}}}^{\sigma }$, expressed as $$\begin{array}{cc}{\tilde{U}}_{\text{T}}^{\sigma }={\displaystyle \sum _j{\tilde{U}}_{\text{T;}j}^{\sigma }{\widehat{Y}}_j\left(\widehat{\theta },\widehat{\phi }\right),}& \widetilde{\text{Ξ}}\sigma ={\displaystyle \sum _j{\widetilde{\text{Ξ}}}_{\text{T;}j}^{\sigma }{\widehat{Y}}_j\left(\widehat{\theta },\widehat{\phi }\right).}\end{array}$$(45)

As highlighted in earlier studies (e.g. Matsuyama 2014; Matsuyama et al. 2018; Auclair-Desrotour et al. 2019a), the combined contribution of ocean loading, self-attraction variation, and deformation of solid regions may be formulated as $${\widetilde{\text{Ξ}}}_j^{\sigma }={\gamma }_{\text{D;}{l}_j}^{\sigma }\left({\displaystyle \sum _n\langle {\widehat{Y}}_j,{\text{Φ}}_n\rangle {\lambda }_n{p}_n}\right)-{\gamma }_{\text{T};l_j}^{\sigma }{\tilde{U}}_{\text{T};l_j}^{\sigma },$$(46)

where we have made use of the scalar product defined by Eq. (C.12) and introduced the solid deformation factors, ${\gamma }_{\text{D};l}^{\sigma }$ and ${\gamma }_{\text{T};l}^{\sigma }$ (see e.g. Hendershott 1972), $${\gamma }_{\text{D};l}^{\sigma }\equiv 1+k_l^{\sigma }-h_l^{\sigma },$$(47) $${\gamma }_{\text{D};l}^{\sigma }\equiv 1-\left(1+k_{\text{L};l}^{\sigma }-h_{\text{L};l}^{\sigma }\right)\frac{3}{\left(2l+1\right)}\frac{{\rho }_{\text{w}}}{\overline{\rho }}.$$(48)

The factors ${\gamma }_{\text{D};l}^{\sigma }$ and ${\gamma }_{\text{T};l}^{\sigma }$ are independent of the order m since they are expressed as functions of the Love numbers of the solid part, which is commonly assumed to be spherically symmetric in the tidal theory. The parameters $k_l^{\sigma }$ and $h_l^{\sigma }$ in the first expression are known as the tidal gravitational and displacement Love numbers, respectively. They describe the response of the solid body to a gravitational tidal forcing, which takes the form of a self-attraction variation and a surface displacement. By analogy, the loading Love numbers $k_{\text{L};l}^{\sigma }$ and $h_{\text{L};l}^{\sigma }$ account for the gravitational and mechanical responses of solid regions to the cumulated gravitational and pressure forces generated by the oceanic mass redistribution, respectively. The second solid deformation factor, ${\gamma }_{\text{T};l}^{\sigma }$, also depends on the ratio of seawater density to the mean density of the solid regions, $\overline{\rho }$. We remark that ${\gamma }_{D;l}^{\sigma }$ and ${\gamma }_{\text{T};l}^{\sigma }$ are sometimes considered as real factors (e.g. Matsuyama 2014; Motoyama et al. 2020), which corresponds to an adiabatic elastic response of the solid part, where dissipative mechanisms are ignored. However, the correspondence principle established by Biot (1954) makes it possible to treat dissipative anelastic cases similarly as long as the anelasticity is linear. As a consequence, the expressions given by Eqs. (47) and (48) can be extended to any rheological model including dissipative processes. In that case, the solid Love numbers are complex transfer functions describing the visco-elastic response of solid regions when subjected to a harmonic tidal force (Remus et al. 2012; Auclair-Desrotour et al. 2019a; Farhat et al. 2022b).

Using the expression given by Eq. (46) and proceeding to a change of basis functions between the oceanic eigen functions and the SPHs associated with the coordinate system ($\widehat{\theta }$, $\widehat{\phi }$), we expand the integrals depending on the forcing in Eqs. (38) and (39) in series of p_j and components of the tidal gravitational potential, $$\begin{array}{l}{\displaystyle {\int }_{𝒪}\overline{{\text{Φ}}_k}\widetilde{\text{Ξ}}\text{d}S={\displaystyle \sum _n\left[{\displaystyle \sum _j\langle {\Phi }_k,{\widehat{Y}}_j\rangle {\gamma }_{\text{D};l_j}^{\sigma }\langle {\widehat{Y}}_j,{\text{Φ}}_n\rangle }\right]{\lambda }_n{p}_n}}\hfill \\ -{\displaystyle \sum _j\langle {\text{Φ}}_k,{\widehat{Y}}_j\rangle {\gamma }_{\text{T};l_j}^{\sigma }{\tilde{U}}_{\text{T};j}^{\sigma },}\hfill \end{array}$$(49) $$\begin{array}{l}{\displaystyle {\int }_{𝒪}{e}_r\cdot }\left(\overline{\tilde{\nabla }{\text{Ψ}}_k}\times \tilde{\nabla }\widetilde{\text{Ξ}}\right)\text{d}S={\displaystyle \sum _n\left[{\displaystyle \sum _{j,q}{\upsilon }_{k,q}\langle Y_q,{\widehat{Y}}_j\rangle {\gamma }_{\text{D};l_j}^{\sigma }\langle {\widehat{Y}}_j,{\text{Φ}}_n\rangle }\right]{\lambda }_n{p}_n}\hfill \\ -{\displaystyle \sum _{j,q}{\upsilon }_{k,q}\langle Y_q,{\widehat{Y}}_j\rangle {\gamma }_{\text{T};l_j}^{\sigma }{\tilde{U}}_{\text{T};j}^{\sigma }},\hfill \end{array}$$(50)

where υ_k,q designates the coefficients defined as $${\upsilon }_{k,q}\equiv {\displaystyle {\int }_{𝒪}{e}_r\cdot \left({\overline{\tilde{\nabla }\text{Ψ}}}_k\times \tilde{\nabla }{Y}_q\right)}\text{\hspace{0.17em}d}S.$$(51)

These coefficients are all zero as shown in Appendix G, which implies that $${\displaystyle {\int }_{𝒪}{e}_r\cdot \left({\overline{\tilde{\nabla }\text{Ψ}}}_k\times \tilde{\nabla }\widetilde{\text{Ξ}}\right)}\text{\hspace{0.17em}d}S=0.$$(52)

We remark that the quadrupolar component of the tidal gravitational potential (the component associated with the ${\widehat{Y}}_2^2$ SPH) is far greater than the components of higher degrees if the size of the planet is small compared with the planet-perturber distance. These components can thus be ignored in standard two-body systems. The series in the formulation of the tidal gravitational given by Eq. (45) then reduces to one term only. This removes the summation over j in the second term of the right-hand member of Eq. (49).

2.7 The tidal solution

As a last step, the tidal equations given by Eqs. (38) and (39) are written down in the Fourier domain. To compute the solution numerically, all the sets of basis functions are truncated, meaning that the infinite spaces of functions they describe are approximated by finite spaces of functions. As the tidal gravitational potential varies over planetary length scales, the tidal response is essentially described by the eigenmodes of largest wavelengths. As a consequence, the truncation does not alter much the solution as long as the number of eigenmodes of the set is sufficiently large. The numbers of basis functions for the sets {Φ_j}, {Ψ_j}, and $\left\{{\widehat{Y}}_j\right\}$ are denoted by M, N, and K, respectively, while the vectors describing the tidal response and the tidal gravitational potential in these sets are expressed as $$\Phi ={\left[p_1,\dots ,p_j,\dots ,p_M\right]}^{\text{T}},$$(53) $$\Psi ={\left[p_{-1},\dots ,p_{-j},\dots ,p_{-N}\right]}^{\text{T}},$$(54) $${\tilde{U}}_{\text{T}}={\left[{\tilde{U}}_{\text{T};1}^{\sigma },\dots ,{\tilde{U}}_{\text{T};j}^{\sigma },\dots ,{\tilde{U}}_{\text{T};K}^{\sigma }\right]}^{\text{T}},$$(55)

with ^T designating the transpose of a matrix. Converted into an algebraic form, the forcing term given by Eq. (49) is written as $${\displaystyle {\int }_{𝒪}\overline{{\text{Φ}}_k}\widetilde{\text{Ξ}}\text{d}S}=E_{M,k}^T\left(A_{\text{Φ}}\Phi +R_{\text{Φ}}{\tilde{U}}_{\text{T}}\right),$$(56)

where we have introduced the unit vectors E_M,k defined as $$E_{M,k}\equiv {\left[{\delta }_{1,k},\dots ,{\delta }_{j,k},\dots ,{\delta }_{M,k}\right]}^{\text{T}},$$(57)

and the matrices Α_Φ, and R_Φ, $$\begin{array}{cc}{A}_{\text{Φ}}={\overline{P_{\widehat{Y},\text{Φ}}}}^{\text{T}}{\Gamma }_{\text{D}}{P}_{\widehat{Y},\text{Φ}}\Lambda ,& R_{\text{Φ}}=-\end{array}{\overline{P_{\widehat{Y},\text{Φ}}}}^{\text{T}}{\Gamma }_{\text{T}}.$$(58)

In the above expressions, $P_{\widehat{Y},\text{Φ}}$ designates the transition matrix from the SPHs ${\widehat{Y}}_k$ to the eigenfunctions Φ_k. The other matrices are defined as $${\Gamma }_{\text{D}}\equiv \left[\begin{array}{ccccc}{\gamma }_{\text{D};l_1}^{\sigma }& 0& \dots & \dots & 0\\ 0& \ddots & \ddots & & ⋮\\ ⋮& \ddots & {\gamma }_{\text{D};l_j}^{\sigma }& \ddots & ⋮\\ ⋮& & \ddots & \ddots & 0\\ 0& \dots & \dots & 0& {\gamma }_{\text{D};l_K}^{\sigma }\end{array}\right],$$(59) $${\Gamma }_{\text{T}}\equiv \left[\begin{array}{ccccc}{\gamma }_{\text{T};l_1}^{\sigma }& 0& \dots & \dots & 0\\ 0& \ddots & \ddots & & ⋮\\ ⋮& \ddots & {\gamma }_{\text{T};l_j}^{\sigma }& \ddots & ⋮\\ ⋮& & \ddots & \ddots & 0\\ 0& \dots & \dots & 0& {\gamma }_{\text{T};l_K}^{\sigma }\end{array}\right],$$(60) $$\Lambda \equiv \left[\begin{array}{ccccc}{\lambda }_1& 0& \dots & \dots & 0\\ 0& \ddots & \ddots & & ⋮\\ ⋮& \ddots & {\lambda }_j& \ddots & ⋮\\ ⋮& & \ddots & \ddots & 0\\ 0& \dots & \dots & 0& {\lambda }_M\end{array}\right],$$(61) $$N\equiv \left[\begin{array}{ccccc}{v}_1& 0& \dots & \dots & 0\\ 0& \ddots & \ddots & & ⋮\\ ⋮& \ddots & v_j& \ddots & ⋮\\ ⋮& & \ddots & \ddots & 0\\ 0& \dots & \dots & 0& v_N\end{array}\right].$$(62)

Also, we introduce the matrix of gyroscopic coefficients $$B\equiv \left[\begin{array}{ll}{\Lambda }^{-1}{B}_{\text{Φ},\text{Φ}}\hfill & {\Lambda }^{-1}{B}_{\text{Φ},\text{Ψ}}\hfill \\ N^{-1}{B}_{\text{Ψ},\text{Φ}}\hfill & N^{-1}{B}_{\text{Ψ},\text{Ψ}}\hfill \end{array}\right],$$(63)

where the block matrices Β_φ,φ, Β_φ,ψ, Β_ψ,φ, and Β_ψ,ψ are defined in terms of the β_k,j defined in Eqs. (40)–(43) as $$B_{\text{Φ},\text{Φ}}\equiv \left[\begin{array}{ccccc}{\beta }_{1,1}& \dots & {\beta }_{1,j}& \dots & {\beta }_{1,M}\\ ⋮& & ⋮& & ⋮\\ {\beta }_{k,1}& \dots & {\beta }_{k,j}& \dots & {\beta }_{k,M}\\ ⋮& & ⋮& \ddots & ⋮\\ {\beta }_{M,1}& \dots & {\beta }_{M,j}& \dots & {\beta }_{M,M}\end{array}\right],$$(64) $$B_{\text{Φ},\text{Ψ}}\equiv \left[\begin{array}{ccccc}{\beta }_{1,-1}& \dots & {\beta }_{1,-j}& \dots & {\beta }_{1,-N}\\ ⋮& & ⋮& & ⋮\\ {\beta }_{k,-1}& \dots & {\beta }_{k,-j}& \dots & {\beta }_{k,-N}\\ ⋮& & ⋮& & ⋮\\ {\beta }_{M,-1}& \dots & {\beta }_{M,-j}& \dots & {\beta }_{M,-N}\end{array}\right],$$(65) $$B_{\text{Φ},\text{Ψ}}\equiv \left[\begin{array}{ccccc}{\beta }_{-1,1}& \dots & {\beta }_{-1,j}& \dots & {\beta }_{-1,M}\\ ⋮& & ⋮& & ⋮\\ {\beta }_{-k,1}& \dots & {\beta }_{-k,j}& \dots & {\beta }_{-k,M}\\ ⋮& & ⋮& & ⋮\\ {\beta }_{-N,1}& \dots & {\beta }_{-N,j}& \dots & {\beta }_{-N,M}\end{array}\right]=-{\overline{B_{\text{Φ},\text{Ψ}}}}^{\text{T}},$$(66) $$B_{\text{Ψ},\text{Ψ}}\equiv \left[\begin{array}{ccccc}{\beta }_{-1,-1}& \dots & {\beta }_{-1,-j}& \dots & {\beta }_{-1,-N}\\ ⋮& & ⋮& & ⋮\\ {\beta }_{-k,-1}& \dots & {\beta }_{-k,-j}& \dots & {\beta }_{-k,-N}\\ ⋮& & ⋮& & ⋮\\ {\beta }_{-N,-1}& \dots & {\beta }_{-N,-j}& \dots & {\beta }_{-N,-N}\end{array}\right].$$(67)

Finally, the matrix accounting for the gravitational and pressure effects of the oceanic surface displacement, A, the matrix describing the coupling between the forcing tidal gravitational potential and the oceanic eigenmodes, R, and the solution vector, X, are respectively expressed as $$A\equiv \left[\begin{array}{cc}{A}_{\text{Φ}}& 0\\ 0& 0\end{array}\right],R\equiv \left[\begin{array}{c}{R}_{\text{Φ}}\\ 0\end{array}\right],X\equiv \left[\begin{array}{c}\Phi \\ \Psi \end{array}\right].$$(68)

The temporal tidal equations given by Eqs. (38) and (39) thus lead to the solution $$X=H\left(\tilde{\sigma },{\tilde{\sigma }}_{\text{R}},{\tilde{\sigma }}_{\text{G}},{\theta }_{\text{c}},{\widehat{\theta }}_{\text{S}};{\gamma }_{\text{D};l}^{\sigma },{\gamma }_{\text{T};l}^{\sigma }\right){\tilde{U}}_{\text{T}},$$(69) $$H\equiv {\tilde{\sigma }}_{\text{G}}^2{\left[\tilde{\sigma }\left(\tilde{\sigma }-i{\tilde{\sigma }}_{\text{R}}\right)I-i\tilde{\sigma }B-{\tilde{\sigma }}_{\text{G}}^{2}A\right]}^{-1}R,$$(70)

with I being the identity matrix. In the above expression, the oceanic tidal response to the tidal perturbation in the Fourier domain is expressed as a complex frequency-dependent matrix (H), which links the tidal gravitational potential (${\tilde{U}}_{\text{T}}$) to the potential and stream functions that describe the horizontal tidal displacement (X).

2.8 Tidal Love numbers, torque, and power

To close this theoretical section, it remains to introduce the parameters that quantify the tidally dissipated energy and its effect on the evolution of planetary systems. The tidal mass redistribution itself is quantified by self-attraction variations, which is derived from the gravitational potential of the tidally distorted body, U_D, defined at its surface as $$U_{\text{D}}=V_0^2\Re \left\{{\displaystyle \sum _{\sigma }{\displaystyle \sum _{l=0}^{+\infty }{\displaystyle \sum _{m=-l}^l{\tilde{U}}_{\text{D};l}^{m,\sigma }{\widehat{Y}}_l^m\left(\widehat{\theta },\widehat{\phi }\right){\text{e}}^{i\sigma t}}}}\right\}.$$(71)

The components ${\tilde{U}}_{\text{D};l}^{m,\sigma }$ of this gravitational potential are expressed as $${\tilde{U}}_{\text{D};l}^{m,\sigma }=k_l^{\sigma }{\tilde{U}}_{\text{T};l}^{m,\sigma }+\left(1+k_{\text{L};l}^{\sigma }\right)\frac{3}{2l+1}\frac{{\rho }_{\text{w}}}{\overline{\rho }}{\tilde{\zeta }}_l^{m,\sigma },$$(72)

where ${\tilde{U}}_{\text{T};l}^{m,\sigma }$ is the corresponding component of the forcing tidal potential introduced in Eq. (44), and the corresponding component of the normalised surface elevation, $${\tilde{\zeta }}_l^{m,\sigma }={\displaystyle \sum _{k=1}^{\infty }\langle {\widehat{Y}}_l^m,{\text{Φ}}_k\rangle }{\lambda }_k{p}_k.$$(73)

The first term of Eq. (72) is the self-attraction variation due to the distortion of the solid part generated by the gravitational tidal force, while the second term is the contribution of the ocean loading, which includes both the potential generated by the oceanic mass redistribution and the potential resulting from the distortion of the solid part induced by ocean loading.

The complex Love numbers are the transfer functions accounting for the intrinsic response of an extended body undergoing a tidal gravitational forcing (e.g. Ogilvie 2014). They are defined, for each SPH, as the ratio of the gravitational potential of the tidal response to the forcing tidal potential evaluated at the body’s surface, $$k_{\text{D};l}^{m,\sigma }\left(\tilde{\sigma },{\tilde{\sigma }}_{\text{R}},{\tilde{\sigma }}_{\text{G}},{\theta }_{\text{c}},{\widehat{\theta }}_{\text{S}};{\gamma }_{\text{D};l}^{\sigma },{\gamma }_{\text{T};l}^{\sigma }\right)\equiv \frac{{\tilde{U}}_{\text{D};l}^{m,\sigma }}{{\tilde{U}}_{\text{T};l}^{m,\sigma }},$$(74)

which, using the expression of the tidal potential given by Eq. (72), yields $$k_{\text{D};l}^{m,\sigma }=k_l^{\sigma }+\left(1+k_{\text{L};l}^{\sigma }\right)\frac{3}{2l+1}\frac{{\rho }_{\text{w}}}{\overline{\rho }}\frac{{\tilde{\zeta }}_l^{m,\sigma }}{{\tilde{U}}_{\text{T};l}^{m,\sigma }}.$$(75)

Considering the above equation, we note that the Love number depends both on the order m and the degree l of the associated SPH in the general case contrary to the solid Love numbers, which only depend on the degrees (see Eqs. (47) and (48)). This is due to the fact that the oceanic tidal response is not spherically symmetric owing to the geometry of coastlines and Coriolis effects.

The long-term effect of tides on the evolution of the planet-perturber system is quantified by the generated tidal torques. The main contributor to this evolution is the time-averaged tidal torque exerted about the spin axis of the planet (Zahn 1966), $${𝒯}_z\equiv \langle {\displaystyle {\int }_{𝒱}{\partial }_{\widehat{\phi }}{U}_{\text{T}}{\delta }_{\rho }\text{d}V}\rangle ,$$(76)