© The Authors 2026

1 Introduction

Tides, both those raised on the central planet and those on natural satellites, are key drivers of the long-term evolution of planetary systems. The amount of energy dissipated due to satellite tides¹ governs the tidal heating rate of the moons’ interiors, and therefore plays a crucial role in their thermal evolution (e.g., Peale & Cassen 1978; Peale et al. 1979). Additionally, the combined effects of satellite and planet tides define the moons’ orbital migration and circularisation rates (e.g., Kaula 1964; Goldreich & Soter 1966), and influence the system’s rotational dynamics (e.g., Goldreich & Gold 1963; MacDonald 1964; Efroimsky & Williams 2009). The present orbit, rotational state, and interior properties of the satellite, on the other hand, determine how much energy is currently dissipated due to tides (see Tobie et al. 2025, and references therein). The system’s dynamical and interior evolution is therefore driven by this intricate feedback between the moons’ orbit, rotation, and tidal response.

Reconstructing the orbits of natural satellites, for instance in ephemeris determination studies, thus requires accounting for the effects of tides on the moons’ dynamics, but also provides us with a natural way to quantify tidal dissipation, both within the central planet or its synchronous satellite(s). Our present knowledge of tidal dissipation parameters beyond the Earth–Moon system primarily comes from astrometry-based constraints on the secular evolution of the moons’ orbits (e.g., Lainey et al. 2009; Lainey 2016; Lainey et al. 2020; Jacobson 2022). From a dynamical modelling perspective, these studies, which investigate long timescales compared to the satellite’s orbital period, only require the averaged long-term effects of planet and satellite tides to be accurately accounted for. Critically, this does not necessitate directly modelling the tidal deformation of the satellite’s gravity field, nor the instantaneous tidal effects on the system’s dynamics. A separate tidal acceleration added to the satellite’s equations of motion is instead sufficient to capture the long-term signature of tidal dissipation in the moons’ dynamics (e.g., Lainey et al. 2007; Lari 2018, see Section 3.4), and therefore perfectly suitable for the determination of ephemerides and dissipation parameters based on direct observations of the moons’ orbits.

However, incorporating tracking data from planetary missions into such ephemeris estimations requires a shift in our modelling paradigm. Spacecraft tracking data are sensitive to tidal signatures both in the spacecraft’s and in the moons’ dynamics (e.g., Magnanini et al. 2024; Fayolle 2025). Whenever in orbit or performing a flyby, the spacecraft’s trajectory is affected by the gravitational response of the tidally deforming moon, but is also influenced by the moon’s own orbit around the central planet, and therefore by tidal effects present in the moon’s dynamics. Accounting for the combined signatures of these different tidal effects in the spacecraft data in a fully consistent manner, critical for a robust estimation of tidal parameters, ideally requires using a single unified model (Fayolle et al. 2022).

However, as mentioned above, tidal effects on the moons’ orbits are typically modelled as an additional perturbing acceleration based on the averaged effect of tides, circumventing the need to model the time-varying tidal deformation of the moon’s gravitational potential. Nonetheless, explicit modelling of the latter becomes necessary to incorporate tides in the spacecraft’s dynamics. In practice, this is typically achieved by introducing gravity coefficient variations, which are defined by the moon’s instantaneous state and orientation (Petit et al. 2010; Iess et al. 2012; Goossens et al. 2024; Park et al. 2025). In principle, directly accounting for satellite tides in the moon’s (time-varying) gravity field in such a way also allows us to model their feedback on the moon’s own dynamics, thus offering a single, coherent way to encapsulate both effects. However, unique challenges arise when attempting to consistently model the impact on the orbital dynamics of the deformation of a synchronous satellite’s gravity field due to tides, and most critically the effects of such gravitational deformation on the satellite’s own orbit and rotation (Fayolle 2025).

In particular, the uniqueness of the challenges arising in the spin-orbit resonance case stems from the leading tidal forcing modes becoming commensurable with the orbital and rotational frequencies (and its harmonics) (e.g., Efroimsky & Williams 2009; Efroimsky 2012). This implies stringent consistency requirements for the satellite’s orbit, rotation, and gravity field tidal deformation, to ensure that this commensurability is maintained at the required accuracy level. Even a small inconsistency in the delicate interplay between translational and rotational dynamics when modelling synchronous satellite tides would indeed lead to unphysical dynamics, either breaking the spin-orbit resonance or mis-estimating the moon’s orbital migration rate (see Magnanini et al. 2026).

Starting from Kaula’s tidal potential expansion (Kaula 1964), numerous past analytical studies have provided a thorough characterisation of satellite tides, including (among others) a theoretical quantification of the amount of energy they dissipate (e.g., Efroimsky & Makarov 2014; Frouard & Efroimsky 2017; Efroimsky 2018), and of their effects on the moon’s orbit (Boué & Efroimsky 2019) and rotation (Efroimsky & Williams 2009; Efroimsky 2012; Williams & Efroimsky 2012; Makarov & Efroimsky 2013). These analytical developments provide us with invaluable insights into the expected manifestations, effects, and signatures of satellite tides. They, however, cannot be directly incorporated in dynamical models of synchronous satellites when fitting real ground-based and/or space-based data. When propagating the moon’s translational dynamics (as required in ephemeris studies and tidal parameter estimation analyses), the exact tidal forcings are no longer known a priori. They are instead influenced by the satellite’s integrated orbit, and are thus an outcome of the numerical propagation. The orbit-rotation-tides interplay then becomes sensitive to the satellite’s (numerically integrated) dynamical response to instantaneous perturbations and, most critically, to any inconsistency in the modelling of that response. This directly conflicts with the strict orbit-rotation-tidal response consistency requirement discussed above for the spin-orbit resonance case. What is required to meet the needs presented by (future) Earth- and space-based tracking data is a model for the tidal gravity field variation of a synchronously rotating satellite that can automatically and consistently adapt itself to variations in orbit and rotation during a numerical propagation of the dynamics.

To address the aforementioned challenges when modelling the gravitational deformation of a synchronous satellite undergoing tidal forcing, we therefore propose a coupled, concurrent propagation of the satellite’s orbit, rotation, and gravity field deformation. This approach follows the framework adopted, for instance, in Correia et al. (2014) and Boué et al. (2016) to study exo-planetary systems. The tidal deformation of the satellite’s gravity field, represented by its spherical harmonics coefficients, is described by an ordinary differential equation and numerically integrated alongside the translational and rotational equations of motion. Such a unified modelling approach automatically accounts for all dynamical couplings and feedbacks at play, while ensuring that the tidally induced gravity deformation remains consistent with the satellite’s orbit, rotation, and assumed rheology, at any instant of the propagation.

As mentioned above, ensuring the consistency between the orbit, tides, and rotation is extremely challenging with classical approaches where they are typically modelled separately, in a decoupled manner. With our coupled model, on the other hand, the self-consistency between the orbital, rotational, and tidal dynamics is naturally maintained throughout the propagation. The ability of our model to produce realistic dynamics, however, becomes contingent upon the determination of a suitable and consistent initial state (to be addressed in Section 5.2). This is a critical challenge when adopting such a coupled modelling strategy: any inconsistency upon initialisation would affect the orbit-rotation-tide interplay and result in a self-consistent yet physically unrealistic solution.

Our paper describes this extended coupled dynamical model and its implementation in a simplified two-dimensional case. As a conceptual validation to motivate further developments for ephemerides applications, we used this simplified propagation test case to show the applicability of such a coupled approach for planet – synchronous satellite systems. In particular, this paper demonstrates our model’s self-consistent representation of the tides-orbit-rotation interplay, while discussing the limitations of currently typical decoupled approaches for the particular case of tides in synchronous satellites. To this end, we leveraged insights from the tidal potential theory to justify the need for our coupled approach, and provided an analytical framework against which to compare our results. We emphasise that the present work presents a proof-of-concept implementation of this unified, fully consistent propagation, which can be used to scrutinise its behaviour and assess its suitability, but is not yet applicable to real data fitting.

Classical models (i.e., tidal effects incorporated as a perturbing acceleration for the moon’s dynamics or as time variations of the moon’s gravity field for the spacecraft’s dynamics) are perfectly suitable when separately considering tidal signatures either in the moon’s or in the spacecraft’s orbit, respectively. Typical natural satellite ephemerides applications (e.g., Lainey et al. 2009, 2012, 2020) and geodetic parameter estimations (e.g., Williams et al. 2014; Lemoine et al. 2013; Durante et al. 2019; Park et al. 2025) are therefore largely unaffected by the modelling complications outlined above. The latter only arise when concurrently modelling and/or extracting the combined tidal signatures in the two types of dynamics, and only for a synchronous satellite.

However, this concurrent modelling of tidal effects (and the associated challenges mentioned above) has become essential for current and upcoming analyses, given the unprecedented size, diversity, and quality of the available datasets. For the upcoming Juice (Jupiter Icy Moons Explorer) (Grasset et al. 2013; Van Hoolst et al. 2024) and Europa Clipper (Pappalardo et al. 2024) missions, tidal dissipation signatures in the moons’ own dynamics and in the spacecraft’s orbit are expected to be both detectable from the spacecraft tracking data, given the expected accuracy level. Even in cases where the dissipation occurring within the moon is too weak for its effect to be visible in spacecraft tracking data, the unified propagation model proposed in this work can still be desirable. By naturally accounting for all dynamical inter-dependencies, such a model prevents moon ephemeris and spacecraft orbit determination errors from affecting estimates of geodetic or rotation parameters. Decoupled analysis of past missions such as Cassini (Durante et al. 2019) have identified possible inconsistencies between orbit, rotation, and gravity solutions (including tidal response) as potential sources of error or discrepancy. Adopting a fully coupled and consistent dynamical model, on the contrary, would facilitate physically realistic and statistically robust estimation solutions. In addition to Juice and Europa Clipper, other current and upcoming missions such as MMX (Martian Moons eXploration) or Hera, but also future icy moon missions (NASA’s Uranus Orbiter and Probe (UOP) or ESA’s L4 concept mission to Enceladus) could benefit from the coupled propagation framework proposed here.

Progressively building the complete analytical framework necessary for the set-up and interpretation of our coupled dynamical model, we start by briefly describing the translational and rotational dynamics in Section 2. Section 3 then describes tides in a planet – synchronous satellite system, including an overview of relevant aspects of tidal theory, expected tidal effects on the satellite’s orbit and rotation, as well as existing modelling strategies and their potential limitations. Section 4 proposes a different tidal modelling strategy and introduces the deformation model used to numerically propagate the gravity coefficient variations. Section 5 then summarises the complete formulation of our coupled orbital-rotational-tidal model and describes the propagation set-ups for both the Earth–Moon and Mars–Phobos systems, which we use as test cases for our model. We concurrently propagated the satellite’s orbital dynamics, rotational motion, and tidal deformation in both systems to verify the behaviour of our coupled model, and its modelling of tidal dissipation effects in particular. In Section 6, we present the results of these coupled propagations and we investigate them in Section 7 in light of relevant insights from the tidal theory outlined in Section 3. Finally, the conclusions, implications, and foreseen future developments are discussed in Section 8.

2 Translational-rotational dynamics

As a necessary background for the tidal theory presented in Section 3 and its application in the present context, we first provide an overview of the well-known equations of motion governing the translational and rotational dynamics of a planet–satellite system. We restrict our models to the gravitational interactions between the planet and satellite. Since our work focusses on modelling the effects of satellite tides on the system’s dynamics, perturbations due to third-body interactions and non-conservative forces are not considered here.

2.1 Translational dynamics

Following Lainey et al. (2004), the inertial (i.e., non-rotating frame origin) acceleration of body i with respect to body j due to its gravitational interaction with j, noted ${\ddot{r}}_{ji}^{(\to \iota )}$, can be decomposed into point mass and extended body interactions: $${\ddot{r}}_{ji}^{(\to \iota )}={\ddot{r}}_{ji}^{(\stackrel{}{¯}\to \overline{\iota })}+{\ddot{r}}_{ji}^{(\stackrel{}{^}\to \overline{\iota })}+{\ddot{r}}_{ji}^{(\stackrel{}{¯}\to \widehat{\iota })}+{\ddot{r}}_{ji}^{(\stackrel{}{^}\to \widehat{\iota })},$$(1)

where $\overline{k}$ and $\widehat{k}$ respectively designate the point mass and extended body components of the gravitational potential of body k. The final term corresponds to figure-figure interactions, and can typically be safely neglected in planet–satellite dynamics (e.g., Lainey et al. 2004; Dirkx et al. 2016). We continue to do so in this paper, but high-accuracy applications for (for instance) Pho-bos may require their reintroduction (Dirkx et al. 2019). Using Newton’s third law, Eq. (1) can be rewritten as $${\ddot{r}}_{ji}^{(\to \iota )}={\ddot{r}}_{ji}^{(\stackrel{}{¯}\to \overline{\iota })}+{\ddot{r}}_{ji}^{(\stackrel{}{^}\to \overline{ı})}-\frac{m_j}{m_i}{\ddot{r}}_{ij}^{(\widehat{ı}\to \stackrel{}{¯})},$$(2)

where the last term accounts for the effect of body i’s extended gravitational potential on body j’s point mass.

The acceleration exerted on a point-mass body i by body j in an inertial frame is given by $${\ddot{r}}_{ji}^{(j\to \overline{ı})}=R_{I/j}{\nabla }_j{U}_j\left(r_{ji}\right),$$(3)

where r_ji is the position of body i with respect to j, and R_I/j the rotation matrix from the body j-fixed reference frame to the inertial frame I. The gravitational potential of j in the body-fixed frame can be expressed as the following spherical harmonics expansion: $$\begin{array}{l}{U}_J(r)=U_{\stackrel{}{¯}}(r)+U_{\stackrel{}{^}}(r)=\frac{{\mu }_j}{r}+\frac{{\mu }_j}{R_J}{\displaystyle \sum _{l=1}^{\infty }{\left(\frac{R_j}{r}\right)}^{l+1}}\\ {\displaystyle \sum _{m=0}^l{P}_{lm}}(\sin \varphi )\left(C_{lm}^j\cos m\lambda +S_{lm}^j\sin m\lambda \right).\end{array}$$(4)

Here r, λ, ϕ are the spherical coordinates of the point at which the potential is evaluated expressed in the body j-fixed reference frame. $C_{lm}^j$ and $S_{lm}^j$ are body j’s spherical harmonics gravity field coefficients at degree l and order m. The values of the spherical harmonics coefficients $C_{lm}^j$ and $S_{lm}^j$ are tied to the body-fixed frame definition. In other words, one can freely modify the orientation of the body-fixed frame without affecting the physical modelling of the dynamical model at hand, provided that the spherical harmonics coefficients are rotated accordingly². We later exploit this freedom in the choice of body-fixed frame and gravity coefficients definition to ensure that our propagation setup follows typical body-fixed frame conventions for synchronous satellites (see Section 5).

In (most of) the rest of this work, we limit ourselves to the two-dimensional case (ϕ = 0), and truncate the potential U_ȷ at degree two. The former assumption implies perfect alignment between the equatorial planes of both the perturbing and perturbed bodies and their mutual orbital plane. This is a reasonable simplification for our analysis, especially given the (very) small inclination and obliquity of natural satellites’ orbits, as it does not hinder our investigation of the coupled model’s behaviour and in particular its ability to reproduce key tidal dissipation features in the system’s dynamics (see Section 5). Focussing on the degree-two terms, on the other hand, is motivated by the fact that tidal effects are strongly dominated by the body’s degree-two response to gravitational forcing. Under these two assumptions, combining Eqs. (3) and (4) leads to $$\begin{array}{l}{\ddot{r}}_{ji}^{(\to \overline{ı})}=R_{I/j}(-\frac{{\mu }_j}{r_{ji}^2}{\widehat{r}}_{ji}+\frac{3}{2}\frac{{\mu }_j{R}_j^2}{r_{ji}^4}{C}_{20}^j{\widehat{r}}_{ji}\\ -9\frac{{\mu }_j{R}_j^2}{r_{ji}^4}\left(C_{22}^j\cos 2\lambda +S_{22}^j\sin 2\lambda \right){\widehat{r}}_{ji}\\ -6\frac{{\mu }_j{R}_j^2}{r_{ji}^4}\left(C_{22}^j\sin 2\lambda -S_{22}^j\cos 2\lambda \right){\widehat{t}}_{ji}),\end{array}$$(5)

where ${\widehat{r}}_{ji}$ and ${\widehat{t}}_{ji}$ designate the radial and tangential unit vectors of body i’s position in the body j-fixed reference frame.

Finally, as the propagation takes place in a non-inertial body j-centred frame (inertial orientation, but non-inertial origin), the acceleration of body i in that frame, noted ${\ddot{r}}_i^{(\to \iota )}$, should also include the inertial acceleration of body j due to body i, as follows: $${\ddot{r}}_i^{(\to \iota )}={\ddot{r}}_{ji}^{(\to \iota )}-{\ddot{r}}_{ij}^{(\iota \to )}.$$(6)

The two inertial acceleration terms on the right-hand side can be obtained by combining Eqs. (2), (3), and (5).

The presence of the rotation matrix R_I/j and body-fixed longitude λ in Eq. (5) highlights the need to model the bodies’ rotations even when focussing on translational dynamics. In practice, typical orbital studies of planet–moon systems numerically integrate the translational equations of motion only, and rely on kinematic rotation models. The latter describe the per-turber and/or perturbed bodies’ rotations either as a function of time or via an analytical representation of the body’s orientation based on its current translational state.

In such analytical models, the long-axis of the satellite is typically considered as pointing towards the orbit’s empty focus. For convenience of implementation in orbital dynamics-focussed studies, the determination of the empty focus’ location α with respect to the perturber commonly relies on the following approximation (e.g., Lainey et al. 2019) $$\alpha =-2e\frac{r·\stackrel{.}{r}}{|r\times \stackrel{.}{r}|}+O\left(e^2\right).$$(7)

The limitations of this particular approximation and of kinematic rotation models in general in terms of orbit-rotation consistency when modelling tides will be elaborated upon in Section 3.6. As a first step towards a more consistent modelling alternative, the next section presents the rotational equations of motion necessary to the concurrent numerical integration of the orbital-rotational dynamics.

2.2 Rotational dynamics

The orientation of a body i, in this work described using the quaternion q_i representation, obeys the following (for detail on quaternion definition, see Fukushima 2008): $${\stackrel{.}{q}}_i=Q\left(q_i\right){\omega }_i.$$(8)

Here ω_i is body i’s angular velocity vector in the body-fixed frame of body i and Q is defined as $$Q(q)=\frac{1}{2}\left(\begin{array}{ccc}-q_1& -q_2& -q_3\\ q_0& -q_3& q_2\\ q_3& q_0& -q_1\\ -q_2& q_1& q_0\end{array}\right).$$(9)

The evolution of the rotational rate is, in turn, given by the Euler equation: $$\frac{\text{d}\left(I_i{\omega }_i\right)}{\text{d}t}+{\omega }_i\times \left(I_i{\omega }_i\right)={\displaystyle \sum _k{\Gamma }_k},$$(10)

the right-hand side being the sum of all external torques Γ_k acting on body i. I_i is its inertia matrix, which relates to the gravity field coefficients as $$I_i=m_i{R}_i^2\left(\begin{array}{ccc}\frac{C_{20}^i}{3}-2C_{22}^i& -2S_{22}^i& -C_{21}^i\\ -2S_{22}^i& \frac{C_{20}^i}{3}+2C_{22}^i& -S_{21}^i\\ -C_{21}^i& -S_{21}^i& -\frac{2C_{20}^i}{3}\end{array}\right)+{\overline{I}}_i{1}_{3\times 3},$$(11)

where ${\overline{I}}_i$ and 1_3×3 respectively represent the body’s mean moment of inertia and the three-by-three identity matrix.

In the simple two-dimensional case that we here consider, the rotation axis and principal axis of inertia are aligned, resulting in $C_{21}^i=S_{21}^i=0$. When considering tidal deformation, the gravity coefficients and thus the inertia matrix are no longer constant but time-dependent (see more details in Section 4). The evolution of the body’s angular velocity is thus given by $$I_i{\stackrel{.}{\omega }}_i+{\stackrel{.}{I}}_i{\omega }_i+{\omega }_i\times \left(I_i{\omega }_i\right)={\displaystyle \sum _k{\Gamma }_k}.$$(12)

We moreover restrict our dynamical model to the mutual gravitational interaction between body i and the central body j (again ignoring figure-figure interactions), such that the only torque acting on body i is the one exerted by j on body i’s extended gravity field: $${\Gamma }_g=-m_j{r}_{ji}\times {\nabla }_i{U}_i\left(-r_{ji}\right),$$(13)

which after plugging in the expression for body i’s degree-two gravitational potential eventually becomes $${\Gamma }_g=\frac{6m_j{\mu }_i}{R_i}{\left(\frac{R_i}{r_{ji}}\right)}^3\left(C_{22}^i\sin 2\lambda -S_{22}^i\cos 2\lambda \right)\widehat{h},$$(14)

with $\widehat{h}$ the orbit normal unit vector.

3 Tides in a planet – synchronous satellite system

This section provides an overview of the modelling of tides in planetary system’s dynamics, with a specific focus on their effects on synchronous satellites’ orbits and rotations. Section 3.1 first describes where tidal deformation comes into play in the orbital and rotational dynamics. Section 3.2 then presents Kaula’s formulation for the tidal potential, including its expansion as a superimposition of different forcing modes. Derived expressions for the tidal torque in the specific spin-orbit resonance case are presented in Section 3.3. Section 3.4 finally introduces a simplified force formulation typically used to model the effects of tides on the (synchronous) orbits of natural satellites (Section 3.5). Having provided the necessary tidal modelling theory, Section 3.5 describes the effects of tides on the satellite’s dynamics, which Section 3.6 builds upon to discuss the physical implications of different modelling strategies. This is critical to (i) identify the complications that arise in the 1:1 spin-orbit resonance case; (ii) understand the need for a coupled integration of tides, orbits, and rotations. The analytical results and insights derived in the following section will moreover be key to the validation of our numerical propagation (see details in Section 5.4).

3.1 Incorporation of tides in the translational-rotational dynamics

The forcing exerted by the tide-raising body causes the gravitational potential of the perturbed body to deform over time. This time variation can be directly incorporated in the coefficients of the spherical harmonics expansion of the deforming body’s gravitational potential (Eq. (4)). In the rest of this work, we use ∆C, S_lm(t) to represent the (time-dependent) variations of the gravity coefficients due to tides (to be added to the body’s static gravity field, see details in Section 5.2.2). The influence of tidal deformation on the system’s translational and rotational dynamics then naturally follows from the time-varying gravity coefficients entering the expressions for the gravity acceleration and torque, respectively given by Eqs. (5) and (14). The inertia matrix I also varies along with the perturbed body’s gravity coefficients (Eq. (11)), and this time dependency (and subsequent non-zero time derivative $\dot{I}$) should also be accounted for in the Euler equations (Eq. (10)).

Introducing time variations in the coefficients defining the gravitational potential’s spherical harmonics expansion is a common way to incorporate tidal effects in translational and/or rotational dynamics. Basic expressions for these time dependent coefficients are usually limited to a single (dominating) frequency (e.g., Petit et al. 2010). The tidal potential theory developed in the following subsection includes a rigorous derivation of more complete, multi-frequency expressions for the tidally driven variations of the perturbed body’s gravity coefficients (see Section 3.2.2).

3.2 Tidal potential and mode decomposition

In the following, we describe the tides raised by a perturbing body (represented by the asterisk superscript^⋆ following the convention introduced by Kaula 1964) at a surface point on a perturbed body (i.e., body subject to tidal deformation). The coordinates of this surface point R = (R, ϕ, λ) are expressed as spherical coordinates in the perturbed body-fixed reference frame. Unless otherwise indicated, all quantities with no asterisk superscript implicitly refer to the perturbed body.

3.2.1 Tidal perturbing potential

At a surface point R, the perturber raises the following tidal potential (e.g. Kaula 1964): $$\begin{array}{l}{W}_t\left(R,r^{\star }\right)=-\frac{{\mu }^{\star }}{r^{\star }}{\displaystyle \sum _{l=2}^{\infty }{\left(\frac{R}{r^{\star }}\right)}^l}{\displaystyle \sum _{m=0}^l\frac{(l-m)!}{(l+m)!}}\left(2-{\delta }_{0m}\right)\\ P_{lm}(\sin \varphi )P_{lm}\left(\sin {\varphi }^{\star }\right)\cos m\left(\lambda -{\lambda }^{\star }\right).\end{array}$$(15)

Here R is the radius of the perturbed body and µ^⋆ the gravitational parameter of the perturber. r^⋆ = (r^⋆, ϕ^⋆, λ^⋆) represents the position of the perturber with respect to the perturbed body, expressed using spherical coordinates in the perturbed body-fixed frame. P_lm are the associated Legendre polynomials of degree l and order m.

Exploiting Kaula’s potential expansion of spherical coordinates to orbital elements (Kaula 1961), the tide-inducing potential given by Eq. (15) can be decomposed as follows (e.g., Kaula 1964): $$\begin{array}{l}{W}_t\left(r,r^{\star }\right)=-\frac{{\mu }^{\star }}{R}{\displaystyle \sum _{l=2}^{\infty }{\left(\frac{R}{a^{\star }}\right)}^{l+1}}{\displaystyle \sum _{m=0}^l\frac{(l-m)!}{(l+m)!}}\left(2-{\delta }_{0m}\right)\\ P_{lm}(\sin \varphi ){\displaystyle \sum _{p=0}^l{F}_{lmp}}\left(i^{\star }\right){\displaystyle \sum _{q=-\infty }^{\infty }{G}_{lpq}}\left(e^{\star }\right)\cos \left({\beta }_{lmpq}^{\star }-m\lambda \right).\end{array}$$(16)

a^⋆, e^⋆, and i^⋆ respectively designate the averaged semi-major axis, eccentricity, and inclination of the perturbing body, and F_lmp and G_lpq are known as inclination and eccentricity functions (see Table C.1). For planet–satellite tides (which this study focusses on), these orbital elements describe the mutual orbit of the central planet and its satellite irrespective of which object is the perturbing or the perturbed body. For the sake of notation conciseness, we can therefore drop the asterisk superscript for the orbital elements. The combinations are here noted β_lmpq (as in Frouard & Efroimsky 2017) and are defined as follows (Kaula 1964)³: $${\beta }_{lmpq}\equiv (l-2p)\tilde{\omega }+(l-2p+q)M+m\Omega -m\theta ,$$(17)

where $\tilde{\omega }$, Ω, and M are the averaged argument of periapsis, right ascension, and mean anomaly of the perturber. θ designates the rotation angle of the perturber around its principal axis of inertia. These combinations are associated with the following forcing modes: $$\begin{array}{l}{\omega }_{lmpq}={\dot{\beta }}_{lmpq}\\ \equiv (l-2p)\dot{\tilde{\omega }}+(l-2p+q)\dot{M}+m\dot{\Omega }-m\dot{\theta }.\end{array}$$(18)

As addressed in detail in Frouard & Efroimsky (2017), ω_lmpq can only be considered as forcing modes if they are secular (i.e., their own time variations happen at much slower rates, $\left|{\dot{\omega }}_{lmpq}/{\omega }_{lmpq}\right|\ll \left|{\omega }_{lmpq}\right|)$. While fulfilled for most bodies, the aforementioned condition presents an additional complication in the case of tides raised on a natural satellite in 1:1 spin-orbit resonance (as will be discussed below).

The rotation angle of a synchronous satellite can be expressed as $$\theta =M+\tilde{\omega }+\Omega +\gamma +\pi ,$$(19)

where γ represents the satellite’s instantaneous (longitudinal) physical libration, i.e., the deviation of the satellite’s orientation from exact synchronicity due to perturbations (caused by non-zero eccentricity and/or inclination, or third body perturbations). The orientation of the satellite-fixed frame here follows classical conventions for synchronous satellites (i.e., long axis pointing towards the central planet at apo- and periapsis), as illustrated in Figure 1, and will be further discussed in Section 5.2.2. In our simplified, two-dimensional planet–moon dynamical system, the libration spectrum is limited to harmonics of the orbital frequency. The physical longitudinal libration is then dominated by a one-per-orbit variation of amplitude $\mathcal{A}$ and can therefore be expressed as $$\gamma \approx -\mathcal{A}\sin M.$$(20)

Other harmonics would appear when considering third body and out-of-plane perturbations. However, the commensurability between the tidal forcing modes, and the librational modes corresponding to harmonics of the orbital frequency makes this particular subset of librational forcings especially relevant for tidal modelling. For a rigid and solid body, the once-per-orbit libration amplitude is given by (e.g., Willner et al. 2010) $$\mathcal{A}=6e\frac{\sigma }{3\sigma -1},$$(21)

where σ is function of the principal moments of inertia (diagonal elements of I): $$\sigma =\frac{I_{yy}-I_{xx}}{I_{zz}}.$$(22)

The minus sign in front of Eq. (20) is purely a convention choice, and the opposite libration definition is also often found in the literature (e.g., Frouard & Efroimsky 2017; Efroimsky 2018). Introducing Eqs. (19) and (20) into (18), the quantities ω_lmpq then become $$\begin{array}{l}(1:1){\omega }_{lmpq}\equiv (l-m-2p)\dot{\tilde{\omega }}+(l-m-2p+q)\dot{M}\\ +mn\mathcal{A}\cos M,\end{array}$$(23)

where the last libration-driven term is no longer secular and instead varies at a once-per-orbit rate.

Resolving this requires expanding Eq. (16) to generalise its applicability to librating bodies, by exploiting Bessel functions to decompose the librational variation(s) onto secular components, using the expansion $$\begin{array}{l}\left\{\begin{array}{c}\cos \\ \sin \end{array}\right\}\left({\beta }_{lmpq}+m\mathcal{A}\sin M\right)\\ ={\displaystyle \sum _{s=-\infty }^{\infty }{J}_s}(-m\mathcal{A})\left\{\begin{array}{c}\cos \\ \sin \end{array}\right\}\left({\beta }_{lmpq}-sM\right),\end{array}$$(24)

where $J_s(-m\mathcal{A})$ are Bessel functions of order s. Applying the above to Eq. (16), an augmented potential formulation accounting for libration(s) can be obtained, as first derived in Frouard & Efroimsky (2017): $$\begin{array}{l}{W}_t\left(r,r^{\star }\right)=-\frac{{\mu }^{\star }}{R}{\displaystyle \sum _{l=2}^{\infty }{\left(\frac{R}{a}\right)}^{l+1}}{\displaystyle \sum _{m=0}^l\frac{(l-m)!}{(l+m)!}}\left(2-{\delta }_{0m}\right)\\ P_{lm}(\sin \varphi ){\displaystyle \sum _{p=0}^l{F}_{lmp}}(i){\displaystyle \sum _{q=-\infty }^{\infty }{G}_{lpq}}(e)\\ {\displaystyle \sum _{s=-\infty }^{\infty }{J}_s}(-m\mathcal{A})\cos \left({\beta }_{lmpqs}-m\lambda \right).\end{array}$$(25)

The expanded combinations β_lmpqs and related quantities ω_lmpqs (which are now secular and therefore valid forcing modes) are given by $${\beta }_{lmpqs}\equiv (l-m-2p)\tilde{\omega }+(l-m-2p+q-s)M$$(26) $$\begin{array}{l}{\omega }_{lmpqs}\equiv (l-m-2p)\dot{\tilde{\omega }}+(l-m-2p+q-s)\dot{M}\\ \approx (l-m-2p+q-s)\dot{M}\end{array}$$(27)

In the absence of librations (i.e., $\mathcal{A}=0$) and outside the specific spin-orbit resonance case, Eqs. (25)–(27) only contain the s = 0 term and simplify back to the initial expansion (Eqs. (16)–(18)).

Fig. 1 Schematic representation of the pointing direction of the satellite’s long axis θ and thus of the orientation of the satellite-fixed frame, including the contribution of the longitudinal physical libration γ. The satellite is represented at periapsis (M = 0), at apoapsis (M = π), and at a arbitrarily chosen location along its orbit. In this work, we define the satellite-fixed reference frame such that the satellite’s long axis points towards the planet at apo- and periapsis.

3.2.2 Mode decomposition of the tidal response

According to the linear theory of tides, the perturbed body’s response to a given tidal forcing mode is independent to forcings at other frequencies (provided that the body does not have lateral heterogeneities, see Rovira-Navarro et al. 2024), and can be described by a phase-lagged linear deformation. A purely linear deformation describes the body’s response to the static tide, while the phase lag accounts for the dynamical tide (i.e., not purely elastic response, see for example Efroimsky & Lainey 2007). Both the amplitude and phase of the deformation are frequency-dependent, as the body’s response to forcing (in terms of both how large the resulting deformation is, and how long it takes for the body to deform) depends on the frequency at which the forcing occurs, a behaviour described by the body’s rheology (see Section 4). The perturbed body’s response to the perturbing potential given by Eq. (25) therefore raises the following tide-induced potential, here evaluated at an external point r = (r, ϕ, λ) with r > R: $$\begin{array}{l}{U}_t\left(r,r^{\star }\right)=-\frac{{\mu }^{\star }}{R}{\displaystyle \sum _{l=2}^{\infty }{\left(\frac{R}{r}\right)}^{l+1}}{\left(\frac{R}{a}\right)}^{l+1}{\displaystyle \sum _{m=0}^l\frac{(l-m)!}{(l+m)!}}\\ \left(2-{\delta }_{0m}\right)P_{lm}(\sin \varphi ){\displaystyle \sum _{p=0}^l{F}_{lmp}}(i){\displaystyle \sum _{q=-\infty }^{\infty }{G}_{lpq}}(e)\\ {\displaystyle \sum _{s=-\infty }^{\infty }{J}_s}(-m\mathcal{A})k_l\left({\omega }_{lmpqs}\right)\cos \left({\beta }_{lmpqs}-{ϵ}_{lmpqs}-m\lambda \right).\end{array}$$(28)

Here the subscript t indicates that U_t is the tide-induced perturbation of the gravitational potential, and not the total gravitational potential U as in Eq. (28). The (frequency-dependent) static Love numbers k_l(ω_lmpqs) and phase lags ϵ_lmpqs describe the perturbed body’s gravity field deformation caused by the static and dynamical tides, respectively. The phase lags ϵ_lmpqs are often found alternatively parametrised by (frequency-dependent) tidal quality factors Q(ω_lmpqs), defined as follows (e.g., Efroimsky 2012): $$1/Q\left({\omega }_{\text{lmpqs\hspace{0.17em}}}\right)=\sin \left|{ϵ}_{\text{lmpqs\hspace{0.17em}}}\right|.$$(29)

In the simplified two-dimensional, degree-two response framework, the above expression can be further simplified as $$\begin{array}{l}{U}_{t,l=2}\left(r,r^{\star }\right)=-\frac{{\mu }^{\star }}{R}{\left(\frac{R}{r}\right)}^3{\left(\frac{R}{a}\right)}^3{\displaystyle \sum _{m=0}^2\frac{(2-m)!}{(2+m)!}}\left(2-{\delta }_{0m}\right)\\ P_{2m}(0){\displaystyle \sum _{p=0}^2{F}_{2mp}}(0){\displaystyle \sum _{q=-\infty }^{\infty }{G}_{2pq}}(e){\displaystyle \sum _{s=-\infty }^{\infty }{J}_s}(-m\mathcal{A})\\ k_2\left({\omega }_{2mpqs}\right)\cos \left({\beta }_{2mpqs}-{ϵ}_{2mpqs}-m\lambda \right).\end{array}$$(30)

In the following and for the sake of conciseness, we drop the subscript l = 2. The above formulation and all subsequent developments can however be easily expanded to higher degrees if required.

All inclination functions F_lmp(i) are O(i) except for $F_{201}(i)=\frac{1}{2}+O(i)$ and F₂₂₀(i) = 3 + O(i) (e.g., Efroimsky & Williams 2009; Boué & Efroimsky 2019). This allows us to further simplify the above expression by only keeping the terms lmpqs = 201qs (order zero response U_t,m=0) and lmpqs = 220qs (order two response U_t,m=2): $$U_t\left(r,r^{\star }\right)=U_{m=0}\left(r,r^{\star }\right)+U_{m=2}\left(r,r^{\star }\right),$$(31) $$\begin{array}{l}{U}_{t,m=0}\left(r,r^{\star }\right)=-\frac{1}{4}\frac{{\mu }^{\star }}{R}{\left(\frac{R}{r}\right)}^3{\left(\frac{R}{a}\right)}^3{\displaystyle \sum _{q=-\infty }^{\infty }{G}_{21q}}(e)\\ {\displaystyle \sum _{s=-\infty }^{\infty }{J}_s}(0)k_2\left({\omega }_{201qs}\right)\cos \left({\beta }_{201qs}-{ϵ}_{201qs}\right)\end{array}$$(32) $$\begin{array}{l}{U}_{t,m=2}\left(r,r^{\star }\right)=-\frac{3}{4}\frac{{\mu }^{\star }}{R}{\left(\frac{R}{r}\right)}^3{\left(\frac{R}{a}\right)}^3{\displaystyle \sum _{q=-\infty }^{\infty }{G}_{20q}}(e)\\ {\displaystyle \sum _{s=-\infty }^{\infty }{J}_s}(-2\mathcal{A})k_2\left({\omega }_{220qs}\right)\cos \left({\beta }_{220qs}-{ϵ}_{220qs}-2\lambda \right).\end{array}$$(33)

The dominant forcing modes are defined as (Eq. (18)) $${\omega }_{201q0}=q\dot{M}-2\dot{\theta },$$(34) $${\omega }_{220q0}=(2+q)\dot{M}-2\dot{\theta }.$$(35)

or, in the specific spin-orbit resonance case (i.e., for tides raised on a synchronous satellite) (Eq. (27)) as $${\omega }_{201qs}={\omega }_{220qs}=(q-s)\dot{M}=(q-s)n.$$(36)

The corresponding eccentricity functions G_21q(e) and G_20q(e) are provided in Table C.1. The forcing mode ω22000 is the only term in O(e⁰) and therefore typically dominates the body’s response. If the perturber is in a 1:1 spin-orbit resonance, however, this leading mode vanishes (q = s = 0 leading to ω₂₂₀₀₀ = 0 in Eq. (36), while ${\omega }_{22000}=2\dot{M}-2\dot{\theta }\ne 0$ outside of the spin-orbit resonance). This is at the core of the complications arising when modelling the tides raised on a synchronous satellite, as will be discussed in Section 3.5.

Considering Eqs. (31)–(33) and the expression for the perturbed body’s static gravitational potential (Section 2.1, Eq. (4)), the tide-induced potential can alternatively be expressed as time-dependent variations to be added to the perturbed body’s static gravity field coefficients. For the degree-two response and in the two-dimensional case, the expected variations of the (unnor-malised) cosine and sine coefficients ∆C_2m and ∆S_2m are given by $$\begin{array}{l}\Delta C_{20}=\frac{1}{2}\frac{{\mu }^{\star }}{\mu }{\left(\frac{R}{a}\right)}^3{\displaystyle \sum _{q=-\infty }^{\infty }{G}_{20q}}(e)\\ {\displaystyle \sum _{s=-\infty }^{\infty }{J}_s}(0)k_2\left({\omega }_{201qs}\right)\cos \left({\beta }_{201qs}-{ϵ}_{201qs}\right),\end{array}$$(37) $$\begin{array}{l}\Delta C_{22}=\frac{1}{4}\frac{{\mu }^{\star }}{\mu }{\left(\frac{R}{a}\right)}^3{\displaystyle \sum _{q=-\infty }^{\infty }{G}_{20q}}(e)\\ {\displaystyle \sum _{s=-\infty }^{\infty }{J}_s}(-2\mathcal{A})k_2\left({\omega }_{220qs}\right)\cos \left({\beta }_{220qs}-{ϵ}_{220qs}\right),\end{array}$$(38) $$\begin{array}{l}\Delta S_{22}=\frac{1}{4}\frac{{\mu }^{\star }}{\mu }{\left(\frac{R}{a}\right)}^3{\displaystyle \sum _{q=-\infty }^{\infty }{G}_{20q}}(e)\\ {\displaystyle \sum _{s=-\infty }^{\infty }{J}_s}(-2\mathcal{A})k_2\left({\omega }_{220qs}\right)\sin \left({\beta }_{220qs}-{ϵ}_{220qs}\right),\end{array}$$(39)

(see Section 3.1 on how to introduce these coefficient variations into the translational and rotational equations of motion).

The potential expansion leading to the above expressions for the deforming body’s gravity coefficients still relies on a number of simplifying assumptions which are important to keep in mind. First, Kaula’s expansion assumes that the semi-major axis a, eccentricity e, and inclination i of the perturber are constant (for expansion to time-varying eccentricity, see e.g., Walker & Rhoden 2022). Moreover, as highlighted above for the spin-orbit resonance case in the presence of librations, the quantities ω_lmpqs, defined by the orbital elements and rotation rates $\dot{\tilde{\omega }},\dot{\Omega },\dot{M}$, and $\dot{\theta }$ (Eq. (18)), can only be considered as forcing modes if they are secular (i.e., exhibit very slow variations). The two aforementioned conditions can be problematic for satellites trapped in mean motion resonances (MMRs) (e.g., Jupiter’s Galilean moons) because of the (rapid) periodic orbital element variations that the resonances induce. Finally, our present libration model only includes a single harmonic at the satellite’s orbital frequency (Eq. (20)). Further expansion of Eq. (25) to account for additional harmonics in the satellite’s librational response are discussed in Section 7.1 and Appendix D, but would come at the cost of increased complexity and computational load upon evaluation of the tidal potential.

3.3 Tidal torque

In this section, we focus on the tidal torque caused by tides raised on a synchronous satellite. This specific case requires more attention due to the vanishing of the leading tidal mode ω22000 and contribution of the physical libration(s), as high-lighted in Section 3.2. The torque induced by satellite tides is moreover essential to our understanding of the effect of satellite tides on the satellite’s own orbit (see Section 3.5). Similar derivations would provide equivalent expressions for the simpler case of planet tides where the leading tidal mode ω₂₂₀₀₀ does not vanish.

Starting from the libration-compatible expansion of the tidal potential given by Eq. (30), one obtains the following expression for the tidal torque, as done in Frouard & Efroimsky (2017). Because the torque only acts along the orbit normal, we directly provide its out-of-plane component: $$\begin{array}{l}{\Gamma }_t=-m^{\star }\left(r\times \nabla U_t\left(r,r^{\star }\right)\right)·\widehat{h}\\ =-\frac{G{m}^{\star 2}}{a}{\left(\frac{R}{a}\right)}^5{\displaystyle \sum _{m=1}^{2}2}m\frac{(2-m)!}{(2+m)!}{\displaystyle \sum _{p=0}^2{F}_{2mp}}(i){\displaystyle \sum }_{q=-\infty }^{\infty }\end{array}$$(40) $$\begin{array}{l}{G}_{2pq}(e){\displaystyle \sum _{p^{\prime }=0}^2{F}_{2m{p}^{\prime }}}(i){\displaystyle \sum _{q^{\prime }=-\infty }^{\infty }{G}_{2p^{\prime }{q}^{\prime }}}(e){\displaystyle \sum _{s=-\infty }^{\infty }{\displaystyle \sum _{s^{\prime }=-\infty }^{\infty }{J}_s}}(-m\mathcal{A})\\ J_{s^{\prime }}(-m\mathcal{A})k_2\left({\omega }_{2mpqs}\right)\sin \left({\beta }_{2mpqs}-{\beta }_{2m{p}^{\prime }{q}^{\prime }{s}^{\prime }}-{ϵ}_{2mpqs}\right),\end{array}$$(41)

with m ≥ 1 as the order zero response does not induce a torque. Keeping only the leading terms lmp = 220 (see Section 3.2.2 and eccentricity functions listed in Table C.1) with p^′ = p = 0, the above eventually becomes $$\begin{array}{l}{\Gamma }_t=-\frac{3}{2}\frac{G{m}^{\star 2}}{a}{\left(\frac{R}{a}\right)}^5{\displaystyle \sum _{q=-\infty }^{\infty }{G}_{20q}}(e){\displaystyle \sum _{q^{\prime }=-\infty }^{\infty }{G}_{20q^{\prime }}}(e)\\ {\displaystyle \sum _{s=-\infty }^{\infty }{\displaystyle \sum _{s^{\prime }=-\infty }^{\infty }{J}_s}}(-2\mathcal{A})J_{s^{\prime }}(-2\mathcal{A})\\ k_2\left({\omega }_{220qs}\right)\sin \left({\beta }_{220qs}-{\beta }_{220q^{\prime }{s}^{\prime }}-{ϵ}_{220qs}\right).\end{array}$$(42)

The secular component of the tidal torque ${\overline{\Gamma }}_t$ can be recovered by only retaining the terms that fulfil the following condition: $${\beta }_{220qs}={\beta }_{220q^{\prime }{s}^{\prime }}\ne 0.$$(43)

Other terms are periodic and cancel out over one orbital period. The second part of the above condition follows from the forcing mode ω_220qs vanishing if β_220qs = 0. Recalling the tidal forcing modes definition in the specific spin-orbit resonance case (Eq. (36)), this translates to⁴ $$q-s=q^{\prime }-s^{\prime }\ne 0.$$(44)

Summing all the terms fulfilling the above condition and simplifying the resulting expression eventually leads to $${\overline{\Gamma }}_t=18\frac{G{m}^{\star 2}{R}^5}{a^6}{k}_2\left[e^2\sin {ϵ}_n+\frac{\mathcal{A}e}{2}\sin {ϵ}_n\right],$$(45)

with ϵ_n the tidal phase lag evaluated at the once-per-orbit forcing frequency ω_lmpqs = n. The first term of the secular torque Γ¯ corresponds to the zero-libration case, and the total tidal torque including the libration’s contribution can thus be re-written as $${\overline{\Gamma }}_t=\left(1+\frac{\mathcal{A}}{2e}\right){\overline{\Gamma }}_{\gamma =0},$$(46)

where ${\overline{\Gamma }}_{t,\gamma =0}$ is the tidal torque in the absence of libration.

3.4 Alternative force formulation

When including tidal effects in the orbital dynamics of a planet–satellite system, a commonly used approach is to translate the tidal potential (Eq. (28)) into a tidal force acting as an additional perturbation on the satellite, both for planet and satellite tides (e.g., Lainey et al. 2007; Lari 2018) (see discussion in Section 3.6 for more details on why such an alternative is necessary). The expression for this tidal force, evaluated at a point r from the perturbed body is the following (see derivation in Section 2.3.2 in Fayolle 2025): $$\begin{array}{l}F\left(r,r^{\star }\right)=-m^{\star }\nabla U_t\left(r,r^{\star }\right)\\ =-\frac{3}{2}G{m}^{\star 2}{R}^5{k}_2{\left(\frac{1}{r^{\star }}\right)}^2{\left(\frac{1}{r}\right)}^3\\ \left[\left(1-5{\left(r^{\star }·r\right)}^2\right)\widehat{r}+2\left(r^{\star }·r\right)\widehat{t}\right].\end{array}$$(47)

To model the effect of tides on the mutual orbit of the perturbed and perturbing bodies, F must be evaluated at the per-turber’s position: for planet tides, the tidal force represents the additional acceleration acting on the satellite due to the tides it raises on its central planet; for satellite tides, it is the reaction to the force acting on the planet due to planet-raised tides on the satellite. To account for the dynamic tide (i.e., not purely elastic response of the perturbed body), an important step in the typical derivation of the tidal force lies in the following relation between the perturber’s position when the tidal potential is raised r, and at the moment of the forcing r^⋆ (Mignard 1980): $$r^{\star }=r-\Delta t(\stackrel{.}{r}-\omega \times r)+O\left(\Delta t^2\right).$$(48)

Here ω is the rotation vector of the perturbed body and ∆t the so-called tidal time lag, i.e., the amount of time required for that body to deform due to the tidal forcing. ∆t relates to the tidal phase lag as $$\Delta t=\frac{T\left|ϵ\right|}{2\pi },$$(49)

with ϵ the phase lag at the main forcing frequency and T the period of this forcing. An important underlying simplification of this tidal force formulation is that it assumes a single frequency forcing and response, as well as a constant time lag, therefore neglecting the influence of the satellite’s eccentricity (so-called constant time lag model). Circling back to the general definition of the forcing modes (Eq. (18)), the nominal leading forcing mode is ${\omega }_{22000}\equiv 2n-2\dot{\theta }$, resulting in the following tidal period: $$T=\frac{2\pi }{\left|{\omega }_{22000}\right|}=\frac{2\pi }{|2\dot{M}-2\dot{\theta }|}.$$(50)

This for instance corresponds the forcing period for the tides raised on the planet. For the specific case of the spin-orbit resonance (Eq. (27)), the leading forcing modes become ω₂₂₀₁₀ = n and ω_220,−1,0 = −n, such that the tidal period is $$T=\frac{2\pi }{\left|{\omega }_{22010}\right|}=\frac{2\pi }{\left|{\omega }_{220,-1,0}\right|}=\frac{2\pi }{n}.$$(51)

Using Eq. (48) to expand (1/r^⋆)³ as a O(∆t²) expression eventually leads to the following expression for the tidal force, after decomposing the above into radial and tangential components and assuming the orbit’s normal and rotation pole are coplanar (for a complete derivation, see Mignard 1980), $$F=-3G{m}^{\star 2}{R}^5{k}_2\frac{1}{r^7}\left(\widehat{r}+\Delta t\frac{3\dot{r}}{r}\widehat{r}-\Delta t(\dot{\theta }-\dot{v})\widehat{t}\right),$$(52)

with $\dot{v}$ the time derivative of the true anomaly (i.e., the orbital velocity).

The typical formulation for the tidal force as given in Eq. (52) relies on a few important simplications:

1. a truncation at O(∆t²) (Eq. (48)), which might not be valid for large eccentricities;

2. a single-frequency forcing and response, as evidenced by the definition of the tidal time lag in Eq. (49) (equivalent to assuming a frequency-independent response to the leading term(s) in Eq. (30)).

While a more rigorous, frequency-dependent expression for the tidal force could be derived from Eq. (30), Eq. (52) has been widely and successfully used to model the effects of both planet and satellite tides in natural ephemerides applications, and proven suitable to reproduce the main features of tidally driven orbital evolution (Lainey et al. 2007, 2009, 2012, 2019, 2020). The effects of planet and satellite tides on their mutual dynamics, and the limitations of this tidal force formulation to model these effects in the specific spin-orbit resonance case will be elaborated upon in Section 3.5.

From the tidal force expression in Eq. (52), one can moreover recover the tidal torque derived from the tidal potential in Eq. (45), including the additional contribution of the physical libration. Recalling that the rotation rate of a synchronous satellite in the presence of one main once-per-orbit longitudinal libration is $$\dot{\theta }=n+\frac{\mathcal{A}}{2e}(n-\dot{v}),$$(53)

the tangential component of the tidal force (from which the tidal torque originates) can be re-written as $$F_{\text{tan}}=3G{m}^{\star 2}{R}^5{k}_2\frac{1}{r^7}\Delta t\left(1+\frac{\mathcal{A}}{2e}\right)(n-\dot{v})\widehat{t}$$(54) $$={\left(1+\frac{\mathcal{A}}{2e}\right)F_{\text{tan}}|}_{\gamma =0},$$(55)

with F_tan|_γ=0 the tangential component of the tidal force in the absence of physical librations. This matches the librational contribution obtained with the tidal potential approach (Eq. (46)).

3.5 Effects of planet and satellite tides on the satellite’s orbit

The loss of energy through tidal dissipation drives the evolution of a natural satellite’s orbit, and more specifically the changing rate of its semi-major axis and eccentricity. The literature provides well-known approximations for the secular evolution of a and e, with an important distinction between planet and satellite tides (e.g., Kaula 1964; Goldreich & Soter 1966; Peale & Cassen 1978; Peale 1999; Malhotra 1991; Lainey et al. 2009; Lari 2018; Emelyanov 2018; Boué & Efroimsky 2019). For the former, i.e., raised on the planet by the satellite, the resulting change in the satellite’s semi-major axis and eccentricity can be approximated by⁵ $$\frac{\text{d}a}{\text{d}t}=3\frac{m^{\star }}{m}\frac{k_2}{Q}{\left(\frac{R}{a}\right)}^{5}na,$$(56) $$\frac{\text{d}e}{\text{d}t}=\frac{57}{8}\frac{m^{\star }}{m}\frac{k_2}{Q}{\left(\frac{R}{a}\right)}^{5}ne.$$(57)

For the case of tides raised by the planet on a synchronous satellite, however, two conflicting secular variations of the moon’s semi-major axis can be found (e.g., Goldreich & Soter 1966; Lainey et al. 2009 vs Emelyanov 2018; Boué & Efroimsky 2019), as reported below⁶, $$(1)\frac{\text{d}a}{\text{d}t}=-21\frac{m^{\star }}{m}\frac{k_2}{Q}{\left(\frac{R}{a}\right)}^{5}na{e}^2,$$(58a) $$(2)\frac{\text{d}a}{\text{d}t}=-57\frac{m^{\star }}{m}\frac{k_2}{Q}{\left(\frac{R}{a_s}\right)}^{5}na{e}^2,$$(58b)

whereas they are consistent for the eccentricity changing rate: $$\frac{\text{d}e}{\text{d}t}=-\frac{21}{2}\frac{m^{\star }}{m}\frac{k_2}{Q}{\left(\frac{R}{a}\right)}^{5}ne.$$(59)

This apparent discrepancy for the effects of satellite tides finds its origin in the commensurability between the moon’s orbital and rotational periods, combined with the need to conserve rotational energy for the spin-orbit resonance to be maintained. A detailed explanation can be found in Magnanini et al. (2026), but the reasoning will be briefly reconstructed here, as it brings insights into the modelling challenges associated with satellite tides.

Evaluating the orbital energy dissipated through satellite tides in the case of exact synchronicity leads to the semi-major axis evolution given in Eq. (58b), a result that is equivalently obtained starting from the simplified force formulation in Eq. (52) or from the complete tidal potential expansion (Eq. (30)) (for a rigorous derivation of the orbital elements evolution from Kaula’s tidal potential expansion, see Boué & Efroimsky 2019; Magnanini et al. 2026). As shown in Section 3.3, however, exact 1:1 spin-orbit synchronicity implies a non-zero secular tidal torque (Eq. (45)). A slightly super-synchronous rotation rate can be computed for which the tidal torque cancels out over one orbit (Levrard 2008; Wisdom 2008), $$\dot{\theta }=\left(1+6e^2\right)n,$$(60)

which can thus seem a more stable rotational state (no residual tidal torque spinning up the satellite’s rotation). For this pseudo-synchronicity case, the amount of tidal energy dissipation is no longer consistent with Eq. (58b), but one instead recovers the more commonly used approximation given by Eq. (58a) (see Section 2.3.4 in Fayolle 2025)⁷.

However, the exactly synchronous rotation of the Solar System’s natural satellites contradicts Eq. (60), and the physical realism of such a pseudo-synchronous equilibrium for terrestrial bodies was indeed refuted (e.g., Makarov & Efroimsky 2013). For natural satellites to maintain exact spin-orbit resonance, the secular tidal torque must nonetheless be compensated, as it would otherwise progressively speed up the satellite’s rotation out of spin-orbit synchronicity.

In practice, this torque balance is maintained via a slight re-orientation of the satellite (e.g. as theorised in Yoder & Peale 1981; Murray & Dermott 1999; Williams & Efroimsky 2012; Makarov & Efroimsky 2013): the body’s triaxiality causes a non-zero secular torque if the principal axis of inertia and the direction of the planet at periapsis are slightly misaligned (see Eq. (14) and Figure 1). To conserve spin-orbit resonance, the averaged torque exerted by the central planet on the non-spherical satellite ${\overline{\Gamma }}_g$ due to this misalignment (Eq. (14)) and the secular torque caused by satellite tides ${\overline{\Gamma }}_t$ (Eq. (45)) should therefore cancel out, $${\overline{\Gamma }}_g+{\overline{\Gamma }}_t=0$$(61) $$\frac{6m^{\star }m{R}^2}{a^3}{S}_{22}=18\frac{G{m}^{\star 2}{R}^5}{a^6}{k}_2{e}^2\sin {ϵ}_n\left(1+\frac{\mathcal{A}}{2e}\right),$$(62)

thus leading to the following equilibrium S₂₂ value (e.g., Magnanini et al. 2026): $$S_{22}=3\frac{{\mu }^{\star }}{\mu }{\left(\frac{R}{a}\right)}^3{k}_2\sin {ϵ}_n{e}^2\left(1+\frac{\mathcal{A}}{2e}\right).$$(63)

This non-zero S₂₂ ensures constant rotational energy, but also causes a secular drift of the satellite’s semi-major axis, which exactly compensates for the discrepancy between Eqs. (58b) and (58a). This counterbalancing re-orientation of the satellite here translates into a non-zero S 22 because we choose to keep our body-fixed convention unchanged (i.e., satellite’s x-axis pointing towards the planet at apo- and periapsis, see Figure 1 and details in Section 5.2.2). As discussed in Section 2.1, it could equivalently be represented by a small rotation of the satellite-fixed frame⁸. Nevertheless, the counteracting triaxial torque caused by the satellite’s re-orientation and its subsequent effects on the dynamics are invariant upon different body-fixed frame conventions. Evaluating Eq. (63) for Io (the only Moon for which we have an estimate for k₂/Q = k₂ sin |ϵ_n| = 0.015, Lainey et al. 2009), the conservation of rotational energy requires a balancing S₂₂ = 1.30 · 10⁻⁹, equivalent to an orientation offset of 6.67 · 10⁻⁵ rad.

The total semi-major axis rate of the satellite is therefore in agreement with Eq. (58a), irrespective of any assumption on exact or pseudo-synchronicity. The sole difference is that in the physically realistic case of fully synchronous satellites, only part of the total secular change in a directly results from the effect of tides on the orbit, while the other part is an indirect effect of the conservation of rotational energy and the resulting non-zero S 22 (and is thus oftentimes overlooked when focussing on orbital evolution).

Finally, it should be noted that Eqs. (58a) and (59) approximate the secular change in the satellite’s semi-major axis and eccentricity caused by satellite tides in the absence of librations. The amount of orbital energy dissipated when accounting for the influence of the main once-per-orbit longitudinal libration is $$\Delta E_{\gamma }=\left(1+\frac{4}{7}\frac{\mathcal{A}}{e}+\frac{1}{7}\frac{{\mathcal{A}}^2}{e^2}\right)\Delta E_{\gamma =0},$$(64)

with ∆E_γ and ∆E_γ=0 representing the energy dissipation with and without libration, respectively (see detailed derivation in Efroimsky 2018). This eventually leads to $$\left\{\begin{array}{l}{\frac{\text{d}a}{\text{d}t}|}_{\gamma }\hfill \\ {\frac{\text{d}e}{\text{d}t}|}_{\gamma }\hfill \end{array}\right\}=\left(1+\frac{4}{7}\frac{\mathcal{A}}{e}+\frac{1}{7}\frac{{\mathcal{A}}^2}{e^2}\right)\left\{\begin{array}{l}{\frac{\text{d}a}{\text{d}t}|}_{\gamma =0}\hfill \\ {\frac{\text{d}e}{\text{d}t}|}_{\gamma =0}\hfill \end{array}\right\}$$(65)

for the secular evolution of a and e caused by tides raised on a librating synchronous satellite, with ${\frac{\text{d}a}{\text{d}t}|}_{\gamma =0}$ and ${\frac{\text{d}e}{\text{d}t}|}_{\gamma =0}$ given by Eqs. (58a) and (59), respectively.

3.6 Implications for satellite tides modelling

As demonstrated in Section 3.5, accurately representing the effect of satellite tides on the dynamical evolution of a synchronous satellite requires the consistent modelling of its orbit and rotation for (i) the main tidal forcing mode ω₂₂₀₀₀ to vanish (Eq. (36)); (ii) the spin-orbit resonance to be maintained (Eqs. (61)–(63)). These two conditions are inter-dependent, and breaking one of them would eventually lead to breaching the other. Maintaining both above conditions is, however, very difficult to achieve in practice when relying on classical kinematic rotation models (i.e., analytical representations of the satellite’s rotation, with the possible incorporation of multi-frequency perturbations, see Section 2.1), as will be discussed in the following.

Due to the vanishing of the main tidal forcing mode ω22000 in Eq. (33), the secular evolution of the satellite’s orbital elements caused by satellite tides is dominated by linear functions of e² (Eqs. (58a) and (59)). This implies that any model simplification that incurs an error in O(e²) is no longer acceptable and would violate the two aforementioned conditions (i.e., ω22000 = 0 while maintaining the resonance). This O(e²) consistency requirement is well illustrated by the conservation of rotational energy issue exposed in Section 3.5 and its implications to both maintain the spin-orbit resonance and predict the right secular evolution of the satellite’s semi-major axis (discrepancy between Eqs. (58a) and (58b)). The equilibrium between the secular tidal torque (Eq. (45)) and the restoring triaxial torque (Eq. (61)), and its consequences (additional contribution to the semi-major axis change rate and conservation of exact- rather than pseudo-synchronicity, Eq. (60)) are indeed all O(e²) effects, whose delicate balance is contingent upon the coherent modelling of the intricate coupling between orbit and rotation.

In practice, this level of orbit-rotation consistency is extremely difficult to achieve with a kinematic representation of the satellite’s rotation. Such analytical rotation models typically assume that the satellite’s long axis points towards the empty focus (Section 2.1). This assumption is actually a O(e²) approximation in itself (e.g., Murray & Dermott 1999), whose practical implementation oftentimes relies on an additional O(e²) simplification (see Eq. (7)). This demonstrates that relying on the orbit’s empty focus to define the satellite’s rotation is not compatible with the O(e²) modelling consistency requirement necessary to accurately model satellite tides, with the important implications for the predicted orbital migration rate evidenced in Section 3.5. In addition to the inadequacy of such an implementation shortcut, a more fundamental flaw in adopting an analytical rotation model stems from the satellite’s rotational inertia, which prevents its orientation from instantly matching the osculating orbit. The satellite’s rotation indeed follows the long-term evolution of the orbit, rather than adapting to immediate perturbations (see in-depth discussion in Martinez et al. 2025). This is critical to maintain spin-orbit resonance, but makes it extremely challenging to define a suitable time-dependent rotation model while concurrently propagating the (perturbed) orbital dynamics of the satellite.

When modelling the effect of satellite tides on the satellite’s own orbit using the tidal force formulation given in Eq. (52), a commonly used approach allows to circumvent this orbit-rotation consistency issue. It exploits the fact that the energy dissipated by librational tides over one orbit amounts to 4/3 of that dissipated by radial tides (Murray & Dermott 1999). This allows us to translate Eq. (52) into a scaled tidal force acting in the radial direction only, while getting rid of the $\dot{\theta }-\dot{v}$ librational term where the dependency on the satellite’s rotation comes into play (e.g., Lari 2018): $$F=-3G{m}^{\star 2}{R}^5{k}_2\frac{1}{r^7}\left(\widehat{r}+\Delta t\frac{7\dot{r}}{r}\widehat{r}\right).$$(66)

This approach is a very efficient and robust workaround to accurately represent satellite tides in an orbital propagation without the need for a fully consistent rotation, which is precisely what makes the tidal force formulation so appealing for orbital dynamics propagation (see detailed discussion in Fayolle 2025).

By design, this force formulation is, however, only suitable to model and/or estimate the averaged effects of satellite tides on the satellite’s own orbit (over one orbital period). It is, by contrast, inapplicable to represent their influence on the dynamics of a neighbouring spacecraft (whether in orbit, or performing a flyby), or for highly accurate propagations of the moon’s instantaneous dynamics. The latter is typically accounted for by modelling tides as additional time-variations of the deformed body’s gravitational potential, as described in Section 3.2 (e.g., Petit et al. 2010; Iess et al. 2012; Goossens et al. 2024; Park et al. 2025), a modelling approach that is extremely rotation-sensitive when applied to tidal effects on the satellite’s dynamics (Magnanini et al. 2026). Using different models to include tides in the moons’ and spacecraft’s dynamics would nevertheless be prone to discrepancies, especially since satellite tides affect the spacecraft in two ways: via the direct effect of the satellite’s deformed gravitational potential, and via the indirect contribution of the moon’s tidally perturbed orbit on the spacecraft’s motion. For planetary missions relying on a series of moons’ flybys (e.g., Cassini, Juice, Europa Clipper), addressing such a modelling challenge is therefore particularly critical, not only to ensure the consistent propagation of the satellites and spacecraft, but also to ultimately allow for the consistent estimation of tidal (and rotational) parameters from their combined effects on the moons’ and spacecraft’s orbits. This is especially important to avoid mismodelling the satellite’s expected orbital migration rate and therefore under- and over-evaluating the tidal dissipation because of any translational-rotation inconsistency (see Section 3.5).

The unique complications arising in the spin-orbit resonance case, combined with the need to consistently model tidal effects on both the satellite’s and spacecraft’s dynamics, call for a unified way to model the complex interplay between the satellite’s orbit, rotation, and tidal deformation. In particular, this requires accounting for the instantaneous tidal deformation of the satellite’s gravitational potential, in a manner that ensures full consistency between the tidal forcing and subsequent response.

4 Deformation model

Adopting the approach followed by Correia et al. (2014) and Boué et al. (2016) for exoplanetary systems, we propose to concurrently propagate the perturbed body’s translational and rotational dynamics, along with the instantaneous deformation of its gravitational potential caused by tidal effects. Such a three-level dynamical model (orbit, rotation, deformation), by automatically accounting for all couplings at play, ensures the full consistency of the propagation results while reproducing the main dissipation-driven dynamical evolution features (Sections 3.3 and 3.5).

While the complete set of equations (including the translational and rotational dynamics from Section 2) will be re-iterated in Section 5.1, the following section focusses on the missing component: i.e., the deformation model adopted for the satellite’s deformation to dynamical forcing. Section 4.1 first introduces the general deformation modelling framework necessary to describe the tidal deformation of the perturbed body as a differential equation. Section 4.2 then provides the exact formulation for the Maxwell rheology used in the rest of this work.

4.1 General formulation

This section presents the model used to describe the perturbed body’s gravity field deformation induced by a perturbing potential W. Following Correia et al. (2014), the total perturbing potential evaluated at a surface point R is here defined as the sum of the tide-raising potential W_t and centrifugal potential W_c, as $$W(R,t)=W_t(R)+W_c(R)$$(67) $$=W_t(R)+\frac{1}{2}{(}^{\omega },$$(68)

with W_t(R) being given by Eq. (15) and ω the perturbed body’s rotational rate vector. The gravitational potential U^′ induced by the body’s deformation under the above perturbing potential is given by the convolution (e.g., Boué et al. 2016) $$U^{\prime }(R,t)={\left(\frac{R}{r}\right)}^3{k}_2(t)*W(R,t),$$(69)

where k₂(t) is the complex degree-two Love number distribution. The total gravitational potential U^′ differs from the tidally induced potential U_t as our coupled model also accounts for the centrifugal deformation. We here restrict ourselves to the degree-two (l = 2) perturbation and response, but the above could be expanded to higher degrees (see e.g., Boué et al. 2016). It is moreover important to keep in mind that the tidally induced potential resulting from Eq. (69) is actually equivalent to the expression already derived in Section 3.2.2 (Eq. (30)), but presents the critical advantage of not relying on a mode decomposition.

Eq. (69) is generally applicable when modelling gravitational deformation under a perturbing potential, the body’s specific response to different forcings based on its interior properties (i.e., rheology) being described by its Love number distribution k₂(t). It remains valid for all eccentricities (thus circumventing the possible O(e²) inconsistencies mentioned in Section 3.6), and does not rely on a Fourier-like series similar to the one used by Kaula (1964), implying that the dynamical problem at hand needs not be decomposable into such an expansion, which might indeed not be the case for satellites in MMRs (see Section 3.2.2).

4.2 Maxwell rheology

A viscoelastic response is typically assumed for natural bodies, with the Maxwell rheology (equivalent to a viscous damper mounted in series with an elastic spring) being widely used among possible rheological models. This model presents the combined advantages of offering a very simple mathematical formulation, while still adequately representing the main features of tidal dissipation (Boué et al. 2016). In the rest of this work, we therefore assume that both the planets and satellites under consideration behave as Maxwell bodies. The implications of such a modelling assumption will be discussed in Section 8. For a Maxwell body, the Fourier transform of the complex Love number distribution is (e.g., Henning et al. 2009) $$k_2(\omega )=k_2^f\frac{1+i{\tau }_e\omega }{1+i\tau \omega },$$(70)

with $k_2^f$ the body’s fluid Love number of degree two, and τ and τ_e the global and elastic (or Maxwell) relaxation times, respectively.

This Fourier transform turns the convolution in Eq. (69) into the following differential equation for the deformation-induced potential U(R) (adapted from Correia et al. 2014; Boué et al. 2016) $$U^{\prime }(R)+\tau {\dot{U}}^{\prime }(R)=k_2^f{\left(\frac{R}{r}\right)}^3\left(W(R)+{\tau }_e\dot{W}(R)\right)$$(71) $$U^{\prime }(R)+\tau {\dot{U}}^{\prime }(R)=U_{eq}(R)+{\tau }_e{\dot{U}}_{eq}(R),$$(72)

with U_eq the equilibrium deformation-induced potential (i.e., gravitational potential due to the body’s fluid response), defined as $$U_{eq}(R)=k_2^f{\left(\frac{R}{r}\right)}^{3}W(R).$$(73)

The dependency on t is made implicit in Eqs. (71)–(73) compared to Eq. (69).

Recalling the spherical harmonics expansion of the gravitational potential (Eq. (4)), the deformation-induced potential can be represented by time-dependent increments of the perturbed body’s gravity field coefficients. Owing to the orthogonality of spherical harmonic basis functions, the coefficient increments follow the differential equation for the deformation-potential (Eq. (72)): $$\Delta C_{2m}+\tau \Delta {\dot{C}}_{2m}=\Delta C_{2m}^{eq}+{\tau }_e\Delta {\dot{C}}_{2m}^{eq},$$(74) $$\Delta S_{2m}+\tau \Delta {\dot{S}}_{2m}=\Delta S_{2m}^{eq}+{\tau }_e\Delta {\dot{S}}_{2m}^{eq}$$(75)

with ∆C_2m, ∆S_2m and $\Delta C_{2m}^{eq},\Delta S_{2m}^{eq}$ the instantaneous and equilibrium coefficients, respectively. The latter can be directly derived from the equilibrium potential U_eq (Correia et al. 2014): $$\Delta C_{20}^{eq}=-k_f\left[\frac{{\dot{\theta }}^2{R}^3}{3\mu }+\frac{1}{2}\frac{{\mu }^{\star }}{\mu }{\left(\frac{R}{r^{\star }}\right)}^3\right],$$(76) $$\Delta C_{22}^{eq}=\frac{k_f}{4}\frac{{\mu }^{\star }}{\mu }{\left(\frac{R}{r^{\star }}\right)}^3\cos 2{\lambda }^{\star },$$(77) $$\Delta S_{22}^{eq}=\frac{k_f}{4}\frac{{\mu }^{\star }}{\mu }{\left(\frac{R}{r^{\star }}\right)}^3\sin 2{\lambda }^{\star }.$$(78)

Here $\dot{\theta }$ corresponds to the rotational rate, as the rotational velocity vector $\omega ={\left(\begin{array}{lll}0\hfill & 0\hfill & \dot{\theta }\hfill \end{array}\right)}^T$ is perfectly aligned with the orbit’s normal in our simplified two-dimensional case. These equilibrium coefficients act as a forcing in the differential equations driving the time evolution of the gravity deformation (Eqs. (74)–(75)). From Eqs. (77)–(78), their time derivatives can be computed as follows: $$\Delta {\dot{C}}_{20}^{eq}=k_2^f\left[\frac{2\dot{\theta }{R}^3}{3\mu }\ddot{\theta }-\frac{3}{2}\frac{{\mu }^{\star }}{\mu }{\left(\frac{R}{r^{\star }}\right)}^3\frac{{\dot{r}}^{\star }}{r^{\star }}\right],$$(79) $$\Delta {\dot{C}}_{22}^{eq}=-\frac{k_2^f}{4}\frac{{\mu }^{\star }}{\mu }{\left(\frac{R}{r^{\star }}\right)}^3\left[3\frac{{\dot{r}}^{\star }}{r^{\star }}\cos 2{\lambda }^{\star }+2{\dot{\lambda }}^{\star }\sin 2{\lambda }^{\star }\right],$$(80) $$\Delta {\dot{S}}_{22}^{eq}=-\frac{k_2^f}{4}\frac{{\mu }^{\star }}{\mu }{\left(\frac{R}{r^{\star }}\right)}^3\left[3\frac{{\dot{r}}^{\star }}{r^{\star }}\sin 2{\lambda }^{\star }-2{\dot{\lambda }}^{\star }\cos 2{\lambda }^{\star }\right].$$(81)

Using the Darwin-Radau approximation, the fluid Love number and relaxation times are all related to the body’s interior properties (Correia et al. 2014): $$k_2^f=5{\left(1+{\left(\frac{5}{2}-\frac{15I_{zz}}{4m{R}^2}\right)}^2\right)}^{-1}-1,$$(82) $${\tau }_e=\frac{\eta }{\mu },$$(83) $$\tau =\frac{19\eta }{2g\rho R},$$(84)

where ρ is the mean density, µ the shear modulus (i.e., rigidity), and η the viscosity. For a spherical body of constant density, the global relaxation time τ becomes $$\tau =\frac{38\pi R^4\eta }{3G{m}^2}.$$(85)

It should be noted that the above definitions for the fluid Love number and relaxation times assume a homogeneous and incompressible body, and can therefore only describe the effective or averaged response of a body’s (stratified) interior. The implications and limitations of such an assumption will be discussed in Section 8. We nonetheless emphasise that they only affect the specific formulation of the differential equations for the gravity coefficients, and not the general framework of our coupled model. From Eq. (70), the (frequency-dependent) modulus and argument of a Maxwell body’s complex Love number can moreover be linked to the quantities k₂(ω) and ϵ(ω) typically used in tidal potential expansions (Section 3.2) as follows: $$k_2(\omega )=k_2^f\sqrt{\frac{1+{\tau }_e^2{\omega }^2}{1+{\tau }^2{\omega }^2}},$$(86) $$ϵ(\omega )=\arctan \left(\frac{\left(\tau -{\tau }_e\right)\omega }{1+\tau {\tau }_e\omega }\right).$$(87)

This allows us to evaluate the body’s Maxwell response to tidal forcings at the different modes ω_lmpqs, enabling a direct comparison between the results of our numerical propagation and the analytical formalism presented in Section 3.2. The incorporation of the deformational Equations (74)–(75) in the rest of our dynamical model will be discussed in Section 5.1.

It should be stressed that Eqs. (76)–(81), once inserted back in Eqs. (74)–(75), provide a specific formulation of the first order differential equation governing the gravity deformation, valid in our two-dimensional, zero obliquity case with limited perturbations. For high-accuracy applications, this simplified framework could however easily be expanded to more complex configurations. More specifically, generalising the above to three-dimensional problems would require foregoing the zero latitude assumption in the tidal potential formulation (Eq. (28)) and subsequently including order 1 coefficients ∆C₂₁, ∆S₂₁ in the propagation ($\Delta C_{21}^{\text{eq}},\Delta S_{21}^{\text{eq}}\ne 0$ if ϕ^⋆ ≠ 0). Eq. (72) can moreover easily be adapted to incorporate extra perturbations (e.g., third-body effects) in the gravitational deformation forcing in addition to the tidal and centrifugal potentials (thus modifying the equilibrium coefficients definition) and can be expanded to higher degrees.

5 Complete coupled model formulation and set-up

This section provides a detailed description of the propagation set-ups used to implement and test our coupled model, as well as an overview of the dynamical features and quantities to compare to classical models. Section 5.1 first provides the complete set of equations of motion to be integrated in our coupled model framework and highlights their interdependencies. Section 5.2 then presents the strategy required to identify a consistent, physically realistic initial state for our coupled dynamics.

To investigate the behaviour of our coupled translational-rotational-tidal approach, we use both the Earth–Moon and Mars–Phobos systems as test cases, paying special attention to the model’s ability to reproduce the dynamical evolution of a planet – synchronous satellite system caused by satellite tides. Section 5.3 describes the specific dynamical environment and propagation settings used for the two planet–satellite systems. We underline that our analyses do not aim to provide a complete and fully realistic dynamical descriptions of the systems under consideration, but rather to analyse the dynamical properties that naturally emerge from our coupled propagation, with special attention paid to tidal dissipation features. The main dynamical features that we expect our unified model to reproduce, outlining the scope of our model validation, are listed in Section 5.4.

All our simulations have been performed using the TU Delft Astrodynamics Toolbox (Dirkx et al. 2022, 2025). For this work, Tudat has been extended with the numerical integration of the gravity field in Eq. (93), building upon the existing functionality for coupled orbital-rotational dynamics of natural satellites (Martinez et al. 2025). Tudat has previously been applied to various studies in the modelling and estimation of natural satellite dynamics (e.g. Dirkx et al. 2016; Fayolle 2025).

5.1 Complete set of equations of motion

We here combine the orbit, rotation, and gravitational deformation models (separately presented in Sections 2.1, 2.2, and 4, respectively) into the complete set of equations forming our coupled model. The complete propagated state is defined as: $$x={\left(\begin{array}{lllll}r\hfill & i\hfill & q\hfill & w\hfill & z\hfill \end{array}\right)}^T,$$(88)

where r, q, and z respectively represent the body’s translational, rotational, and gravitational states, the latter combining the propagated degree-two coefficients⁹ z = (∆C₂₀ ∆C₂₂ ∆S₂₂)^T (see Section 5.2.2 for the coefficients decomposition into static and time-varying components). Under the gravitational influence of a perturber ⋆, the time derivative of this extended state (in the perturber-centred frame) obeys the following set of equations of motion: $$\stackrel{.}{r}=\stackrel{.}{r};$$(89) $$\begin{array}{l}\stackrel{\mathrm{..}}{r}=\frac{{\mu }^{\star }+\mu }{\mu }{R}_{I/(q)}\nabla U(r,q,z)\\ +\frac{{\mu }^{\star }+\mu }{{\mu }^{\star }}{R}_{I/\star }\left(q^{\star }\right){\nabla }_{\star }{U}_{\widehat{\star }}\left(r^{\star }\right);\end{array}$$(90) $$\stackrel{.}{q}=Q(q)\omega ;$$(91) $$\stackrel{.}{w}=I^{-1}\left(-m^{\star }r\times \nabla U(r,q,z)-\stackrel{.}{I}\omega -\omega \times (I\omega )\right);$$(92) $$\dot{z}=\frac{1}{\tau }\left(z^{eq}-z+{\tau }_e{\dot{z}}^{eq}\right).$$(93)

Following the notations introduced in Section 2.1, ∇U and ${\nabla }_{\star }{U}_{\widehat{\star }}$ represent the gradients of the perturbed and perturbing bodies’ gravitational potentials¹⁰ in their respective body-fixed frames (linked to the inertial frame via the rotation matrices R_I/ and R_I/⋆). Both gravitational potentials U and $U_{\widehat{\star }}$ can be expressed using the formulation given by Eq. (4) (expanded to degree two in this work), but the gravity field coefficients of the deforming body are now time-varying, and are actually described by the gravitational state z. Finally, z^eq designates the equilibrium coefficients vector (Eqs. (76)–(78)).

The coupling between the translational, rotational, and tidal dynamics is immediately apparent from the above set of equations, which all depend on the perturber’s body-fixed position and on the gravity coefficients of the deforming body. In particular, the perturbed gravitational potential of the body undergoing tides U(r, q, z) automatically accounts for the contribution of tidal deformation to the translational and rotational dynamics of the system. The orbital elements’ evolution and emergence of a secular tidal torque therefore naturally follow from the propagation of the tidally perturbed equations of motion (assuming a suitable initial condition of the dynamical system, see Section 5.2.1). This is a crucial distinction with respect to Sections 3.3 and 3.5 where they are separately computed from an additional tidal potential.

It is furthermore important to remark that $\stackrel{.}{I}$ depends on the time derivative $\dot{z}$ of the gravity coefficients, by definition of the inertia matrix (Eq. (11)). Conversely, the time derivative of the equilibrium coefficients ${\dot{z}}^{eq}$ depends on the time derivative of the perturber’s body-fixed longitude ${\dot{\lambda }}^{\star }$ (Eqs. (79)–(81)), which in turn is a function of the time derivative of the angular velocity vector $\dot{\omega }$. These interdependencies between the time derivatives of the rotational and gravitational states should be carefully considered when implementing the above equations, and again highlight the need for a concurrent integration.

As a brief summary of the advantages and implications of the unified modelling approach described above with respect to classical decoupled models:

• our concurrent propagation provides the body’s instantaneous gravitational deformation under the perturbing potential(s) (here tidal and centrifugal), that deformation being inherently consistent with the translational and rotational dynamics of the system. This allows tides to be consistently incorporated in the dynamics of both the satellite and a nearby spacecraft, while automatically accounting for the effects of any variation in the satellite’s rotation;

• all instantaneous forcings (at all frequencies) are automatically incorporated in the equilibrium potential U_eq induced by the deformation (Eq. (73));

• this approach does not involve any series expansion of the tidal potential, and therefore does not introduce any assumption regarding the secular rates of the satellite’ orbital elements, or any series truncation error (especially truncations in eccentricity).

5.2 Obtaining a realistic initial state

While our coupled model automatically ensures that the translational, rotational, and tidal dynamics remain fully consistent during the propagation, its ability of produce realistic dynamics for a given known planet–satellite system depends on the definition of a suitable initial state. This is especially critical in a propagation framework such as the one presented in this paper. As opposed to real data fitting, in a propagation-only set-up, there is no observation constraint forcing the system to follow its present dynamical state. In particular, two important points of attention arise which are specific to our fully coupled approach:

• the propagated dynamics should not be corrupted by the artificial introduction of rotational normal modes (Section 5.2.1);

• the mean gravity coefficients that the coupled model converges to should match our present knowledge of the satellite’s gravity field (Section 5.2.2).

The initialisation process adopted in this work (see detail in Section 5.2.3) relies on a long relaxation phase until the system’s dynamics reach a converged, physically consistent steady state (ensuring that the two above conditions are fulfilled, among others). The numerical integration of the tidal deformation makes this relaxation phase crucial: the propagation first undergoes a transient phase during which the gravity coefficients progressively converge to their mean equilibrium values.

5.2.1 Normal modes damping

Integrating the equations of motion for the rotational dynamics gives rise to resonant oscillations of the body-fixed principal axes, referred to as normal modes (e.g., Borderies & Yoder 1990; Rambaux et al. 2010). The frequency at which these normal modes occur is directly linked to the body’s interior properties. Specifically focussing on the orientation of the body’s long axis in our two dimensional case, the longitudinal normal mode frequency is given by (e.g., Rambaux et al. 2012) $${\omega }_{\text{lon\hspace{0.17em}}}\approx n\sqrt{\frac{3\left(I_{yy}-I_{xx}\right)}{I_{zz}}},$$(94)

with I_xx, I_yy, I_zz the diagonal elements of the inertia matrix (which are function of the degree-two gravity coefficients, see Eq. (11)) and n the orbit’s mean motion.

These resonant oscillations are typically damped over planetary system evolution timescales, and no longer affect the rotational dynamics of the Earth–Moon and Mars–Phobos systems under consideration. Numerically propagating the rotational equations of motion from an insufficiently judiciously defined, undamped initial state will, however, automatically introduce the normal modes in the resulting numerical solution of the equations of motion (Martinez et al. 2025). The interdependency between the normal mode’s frequency (Eq. (94)) and the gravity coefficient values makes the obtention of a damped initial state particularly critical when concurrently propagating a body’s rotation and gravitational deformation.

To make ensure that the non-physical oscillations originating from the normal modes are absent from our integrated rotational solution, we applied a similar damping algorithm as done in Rambaux et al. (2012) or Martinez et al. (2025) (among others). This algorithm relies on the addition of an artificial damping torque to Eq. (10), $${\Gamma }_{\text{damping\hspace{0.17em}}}=-\frac{1}{{\tau }_{\text{d}}}I\left(\omega -{\omega }_0\right),$$(95)

with τ_d a damping timescale. This torque forces the body to maintain its (known) mean rotation ω₀ while damping other rotational responses (for more details, see e.g. Rambaux et al. 2012; Martinez et al. 2025). Once a damped initial state (i.e., no longer affected by normal modes) has been obtained, Γ_damping is removed, and the coupled translational-rotational-tidal dynamics are re-propagated from this damped initialisation.

5.2.2 Gravity deformation conventions and set-up

An important aspect of our model set-up relates to the distinction between the static and varying parts of the perturbed body’s gravity field. The fluid number of a uniform body in hydrostatic should theoretically be equal to 3/2. However, both the Moon and Phobos are known to depart from this exact uniform hydrostatic state. We chose to define the fluid number k_f such that it matches our current knowledge of these moons’ C₂₀ values using Eq. (82). Despite this particular set-up choice, the equilibrium deformation described by Eqs. (76)–(78) cannot exactly reproduce our present estimates of the satellite’s degree-two gravity coefficients. Eqs. (76)–(78) indeed cannot (entirely) account for the non-hydrostatic part of the satellite’s gravity field, even with a scaled k_f value. It should nonetheless be noted that, provided that one can find an appropriate Maxwell rheology law (Eq. (70)) to describe the body’s response to forcing, this limitation does not hinder the ability of our coupled model to produce a realistic deformation response and consistently account for the dynamical signatures of said deformation.

However, to ensure that the deforming body’s gravitational potential is in agreement with its known gravity coefficients once the equilibrium deformation (Eq. (73)) is reached, we defined a base static gravity field to which the deformation propagated with Eqs. (74)–(75) is added. The instantaneous gravity coefficients are thus given by $$C_{2m}(t)=C_{2m}^s+\Delta C_{2m}(t),$$(96) $$S_{2m}(t)=S_{2m}^s+\Delta S_{2m}(t).$$(97)

Here ∆C_2m(t) and ∆S 2_m(t) are the coefficient variations numerically integrated with Eqs. (74)–(75), their equilibrium values being given by Eqs. (76)–(78). $C_{2m}^s$ and $S_{2m}^s$, on the other hand, describe the static part of the gravity field (including the non-hydrostatic component). When initialising our model, these static coefficients were computed by subtracting the equilibrium deformation (Eqs. (76)–(78), evaluated at initialisation) from the body’s known coefficients (noted C_2m and S_2m), as follows: $$C_{2m}^s=C_{2m}-\Delta C_{2m}^{eq}$$(98) $$S_{2m}^s=S_{2m}-\Delta S_{2m}^{eq}.$$(99)