Tidal evolution of the Keplerian elements

TL;DR

This paper derives corrected formulas for the rates of Keplerian orbital elements under tidal perturbations, addressing errors in previous literature and providing expressions applicable to systems with comparable masses.

Contribution

It presents a re-derivation of the tidal evolution equations, correcting previous formulas and clarifying the contributions to orbital inclination changes for binary systems.

Findings

Derived new expressions for da/dt, de/dt, di/dt including higher order terms.

Identified and corrected errors in Kaula's classical formulas for systems with comparable masses.

Highlighted the importance of angular momentum conservation in orbital inclination evolution.

Abstract

We address the expressions for the rates of the Keplerian orbital elements within a two-body problem perturbed by the tides in both partners. The formulae for these rates have appeared in the literature in various forms, at times with errors. We reconsider, from scratch, the derivation of these rates and arrive at the Lagrange-type equations which, in some details, differ from the corresponding equations obtained previously by Kaula (1964). We also write down detailed expressions for $da/dt$, $de/dt$ and $di/dt$, to order $e^4$. They differ from Kaula's expressions which contain a redundant factor of $M/(M+M^{\prime}),$ with $M$ and $M^{\prime}$ being the masses of the primary and the secondary. As Kaula was interested in the Earth-Moon system, this redundant factor was close to unity and was unimportant in his developments. This factor, however, must be reinstated when Kaula's theory is applied to a binary composed of partners of comparable masses. We have found that, while it is legitimate to simply sum the primary's and secondary's inputs in $da/dt$ or $de/dt$, this is not the case for $di/dt$. So our expression for $di/dt$ differs from that of Kaula in two regards. First, the contribution due to the dissipation in the secondary averages out when the apsidal precession is uniform. Second, we have obtained an additional term which emerges owing to the conservation of the angular momentum: a change in the inclination of the orbit causes a change of the primary's plane of equator.

Tables

Table 1. Table 1: Symbol key.

Variable	Explanation	Reference
$𝝆$	inertial position of the primary
${𝝆}^{\prime }$	inertial position of the secondary
$𝒓$	position of the secondary in the primary’s equatorial frame	eqn (5)
${𝒓}^{\prime }$	position of the primary in the secondary’s equatorial frame	eqn (5)
$U$	tidal potential of the deformed primary
$𝑭$	tidal force acting on the secondary
	due to the primary’s deformation	eqn (124)
$U^{\prime }$	tidal potential of the deformed secondary
${𝑭}^{\prime }$	tidal force acting on the primary
	due to the secondary’s deformation	eqn (125)
$M$	mass of the primary
$M^{\prime }$	mass of the secondary
$R$	radius of the primary
$R^{\prime }$	radius of the secondary
$\rho$	mean density of the primary
${\rho }^{\prime }$	mean density of the secondary
$\theta$	rotation angle of the primary
${\theta }^{\prime }$	rotation angle of the secondary
$a$	semimajor axis of the mutual orbit
$e$	eccentricity of the mutual orbit
$i$	orbit inclination on the primary’s equator
$i^{\prime }$	orbit inclination on the secondary’s equator
$\Omega$	longitude of the ascending node on the primary’s equator
${\Omega }^{\prime }$	longitude of the ascending node on the secondary’s equator
$\omega$	argument of the pericentre on the primary’s equator
${\omega }^{\prime }$	argument of the pericentre on the secondary’s equator
$ℳ$	mean anomaly
$n\equiv \dot{ℳ}$	anomalistic mean motion
$F_{lmp}\left(i\right)$	inclination functions	eqn (134)
$G_{lpq}\left(e\right)$	eccentricity function	eqn (134)
${\omega }_{lmpq}$	Fourier modes of the tides in the primary	eqn (140)
${\chi }_{lmpq}\equiv \left|{\omega }_{lmpq}\right|$	forcing frequencies excited in the primary
${ϵ}_l={ϵ}_l\left({\omega }_{lmpq}\right)$	tidal phase lags in the primary
$k_l=k_l\left({\omega }_{lmpq}\right)$	dynamical Love numbers of the primary
$K_l\left({\omega }_{lmpq}\right)\equiv k_l\left({\omega }_{lmpq}\right)\sin {ϵ}_l\left({\omega }_{lmpq}\right)$	quality functions of the primary
${\omega }_{lmpq}^{\prime }$	Fourier modes of the tides in the secondary	eqn (142)
${\chi }_{lmpq}^{\prime }\equiv \left|{\omega }_{lmpq}^{\prime }\right|$	forcing frequencies excited in the secondary
${ϵ}_l^{\prime }={ϵ}_l^{\prime }\left({\omega }_{lmpq}^{\prime }\right)$	tidal phase lags in the secondary
$k_l^{\prime }=k_l\left({\omega }_{lmpq}\right)$	dynamical Love numbers of the secondary
$K_l^{\prime }\left({\omega }_{lmpq}\right)\equiv k_l^{\prime }\left({\omega }_{lmpq}\right)\sin {ϵ}_l^{\prime }\left({\omega }_{lmpq}\right)$	quality functions of the primary
$𝒢$	gravitational constant

Equations

$${\bm{a}}\;=\;\ddot{\bf{\mbox{\boldmath${\boldmath\rho}$}}}^{\,\prime}\;-\;\ddot{\bf{\mbox{\boldmath${\boldmath\rho}$}}}\;=\;\frac{\mbox{\boldmath${\boldmath{F}}$}}{M^{\,\prime}}\;-\;\frac{(-\,\mbox{\boldmath${\boldmath{F}}$})}{M}\;=\;\frac{M\,+\,M^{\,\prime}}{M\;M^{\,\prime}}\;\mbox{\boldmath${\boldmath{F}}$}\,\;.$$

$$\beta =\frac{M{M}^{\prime }}{M+M^{\prime }}.$$

$$\text{varmathbb}R\equiv \text{varmathbb}{R}_3(\psi )\text{varmathbb}{R}_1(-\epsilon )\text{varmathbb}{R}_3(\theta ),$$

$$\text{varmathbb}{R}_1(\phi )=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \phi & -\sin \phi \\ 0& \sin \phi & \cos \phi \end{array}\right],\text{varmathbb}{R}_3(\phi )=\left[\begin{array}{ccc}\cos \phi & -\sin \phi & 0\\ \sin \phi & \cos \phi & 0\\ 0& 0& 1\end{array}\right].$$

$$\mathbf{r}={}^{\mathrm{T}}\text{varmathbb}R({\mathbf{\rho }}^{\prime }-\mathbf{\rho }),{\mathbf{r}}^{\prime }={}^{\mathrm{T}}\text{varmathbb}{R}^{\prime }({\mathbf{\rho }}^{\prime }-\mathbf{\rho }),$$

$$T=\frac{1}{2}\mathbf{\Omega }\cdot \text{varmathbb}{I}_{\mathrm{n}}\mathbf{\Omega }+\frac{1}{2}{\mathbf{\Omega }}^{\prime }\cdot \text{varmathbb}{I}_{\mathrm{n}}^{\prime }{\mathbf{\Omega }}^{\prime }+\frac{1}{2}\beta \Vert \dot{\mathbf{r}}+\mathbf{\Omega }\times \mathbf{r}{\Vert }^2,$$

$$V(\mathbf{r},{\mathbf{r}}^{\prime })=V_0(\mathbf{r})+V_1(\mathbf{r},{\mathbf{r}}^{\prime }).$$

$$\mathcal{F}=\mathcal{L}+\mathbf{\Lambda }\cdot (\text{varmathbb}R\mathbf{r}-{\text{varmathbb}R}^{\prime }{\mathbf{r}}^{\prime }).$$

$$\frac{d}{dt}f(\mathbf{\Theta })=\mathbf{\Omega }\cdot \hat{\mathbf{J}}(f)=\mathbf{\Omega }\cdot \text{varmathbb}J\frac{\partial f}{\partial \mathbf{\Theta }},$$

$$\frac{d}{dt}f(\mathbf{\Theta })=\dot{\mathbf{\Theta }}\cdot \frac{\partial f}{\partial \mathbf{\Theta }}.$$

$${\Omega }_X$$

$${\Omega }_Y$$

$${\Omega }_Z$$

$$\text{varmathbb}J=\left[\begin{array}{ccc}-\frac{\sin \theta }{\sin \epsilon }& -\cos \theta & \sin \theta \cot \epsilon \\ -\frac{\cos \theta }{\sin \epsilon }& \sin \theta & \cos \theta \cot \epsilon \\ 0& 0& 1\end{array}\right].$$

$$\frac{d}{dt}(\mathbf{\Lambda }\cdot \text{varmathbb}R\mathbf{r})=\mathbf{\Lambda }\cdot \dot{\text{varmathbb}R}\mathbf{r}=\mathbf{\Lambda }\cdot \text{varmathbb}R{}^{\mathrm{T}}\text{varmathbb}R\dot{\text{varmathbb}R}\mathbf{r}=({}^{\mathrm{T}}\text{varmathbb}R\mathbf{\Lambda })\cdot \left(\mathbf{\Omega }\times \mathbf{r}\right)=\mathbf{\Omega }\cdot \left(\mathbf{r}\times ({}^{\mathrm{T}}\text{varmathbb}R\mathbf{\Lambda })\right).$$

$$\hat{\mathbf{J}}(\mathbf{\Lambda }\cdot \text{varmathbb}R\mathbf{r})=\mathbf{r}\times ({}^{\mathrm{T}}\text{varmathbb}R\mathbf{\Lambda }).$$

$$\frac{d}{dt}\frac{\partial \mathcal{F}}{\partial \mathbf{\Omega }}$$

$$\frac{d}{dt}\frac{\partial \mathcal{F}}{\partial {\mathbf{\Omega }}^{\prime }}$$

$$\frac{d}{dt}\frac{\partial \mathcal{F}}{\partial \dot{\mathbf{r}}}$$

$$\frac{d}{dt}\frac{\partial \mathcal{F}}{\partial {\dot{\mathbf{r}}}^{\prime }}$$

$$\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \mathbf{\Omega }}$$

$$\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial {\mathbf{\Omega }}^{\prime }}$$

$$\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{\mathbf{r}}}$$

$$\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial {\dot{\mathbf{r}}}^{\prime }}$$

$$\mathbf{\Lambda }=\text{varmathbb}{R}^{\prime }\frac{\partial \mathcal{L}}{\partial {\mathbf{r}}^{\prime }}.$$

$$\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \mathbf{\Omega }}$$

$$\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial {\mathbf{\Omega }}^{\prime }}$$

$$\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{\mathbf{r}}}$$

$$\ddot{\mathbf{r}}+\frac{\mathcal{G}(M+M^{\prime })}{r^3}\mathbf{r}=-\dot{\mathbf{\Omega }}\times \mathbf{r}-2\mathbf{\Omega }\times \dot{\mathbf{r}}-\mathbf{\Omega }\times (\mathbf{\Omega }\times \mathbf{r})-\frac{1}{\beta }\left(\frac{\partial V_1}{\partial \mathbf{r}}+{}^{\mathrm{T}}\text{varmathbb}R\text{varmathbb}{R}^{\prime }\frac{\partial V_1}{\partial {\mathbf{r}}^{\prime }}\right).$$

$$\ddot{\mathbf{r}}+\frac{\mathcal{G}(M+M^{\prime })}{r^3}\mathbf{r}=-\dot{\mathbf{\Omega }}\times \mathbf{r}-2\mathbf{\Omega }\times \dot{\mathbf{r}}-\mathbf{\Omega }\times (\mathbf{\Omega }\times \mathbf{r})+\left(\frac{\partial \mathcal{R}}{\partial \mathbf{r}}+{}^{\mathrm{T}}\text{varmathbb}R\text{varmathbb}{R}^{\prime }\frac{\partial \mathcal{R}}{\partial {\mathbf{r}}^{\prime }}\right).$$

$$\mathbf{\Sigma }$$

$${\mathbf{\Pi }}^{\prime }$$

$$\mathbf{p}$$

$$T=\frac{1}{2}(\mathbf{\Sigma }-\mathbf{r}\times \mathbf{p})\cdot \text{varmathbb}{I}_{\mathrm{n}}^{-1}(\mathbf{\Sigma }-\mathbf{r}\times \mathbf{p})+\frac{1}{2}{\mathbf{\Pi }}^{\prime }\cdot \text{varmathbb}{I}_{\mathrm{n}}^{\prime -1}{\mathbf{\Pi }}^{\prime }+\frac{\Vert \mathbf{p}{\Vert }^2}{2\beta }.$$

$$\frac{d\mathbf{\Sigma }}{dt}$$

$$\frac{d{\mathbf{\Pi }}^{\prime }}{dt}$$

$$\frac{d\mathbf{p}}{dt}$$

Full text

We also write down detailed expressions for $\,da/dt\,$, $\,de/dt\,$ and $\,di/dt\,$, to order $\,e^{4}\,$. They differ from Kaula’s expressions which contain a redundant factor of $\,M/(M+M^{\,\prime})\,$, with $\,M\,$ and $\,M^{\,\prime}\,$ being the masses of the primary and the secondary. As Kaula was interested in the Earth-Moon system, this redundant factor was close to unity and was unimportant in his developments. This factor, however, must be removed when Kaula’s theory is applied to a binary composed of partners of comparable masses.

We have found that, while it is legitimate to simply sum the primary’s and secondary’s inputs in $\,da/dt\,$ or $\,de/dt\,$, this is not the case for $\,di/dt\,$. So our expression for $\,di/dt\,$ differs from that of Kaula in two regards. First, the contribution due to the dissipation in the secondary averages out when the apsidal precession is uniform. Second, we have obtained an additional term which emerges owing to the conservation of the angular momentum: a change in the inclination of the orbit causes a change of the primary’s plane of equator.

1 Motivation

The Darwin-Kaula theory of bodily tides is a fundamental development with many ramifications. It provides the means for calculating spin-orbit evolution of planets and moons, including their entrapment in spin-orbit resonances (e.g, Correia, Laskar & de Surgy 2003; Correia & Laskar 2003; Makarov, Berghea & Efroimsky 2012; Noyelles et al. 2014) and the final obliquities (Cunha et al. 2015; Ferraz-Mello et al. 2008). In the situations where tidal heating is intensive, an approach based on this theory gives the key to the thermal histories of celestial bodies (e.g., Peale & Cassen 1978; Efroimsky & Makarov 2014; Efroimsky 2018; Makarov et al. 2018). This theory also enables one to calculate the influence of the lunisolar tides on the orbital motion of an artificial spacecraft (Pucacco & Lucchesi 2018).

In the current paper, we address the orbit evolution within a two-body problem perturbed by the tides in both partners. Specifically, we are interested in the tidal rates $\,da/dt\,$, $\,de/dt\,$, and $\,di/dt\,$ to order $\,e^{4}\,$ (the symbols being used as defined in Table 1). Formulae for these rates appeared in the literature in various forms, usually to a lower precision and sometimes with mistakes. So we compare our expressions with those suggested in some other publications, including the cornerstone work by Kaula (1964).

For the additional tidal potential of a disturbed body, Kaula (1964) developed an expansion valid for an arbitrary rheology (i.e., for an arbitrary frequency-dependence of the quality function $\,k_{l}/Q_{l}\,$). Kaula’s derivation was terse and omitted several steps as self-evident. We accurately fill in these gaps, and point out a step at which Kaula made a tacit approximation $\,M\gg M^{\,\prime}\,$, where $\,M\,$ and $\,M^{\,\prime}\,$ are the masses of the primary and the secondary. Owing to this assumption, Kaula’s expressions for the orbital elements’ rates contain a redundant factor of $\,M/(M+M^{\,\prime})\,$. Since Kaula was concerned with the Earth-Moon system, this approximation made little difference as the factor was close to unity. However, in the case of a binary composed of bodies of comparable masses, this redundant factor must be removed from these expressions.

We also reexamine from scratch Kaula’s derivation of the rates $\,da/dt\,$, $\,de/dt\,$, and $\,di/dt\,$. We find that, while it is legitimate to simply sum the primary’s and secondary’s inputs in $\,da/dt\,$ or $\,de/dt\,$, this is not the case for $\,di/dt\,$. It turns out that in the expression for the primary’s $\,di/dt\,$ the contribution due to the dissipation in the secondary averages out when the apsidal precession is uniform. Also, in that expression we obtain an additional term emerging from the conservation of the angular momentum: a change in the inclination of the orbit causes a change of the primary’s plane of equator. For these two reasons, our formula for $\,di/dt\,$ differs considerably from that of Kaula.

2 Basics

2.1 The two-body problem perturbed by tides

Consider two near-spherical bodies. One, called “planet” or “primary”, has a mass $\,M\,$ and an inertial position  $\rho$ . Another, named “secondary”, has a mass $\,M^{\,\prime}\,$ and is residing in \,\mbox{{\boldmath\rho}}^{\,\prime}\,. We are interested in the orbital evolution of this system, with tides in both partners taken into account. Within this setting, Kaula (1964) expressed the perturbing gravitational potential through the Keplerian elements of the mutual orbit, thus allowing him to describe the evolution of the system by Lagrange’s planetary equations. Three caveats are in order regarding Kaula’s development.

First of all, by definition of a mutual orbit, the position of the secondary is measured with respect to the primary’s centre of mass. Let  ${\boldmath{F}}$  be the force exerted by the primary on the secondary. By virtue of Newton’s third law of motion, the secondary simultaneously exerts a force \,-\,\mbox{\boldmath{\boldmath{F}}}\, on the primary. Hence, the mutual acceleration reads as

[TABLE]

In the limit of the secondary being a test particle ($\,M\,\gg\,M^{\,\prime}\,$), the above relation becomes simply \,\bm{a}\,=\,\mbox{\boldmath{\boldmath{F}}}/M^{\,\prime}\,, which was the tacit approximation accepted by Kaula. We in our study shall rely on the general expression (1) and therefore shall employ the reduced mass $\,\beta\,$ defined as

[TABLE]

Secondly, it should be noted that within Kaula’s formalism the elements $\,(a,\,e,\,i,\,{\cal M},\,\omega,\,\Omega)\,$ of the mutual orbit are defined in an inertially fixed frame coinciding with the primary’s equator at the instant when the equations of motion are computed. To simplify the interpretation of the Keplerian elements, we here assume that they are defined in a frame coprecessing with the primary’s equator. 111 In this frame, the role of the origin of longitude is played by the descending node of the primary’s equator on an inertial plane. The perturbing forces showing up in this setting include the inertial forces associated with the non-Galilean nature of the coprecessing frame. An important feature of these forces is that they depend not only on the positions but also on velocities. As was pointed out in Efroimsky (2005a,b), for such kind of disturbances the Lagrange- and Delaunay-type planetary equations in their standard form render orbital elements which are not osculating. To be more precise, if we wish the orbital element to osculate the orbit defined in a noninertial frame, not only must we amend the disturbing function, but we should also insert in the Lagrange- and Delaunay-type equations additional terms that are not part of the disturbing function. If, however, we choose only to amend the disturbing function, the orbital equations will give us Keplerian orbital elements which will not osculate with the orbit defined in the noninertial frame, but will instead osculate with the orbit as seen in the inertial frame. 222 Such orbital elements are sometimes called  contact elements, in order to distinguish them from their osculating counterparts. For example, the semimajor axis $\,a_{osc}\,$ and eccentricity $\,e_{osc}\,$ osculating in the noninertial frame wherein the orbit is defined are linked via \,\beta\,(\mbox{{\boldmathr}}\times\dot{\textbf{\mbox{\boldmath{\boldmath r}}}})\,=\,\beta\,\sqrt{{\cal G}(M+M^{\,\prime})\,a_{osc}\,(1-e_{osc}^{2})}\,\hat{\bf{w}}_{osc}\, to the relative position  $r$  and the velocity \,\dot{\textbf{\mbox{\boldmath{\boldmath r}}}}\, in that same noninertial frame. At the same time, the contact elements $\,a\,$ and $\,e\,$ (those rendered by the Lagrange- or Delaunay-type equations with only the disturbing function amended) are connected through \,\mbox{{\boldmathr}}\times\mbox{{\boldmathp}}\,=\,\beta\,\sqrt{{\cal G}(M+M^{\,\prime})\,a\,(1-e^{2})}\,\hat{\bf{w}}\, with the relative position  $r$  and the momentum  $p$  which is equal to the reduced mass multiplied by the velocity in the inertial frame: \,\mbox{{\boldmathp}}\,=\,\beta\,(\dot{\textbf{\mbox{\boldmath{\boldmath r}}}}\,+\,\mbox{{\boldmath\bm{\mu}}}\times\mbox{{\boldmathr}})\;. In these formulae,  $\bm{\mu}$  is the precession rate of the noninertial frame relative the inertial one, while $\,\hat{\bf{w}}_{osc}\,$ and $\,\hat{\bf{w}}\,$ denote the unit vector in the direction of the orbital angular momentum, as seen in the noninertial and inertial frames, correspondingly. For comprehensive treatment, see Efroimsky (2005a,b) and references therein.

While the treatment in the coprecessing frame has its advantages, there exists alternatives to it: both the orbital motion and the primary’s spin can be described in the Laplace plane (e.g., Boué, Correia & Laskar 2016; Rubincam 2016), or the primary’s spin can be reckoned from the orbital plane whose Keplerian elements are given in the Laplace plane (e.g., Néron de Surgy & Laskar 1997; Correia, Laskar & Néron de Surgy 2003; Correia & Laskar 2003; Correia & Laskar 2010).

Lastly, a complete description of motion requires not only the aforementioned set of the Keplerian elements relative to the primary’s equator, but also a set of the elements $\,(a^{\,\prime},\,e^{\,\prime},\,i^{\,\prime},\,{\cal M}^{\,\prime},\,\omega^{\,\prime},\,\Omega^{\,\prime})\,$ relative to the secondary’s equator. 333 As we shall see below, in the paragraph after equation (42), the instantaneous ellipses, as seen from the primary’s and secondary’s equator, will have the same shape, with $\,a=a^{\,\prime}\,$ and $\,e=e^{\,\prime}\,$ (and also with $\,{\cal{M}}={\cal{M^{\,\prime}}}\,$). Generally, though, these equalities are not obligatory — see an example in Footnote 2 above. Therefore, the equations of motion must satisfy the relations existing between these two sets.

With all these details taken into account, we now derive the Lagrange-type planetary equations compatible with Kaula’s (1964) formalism.

2.2 Lagrangian formalism

2.2.1 Lagrangian function

Let $\varmathbb{I}_{\mathrm{n}}$ be the primary’s inertia matrix and $\,{\varmathbb{R}}\,$ the rotation matrix, function of the Euler (3,-1,3) angles $\,\Theta=(\psi,\varepsilon,\theta)\,$, describing orientation in an inertial frame. More explicitly,

[TABLE]

where $\varmathbb{R}_{1}$ and $\varmathbb{R}_{3}$ represent the rotation matrices around the first and third axis, respectively:

[TABLE]

The angle $\psi$ is thus the longitude of the descending node of the equator with respect to the inertial frame, $\varepsilon$ the inclination of the spin axis in the same inertial frame, and $\theta$ the rotation angle around the body’s figure axis (see Figure 1). This convention is chosen for $\theta$ to be reckoned from the descending node of the equator, as in Kaula (1964). We similarly denote by $\,\varmathbb{I}_{\mathrm{n}}^{\,\prime}\,$ and $\,{\varmathbb{R}}^{\,\prime}(\,\bm{\Theta}^{\,\prime}=(\psi^{\,\prime},\varepsilon^{\,\prime},\theta^{\,\prime})\,)\,$ the inertia and rotation matrices of the secondary, respectively. The planetocentric position of the secondary will be denoted with $\bm{r}$, when expressed in the body-fixed frame of the primary, or with $\bm{r}^{\,\prime}$, when expressed in the body-fixed frame of the secondary :

[TABLE]

where ${}^{\mathrm{T}}{(\cdot)}$ denotes the transposition operator. Notice that we do not strictly follow Kaula’s convention because, to represent the orbit, we are employing the corotating frames instead of the coprecessing ones. At the end of the derivation, we explain how to switch between these two classes of frames.

For the state of the system to be entirely defined, we also introduce the primary’s and secondary’s angular velocities $\,\bm{\Omega}\,$ and $\,\bm{\Omega}^{\,\prime}\,$, respectively. These vectors are expressed in their respective body-fixed frame.

The Lagrangian $\,{\cal L}=T-V\,$ of the system is a function of $\,(\bm{\Omega},\,\bm{\Omega}^{\,\prime},\,\dot{\bm{r}\,},\,\bm{r},\,\bm{r}^{\,\prime})\,$. Specifically, the kinetic energy $\,T(\bm{\Omega},\,\bm{\Omega}^{\,\prime},\,\dot{\bm{r}},\,\bm{r})\,$ is given by

[TABLE]

where $\,\beta=MM^{\,\prime}/(M+M^{\,\prime})\,$ is the reduced mass of the system. As for the potential energy $\,V(\bm{r},\bm{r}^{\,\prime})\,$, we decompose it into the point mass potential energy $\,V_{0}(\bm{r})=-{\cal G}MM^{\,\prime}/r\,$ and a perturbation $V_{1}(\bm{r},\bm{r}^{\,\prime})$, whose expression will be specified later on:

[TABLE]

To this Lagrangian function, we have to add the constraint $\,{\varmathbb{R}}\bm{r}={\varmathbb{R}}^{\,\prime}\bm{r}^{\,\prime}\,$ which links two expressions of the same quantity $\,\bm{\rho}^{\,\prime}-\bm{\rho}\,$. This constraint will enter the Lagrangian, accompanied with Lagrange multipliers $\,\bm{\Lambda}\in\mathbb{R}^{3}\,$, leading to a new Lagrangian

[TABLE]

Note that this Lagrangian also depends on the additional variables $\,(\bm{\Theta},\,\bm{\Theta}^{\,\prime})\,$.

2.2.2 Spin operator

To derive the Euler-Poincaré-Lagrange equations of motion for this Lagrangian, let us first introduce the spin operator $\hat{\bm{J}}\,\equiv\,\varmathbb{J}\,\partial/\partial\bm{\Theta}\,$, where $\,\varmathbb{J}\,$ is a matrix yet to be defined (e.g., Boué 2017; Boué et al. 2017). On the one hand, by definition of the spin operator, 444 For a detailed introduction in the theory of the spin operator, see Varshalovich et al. (1988). the time derivative of an arbitrary function $\,f(\bm{\Theta})\,$ can be written as

[TABLE]

with $\bm{\Omega}$ the angular velocity expressed in the same frame as $\,\hat{\bm{J}}\,$, which here is the body-fixed frame. On the other, applying the chain rule for the time derivative of $\,f(\bm{\Theta})\,$, we get

[TABLE]

By identification of (9) and (10), we deduce that $\varmathbb{J}$ is the matrix such that $\dot{\bm{\Theta}}={}^{\mathrm{T}}{\varmathbb{J}}\,\bm{\Omega}$, or equivalently, $\bm{\Omega}={}^{\mathrm{T}}{\varmathbb{J}}^{-1}\dot{\bm{\Theta}}$.

Knowing that the components of $\bm{\Omega}$ in the body-fixed frame (rotated with respect to an inertial frame according to (3)) are

[TABLE]

we get

[TABLE]

The matrix $\varmathbb{J}^{\,\prime}$ is equivalently defined for the secondary.

In the following, we also have to determine the image of the function $\varmathbb{R}\mapsto\bm{\Lambda}\cdot\varmathbb{R}\bm{r}$ by the spin operator $\hat{\bm{J}}$ which, by definition, evaluates the variation of a function under infinitesimal rotation of the primary (therefore the rotation matrix $\varmathbb{R}$ alone is affected by the operator $\hat{\bm{J}}$ and $\bm{r}$ shall be taken constant in this calculation). Using ${}^{\mathrm{T}}{\varmathbb{R}}\dot{\varmathbb{R}}\,\bm{v}=\bm{\Omega}\times\bm{v}$ for any vector $\bm{v}\in\mathbb{R}^{3}$, we get

[TABLE]

Identifying (9) with (15), we obtain

[TABLE]

2.2.3 Equations of motion

Using the matrices $\varmathbb{J}$ and $\varmathbb{J}^{\,\prime}$ defined hereinabove, the Euler-Poincaré-Lagrange equations of motion read (e.g., Boué 2017; Boué et al. 2017)

[TABLE]

Notice that the Euler-Lagrange equation (20) is equal to zero. This is due to the fact that the Lagrangian does not depend on $\dot{\bm{r}}^{\,\prime}$. Let us now rewrite the equations of motion (17-20) in terms of the original Lagrangian $\cal L$ :

[TABLE]

with, for the problem studied in this paper, $\partial{\cal L}/\partial\bm{\Theta}=\partial{\cal L}/\partial\bm{\Theta}^{\,\prime}=\bm{0}$. To derive the first two equations, we made use of relation (16).

From the equation of motion (24), we determine the expression of the Lagrange multiplier, namely,

[TABLE]

Substituting this expression in the other equations of motion (21-23), we get

[TABLE]

In the first equation of motion, we made use of the relations $\bm{r}={}^{\mathrm{T}}{\varmathbb{R}}\varmathbb{R}^{\,\prime}\bm{r}^{\,\prime}$ and $(\varmathbb{R}\bm{u}\times\varmathbb{R}\bm{v})=\varmathbb{R}(\bm{u}\times\bm{v})$. In the last equation of motion, we recognise the force $\bm{F}^{\,\prime}\equiv\partial{\cal L}/\partial\bm{r}^{\,\prime}$ written in the body-fixed frame of the primary thanks to the rotation matrix ${}^{\mathrm{T}}{\varmathbb{R}}\varmathbb{R}^{\,\prime}$. Moreover, in the first two equations of motion, we observe the presence of the torque $\bm{r}^{\,\prime}\times\bm{F}^{\,\prime}$ expressed in the body-fixed frame of the primary in equation (26) and in the body-fixed frame of the secondary in equation (27).

Inserting the expressions for the kinetic energy (6) and the potential energy (7) in equation (28), we get

[TABLE]

Were the right-hand side equal to zero, we would have obtained the equation of motion of the classical two-body problem. But here this is not the case. The first three terms of the right-hand side account for the inertial forces of the non-Galilean frame in which $\bm{r}$ is expressed, while the last two terms represent the perturbation induced by tides. It is common to name the quantity $\,-\,V_{1}/\beta\,$ as the perturbing function and to denote it with $\,\cal R\,$. In this notation, equation (29) reads:

[TABLE]

2.3 Hamiltonian formalism

To get the Hamiltonian form of the equations of motion, we apply a Legendre transformation on the Lagrangian. Let $\bm{\Sigma}$, $\bm{\Pi}^{\,\prime}$ and $\bm{p}$ be the generalised momenta given by

[TABLE]

$\bm{\Sigma}=\bm{\Pi}+\bm{\Gamma}$ is the sum of the angular momenta of the primary ($\bm{\Pi}$) and of the orbit ($\bm{\Gamma}=\bm{r}\times\bm{p}$), $\bm{\Pi}^{\,\prime}$ is the angular momentum of the secondary and $\bm{p}$ is the linear orbital momentum with respect to the inertial frame but expressed in the body-fixed frame of the primary.

The Hamiltonian ${\cal H}\equiv\bm{\Sigma}\cdot\bm{\Omega}+\bm{\Pi}^{\,\prime}\cdot\bm{\Omega}^{\,\prime}+\bm{p}\cdot\dot{\bm{r}}-{\cal L}$ can be written as ${\cal H}=T+V$ with

[TABLE]

The equations of motion deduced from the Legendre transformation are

[TABLE]

The kinematic equations of motion are

[TABLE]

Recall that these expressions are general. In our case, the two-body problem perturbed by tides, $\,{\cal H}\,$ is independent of $\,\bm{\Theta}\,$ and $\,\bm{\Theta}^{\,\prime}\,$. So the partial derivatives of the Hamiltonian over these two quantities in equations (35) and (36) become zero.

2.4 Elliptical elements

In Kaula’s work, the equations of motion are written in terms of elliptical elements ($a$, $e$, $i$, ${\cal{M}}$, $\omega$, $\Omega$) reckoned from an inertial frame instantaneously comoving with the primary’s precessing equator. 555 The frame employed by Kaula (1964) should be termed “instantaneously comoving”, not coprecessing, because in his equations of motion the inertial forces were omitted. This set of variables becomes singular at zero inclination or zero eccentricity. We will nevertheless provide the equations of motion in this set of variables for an easier comparison with previous works. Here, we define elliptical elements $(a,e,i,{\cal M},\omega,\bar{\Omega})$ (represented in Figure 2) as a change of variable from the conjugated variables $(\bm{p},\bm{r})$. Therefore, they describe the instantaneous ellipse $\cal E$ constructed from the position vector $\bm{r}$ and its inertial velocity $\dot{\bm{r}}+\bm{\Omega}\times\bm{r}$, both defined in the primary body-fixed frame. Our set of Keplerian elements differs from Kaula’s by the frame in which it is defined. This choice only impacts the longitude of the ascending node, hence the introduction of a new symbol $\bar{\Omega}$. The two longitudes of ascending node are related to each other by

[TABLE]

As for the couple $(\bm{p},\bm{r})$, to the position vector $\bm{r}^{\,\prime}$ we associate the elliptical elements $(a^{\,\prime},e^{\,\prime},i^{\,\prime},{\cal M}^{\,\prime},\omega^{\,\prime},$ $\bar{\Omega}^{\,\prime})$ of the ellipse $\cal E^{\,\prime}$ defined by $\bm{r}^{\,\prime}$ and its inertial velocity, both expressed in the secondary body-fixed frame. Thus, both ellipses $\cal E$ and $\cal E^{\,\prime}$ are the same up to a rotation, i.e., $a=a^{\,\prime}$, $e=e^{\,\prime}$ and ${\cal M}={\cal M}^{\,\prime}$. To get the relation between $(\bar{\Omega}^{\,\prime},i^{\,\prime},\omega^{\,\prime})$ and $(\bar{\Omega},i,\omega)$, we consider the following function with values in $\mathrm{SO}(3)$ :

[TABLE]

We have $\varmathbb{F}={}^{\mathrm{T}}{\varmathbb{R}}\varmathbb{R}^{\,\prime}$. Hence, for fixed orientation of the bodies (as in equations (37) and (40)), we get $d\varmathbb{F}=\bm{0}$ and therefore $d\varmathbb{F}\,{}^{\mathrm{T}}{\varmathbb{F}}=\bm{0}$. Let us denote by $\bm{K}$, $\bm{i}$, $\bm{k}$, $\bm{i}^{\,\prime}$, $\bm{K}^{\,\prime}$ the unit vectors of the rotations of angle $\bar{\Omega}$, $i$, $\omega-\omega^{\,\prime}$, $-i^{\,\prime}$, and $-\bar{\Omega}^{\,\prime}$, respectively. $\bm{K}$ and $\bm{K}^{\,\prime}$ are the primary’s and secondary’s figure axes, respectively, $\bm{i}$ and $\bm{i}^{\,\prime}$ are the directions of the orbit ascending node relative to the primary’s and the secondary’s equatorial plane, respectively, and $\bm{k}$ is the orbit normal (see Figure 2). In the orbital reference frame, we have in particular

[TABLE]

The product $d\varmathbb{F}\,{}^{\mathrm{T}}{\varmathbb{F}}$ belongs to the Lie algebra $\mathrm{so}(3)$. Let $d\varmathbb{X}\equiv\varmathbb{R}_{1}(-i)\varmathbb{R}_{3}(-\bar{\Omega})\,d\varmathbb{F}\,{}^{\mathrm{T}}{\varmathbb{F}}\,\varmathbb{R}_{3}(\bar{\Omega})\varmathbb{R}_{1}(i)$ be its expression in the orbit frame. We have (see Appendix A)

[TABLE]

where the hat over any vector $\bm{v}=(v_{x},v_{y},v_{z})$ denotes the skew-symmetric matrix666The skew-symmetric matrix is so defined that for any two vectors $\bm{a},\bm{b}\in\mathbb{R}^{3}$, their vector product $\bm{a}\times\bm{b}$ is equal to $\hat{\bm{a}}\bm{b}$.

[TABLE]

Applying to $\,d\varmathbb{X}=\bm{0}\,$ the canonical bijection from $\mathrm{so}(3)$ to $\mathbb{R}^{3}$ (i.e., the inverse of the hat application), we get

[TABLE]

We now replace the vectors $\bm{K}$, $\bm{i}$, $\bm{k}$, $\bm{i}^{\,\prime}$, $\bm{K}^{\,\prime}$ by their coordinates and get

[TABLE]

We finally deduce the Jacobian $\varmathbb{A}_{1}$ of the transformation $(\bar{\Omega}^{\,\prime},i^{\,\prime},\omega^{\,\prime})\rightarrow(\bar{\Omega},i,\omega)$, which reads

[TABLE]

and thus,

[TABLE]

2.5 Equations of motion of the Keplerian elements

Planetary equations of motion for $(a,e,i,{\cal M},\omega,\bar{\Omega})$ are deduced from the canonical equations of motion satisfied by Delaunay variables (rescaled by the reduced mass)

[TABLE]

Let $\bm{X}={}^{\mathrm{T}}{(L,G,H)}$, $\bm{x}={}^{\mathrm{T}}{(l,g,h)}$, $\bm{Y}={}^{\mathrm{T}}{(a,e,i)}$ and $\bm{y}={}^{\mathrm{T}}{({\cal M},\omega,\bar{\Omega})}$. We have

[TABLE]

thus

[TABLE]

where $\varmathbb{M}$ is the Jacobian defined as

[TABLE]

The result, written in a matrix form, reads

[TABLE]

These are the classical planetary equations in the form of Lagrange. Let us denote by $\,\varmathbb{B}\,$ the Poisson matrix, i.e., the matrix standing before the gradient of the Hamiltonian in equation (55). In the problem under consideration, an adjustment to this equation has to be made. Owing to the constraint between $\bm{r}$ and $\bm{r}^{\,\prime}$, we have to add the contribution from $\,(\bar{\Omega}^{\,\prime},i^{\,\prime},\omega^{\,\prime})\,$ to the time derivative of the state vector, see equation (37). This is done through the medium of the Jacobian $\,\varmathbb{A}\,$ as follows (see Appendix B)

[TABLE]

where $\,\varmathbb{B}^{\,\prime}\,$ is the equivalent of the matrix $\,\varmathbb{B}\,$ but written as a function of $\,(a,\,e,\,i^{\,\prime})\,$ instead of $\,(a,\,e,\,i)\,$. After some algebra, we arrive at

[TABLE]

2.6 Perturbed two-body problem

Here, we split the Hamiltonian as ${\cal H}={\cal H}_{0}+V_{1}$ with ${\cal H}_{0}=T+V_{0}$, i.e.,

[TABLE]

We shall now write this Hamiltonian in terms of the elliptical elements $(a,e,i,{\cal M},\omega,\bar{\Omega})$. First, we recognise the Keplerian energy of the two-body problem

[TABLE]

Then, $\bm{r}\times\bm{p}$ is the orbital angular momentum, thus

[TABLE]

Moreover, we recall that $\varmathbb{I}_{\mathrm{n}}^{-1}(\bm{\Sigma}-\bm{r}\times\bm{p})$ is equal to the primary’s angular velocity $\bm{\Omega}={}^{\mathrm{T}}{(\Omega_{X},\Omega_{Y},\Omega_{Z})}$. The gradient of the Hamiltonian ${\cal H}_{0}$ is then

[TABLE]

Let $n=({\cal G}(M+M^{\,\prime})/a^{3})^{1/2}$ be the Keplerian mean motion. The equations of motion become

[TABLE]

Using the expressions $L=\beta na^{2}$ and $G=L(1-e^{2})^{1/2}$ and substituting $V_{1}$ by $-\beta\cal R$, these equations of motion read

[TABLE]

Had we defined the Keplerian elements in the coprecessing frame of the primary rather than in its corotating frame, the longitude of the ascending node would have been $\Omega=\bar{\Omega}+\theta$. Then, in all the equations of motion, $\bar{\Omega}$ would have been replaced with $\Omega-\theta$, and in particular the last equation of motion would have become

[TABLE]

Let us now compare equations (78-83) with the well-known Lagrange planetary equations. As the Keplerian elements are here defined with respect to a moving frame, the resulting equations of motion for $i$, $\omega$ and $\bar{\Omega}$ contain driving terms which are functions of the components $(\Omega_{X},\Omega_{Y},\Omega_{Z})$ of the angular velocity $\bm{\Omega}$ of this frame. Furthermore, we observe that the equations of motion are not symmetric in $(i,\omega,\bar{\Omega})$ and in $(i^{\,\prime},\omega^{\,\prime},\bar{\Omega}^{\,\prime})$. This is due to the rotation matrix ${}^{\mathrm{T}}{\varmathbb{R}}\varmathbb{R}^{\,\prime}$ between the frames in which the two sets of Keplerian elements are defined. We nevertheless recover the lost symmetry in the equations of motion when this matrix ${}^{\mathrm{T}}{\varmathbb{R}}\varmathbb{R}^{\,\prime}$ becomes a single rotation around the third axis (i.e., when $i=i^{\,\prime}$ and $\omega=\omega^{\,\prime}$).

2.7 Rotation equations of motion

Equations (78-83) are the equivalent of the Lagrange planetary equations. For completeness, we now provide the explicit equations of motion for the rotation of the two bodies, given by equations (35-36) and (38-39). While these equations should involve the spin operators represented by the matrices $\,\varmathbb{J}\,$ and $\,\varmathbb{J}^{\,\prime}\,$, we discard those terms as the perturbing potential energy is independent of $\,\bm{\Theta}\,$ and $\,\bm{\Theta}^{\,\prime}\,$. The equations also involve the orbital angular momentum operator $\,\hat{\bm{L}}^{\,\prime}=\bm{r}^{\,\prime}\times\partial_{\bm{r}^{\,\prime}}\,$ expressed by a matrix $\,\varmathbb{L}^{\,\prime}\,$ such that

[TABLE]

for all functions $\,f(a^{\,\prime},\,e^{\,\prime},\,i^{\,\prime},\,{\cal M}^{\,\prime},\,\omega^{\,\prime},\,\bar{\Omega}^{\,\prime})\,$. Applying the same approach as in Section 2.2, we arrive at

[TABLE]

Developing the matrix products, we obtain the following dynamical equations of motion:

[TABLE]

and the associated kinematic equations of motion:

[TABLE]

Recall that in these equations, $\,\bm{\Omega}=\varmathbb{I}_{\mathrm{n}}^{-1}(\bm{\Sigma}-\bm{r}\times\bm{p})\,$ and $\,\bm{\Omega}^{\,\prime}={\varmathbb{I}_{\mathrm{n}}^{\,\prime}}^{\,-1}\bm{\Pi}^{\,\prime}\,$.

2.8 The gyroscopic approximation and the approximation of constant inertia matrices

In the equations of motion derived above, the inertial forces are expressed in terms of the components of the rotation vectors $\,\bm{\Omega}\,$ and $\,\bm{\Omega}^{\,\prime}\,$. To simplify these dependencies, we carry out two steps.

• (a) 

We employ the gyroscopic approximation, i.e., assume that the rotation about the axis of maximal inertia is much faster than any change in this axis’ orientation. Mathematically, this implies:

[TABLE]

• (b) 

We neglect the variation of the inertia matrices in the expressions for the angular momenta $\,\bm{\Pi}\,$ and $\,\bm{\Pi}^{\,\prime}\,$. This approximation is nontrivial because, after all, the tidal theory is about deformation. So one always needs to justify accurately why on one occasion the deformations must be taken into account and neglected on other. This justification is provided in Frouard & Efroimsky (2017b, Section 3.1 and Footnote 6).

Let us denote the partners’ principal moments of inertia with $\,(A,\,B,\,C)\,$ and $\,(A^{\,\prime},\,B^{\,\prime},\,C^{\,\prime})\,$. Then the angular momentum of the primary can be written as

[TABLE]

In the afore-explained approximation, its rate becomes

[TABLE]

Similar formulae will be valid for the secondary.

The rate of the angular momentum,

[TABLE]

will, owing to equations (65), (72-77), and (87-88), be equal to

[TABLE]

Naturally, the evolution rate of $\,\bm{\Pi}\,$ has a form analogous to that of $\,\bm{\Pi}^{\,\prime}\,$, equations (90-92). Inserting the exact expressions (102 - 102) and the approximate expressions (103 - 105) in equations (107 - 109), we exclude the components of the angular momentum, to be left with the components of the angular velocity only:

[TABLE]

with $\Omega_{Z}\approx\dot{\theta}$. Similarly, for the secondary we have

[TABLE]

We can now rewrite the equations of motion (78-83) where $\Omega_{X}$ and $\Omega_{Z}$ are substituted by their expressions (110,111). Furthermore, we make the approximations $C\gg(B-A)$ and $\Omega_{Z}\approx\dot{\theta}$ and we replace $\bar{\Omega}$ with $(\Omega-\theta)$. We get

[TABLE]

To close the system, we have to add the equations

[TABLE]

Be mindful that in equation (121) we employed the relation $\,\Omega_{Z}-\dot{\theta}=\dot{\psi}\cos\varepsilon=-\,(\Omega_{X}\sin\theta+\Omega_{Y}\cos\theta)\cot\varepsilon\,$, while in equation  (122) the relation $\dot{\varepsilon}=-\,\Omega_{X}\cos\theta+\Omega_{y}\sin\theta$ was used. Besides, we would emphasise that equations (116-123) have been obtained under the gyroscopic approximation, for which reason it would be illegitimate to consider the limit of $\,\dot{\theta}\,\longrightarrow\,0\,$. Therefore, the presence of $\,\dot{\theta}\,$ in the denominator in equations (118) and (120 - 121) produce no singularity. In the cases where the gyroscopic approximation is invalid, one has to rely on the equations of motion (78-83) instead.

2.9 Comparison with Kaula (1964)

In his 1964 paper, Kaula was mainly interested in the evolution of the semimajor axis, the eccentricity and the inclination. The derivation of the equations (38) in Kaula (1964) should thus only be compared with our formulae (116), (117) and (118). Although the evolution rates of $a$ and $e$ are identical in both approaches up to the definition of the disturbing function (see Section 3.2), the two expressions of the time derivative of the inclination $i$ display important dissimilarities. This difference of behaviour between $(a,e)$ on the one hand and $i$ on the other is due to the fact that only $i$ is affected by a rotation of the reference frame. Let us recall that in Kaula (1964) elliptical elements are defined with respect to a fixed reference frame coinciding with the primary’s equator at the time when the equations of motion are evaluated, whereas our set of equations (116-121) is written in the coprecessing frame of the primary. Therefore, equation (118) contains an inertial force leading to a term in $1/(C\dot{\theta})$ which is absent in Kaula (1964, eqn 38). But these two equations have also distinct dependencies with respect to the primed Keplerian elements. In fact, Kaula erroneously assumed that the planetary and satellite contributions in $di/dt$ could simply be summed up. Within our formalism, this naive (and wrong) assumption is equivalent to omission of the rotation ${}^{\mathrm{T}}{\varmathbb{R}}\varmathbb{R}^{\,\prime}$ in equations (29) and (30). Under this omission, the components of both forces ($-\partial V_{1}/\partial\bm{r}$ and $-\partial V_{1}/\partial\bm{r}^{\,\prime}$) would have been illegitimately summed up, ignoring the fact that they were written in different coordinate systems: by construction, $\partial V_{1}/\partial\bm{r}$ is expressed in the primary’s frame, while $\partial V_{1}/\partial\bm{r}^{\,\prime}$ is expressed in the secondary’s.

3 Tidal potential energy

3.1 General expression

Within the Darwin-Kaula theory (Kaula 1964, Efroimsky 2012, Efroimsky & Makarov 2013), it is taken into account that in the general case the secondary body “feels” the tides which may be generated in the primary not only by the secondary itself but also by some other perturber located in \,\mbox{{\boldmathr}}^{\,*}\,. Then in an arbitrary exterior point  $r$  (which is implied to be the position of the secondary), the tidally deformed planet generates an additional tidal potential \,U(\mbox{{\boldmathr}},\,\mbox{{\boldmathr}}^{\,*})\,, both vectors  $r$  and \,\mbox{{\boldmathr}}^{\,*}\, being planetocentric and parametrised by their Keplerian elements $(a,e,i,{\cal M},\omega,\Omega)$ and $(a^{*},e^{*},i^{*},{\cal M}^{*},\omega^{*},\Omega^{*})$, respectively. In a situation where the secondary coincides with the perturber (and, thereby, is “feeling” the tides it itself is causing in the primary), the potential of the secondary in this field is equal to the value of \,U(\mbox{{\boldmathr}},\,\mbox{{\boldmathr}}^{\,*})\, taken for \,\mbox{{\boldmathr}}^{\,*}=\,\mbox{{\boldmathr}}\,. We, however, shall also need the gradient of the potential. To calculate it, we start out from a general expression with \,\mbox{{\boldmathr}}^{\,*}\neq\mbox{{\boldmathr}}\,, then differentiate with respect to  $r$ , and only thereafter set \,\mbox{{\boldmathr}}^{\,*}\, and  $r$  equal. Hence the tidal force (expressed in the inertial frame) acting on the perturber due to the distortion of the primary is

[TABLE]

Likewise, the tidally deformed secondary generates the additional tidal potential \,U^{\,\prime}(\,-\,\mbox{{\boldmathr}}^{\,\prime},\,\mbox{{\boldmathr}}^{\,\star})\,, where \;-\,\mbox{{\boldmathr}}^{\,\prime}\, is pointing from the centre of mass of the secondary to that of the primary, while \,\mbox{{\boldmathr}}^{\,\star}\, is pointing from the centre of mass of the secondary to that of the fictitious perturber to be identified with the primary body, both vectors being expressed in the body-fixed frame of the secondary. Accordingly, the tidal force (expressed in the inertial frame) acting on the primary due to the tidal distortion of the secondary is 777 Deriving the right-hand side expression of equation (125), we substituted \,\mbox{{\boldmathr}}^{\,\star}\, with \,-\mbox{{\boldmathr}}^{\,\star}\,. This change is acceptable because \,\mbox{{\boldmathr}}^{\,\star}\, does not show up in the final answer anyway — after the differentiation, \,\mbox{{\boldmathr}}^{\,\star}\, must be set equal to \,\mbox{{\boldmathr}}^{\,\prime}\,.

[TABLE]

We endow this force with a prime, because it emerges owing to the distortion of the secondary.

\mbox{\boldmath{\boldmath{F}}}_{0}\, being the Newtonian force (written in the inertial frame), the equations of the orbital motion in the inertial frame are: 888 Up to notation, our equations (126 - 127) agree with equations (97) in Ferraz-Mello et al. (2008). The negative signs of the arguments in the second term in our equation (125) correspond to the $\pi$-rotation in equation (99) in Ferraz-Mello et al. (2008).

[TABLE]

[TABLE]

Together, they render:

[TABLE]

or, equivalently:

[TABLE]

where we employed expressions (124) and (125). Recall that the vectors  $r$  and \,\mbox{{\boldmathr}}^{\;\prime}\, are defined in corotating frames and are related to the inertial vector \,\mbox{{\boldmath\rho}}^{\,\prime}-{\mbox{{\boldmath\rho}}}\, via formulae (5).

The contributions from the two partners enter our expression (LABEL:equation) in a symmetric manner, despite the negative signs of the arguments in the second term on the right-hand side. The easiest way to understand the origin of these negative signs is to imagine a situation where both partners are non-rotating (i.e., maintain a constant orientation with respect to an inertial frame). In this case, both rotation matrices in the definition (5) of  $r$  and \,\mbox{{\boldmathr}}^{\;\prime}\, can be chosen equal to the unity matrix, and we simply have \,\mbox{{\boldmathr}}\,=\,\mbox{{\boldmathr}}^{\;\prime}\,=\,\mbox{{\boldmath\rho}}^{\;\prime}-\,{\mbox{{\boldmath\rho}}}\,. The fact that both these vectors point from the primary to the secondary explains the difference between the arguments’ signs in the two gradients in (LABEL:equation).

Let us now write equation  (LABEL:equation) in the corotating frame of the primary, i.e., with $\,\bm{\rho}-\bm{\rho}^{\,\prime}\,$ replaced by $\,\varmathbb{R}\,\bm{r}\,$. Successive differentiations of the rotation matrix $\,\varmathbb{R}\,$ with respect to time produce the classical inertial forces:

[TABLE]

A direct comparison with equation (29) shows that the total tidal potential energy of the system is

[TABLE]

and that the disturbing function, which should be inserted in the Lagrange- or Delaunay-type planetary equations, is related to the physical potential energy via $\,{\cal{R}}\,=\,-\,{V_{1}}/{\beta}\,$. This gives us:

[TABLE]

where it is implied that in the planetary equations the differentiation of $\,{\cal{R}}\,$ should be carried out before \,\mbox{{\boldmathr}}^{\,*}\, (resp. \,\mbox{{\boldmathr}}^{\,\star}\,) is set equal to  $r$  (resp. \,\mbox{{\boldmathr}}^{\,\prime}\,).

3.2 Comparison with Kaula (1964)

In his developments, however, Kaula (1964) used the approximation

[TABLE]

Thence, Kaula’s expressions for the orbital elements’ tidal rates acquired a redundant factor of $\,M/(M+M^{\,\prime})\,$. Tolerable for the Earth-Moon system (which Kaula was having in mind), this approximation is unacceptable for a binary comprising partners of comparable masses. So, Kaula’s expressions for the rates must be multiplied by $\,(M+M^{\,\prime})/M\,$, to compensate for that oversight. 999 Aside from that, in the general case it is necessary to take into account the tidally-generated change in the orientation of the equator. As we shall see below, this will yield an additional term in the expression for $\,di/dt\,$, see equation (162).

This redundant factor of $\,M/(M+M^{\,\prime})\,$ has become a source of inaccuracy in many publications. At the same time, the overall factor is given correctly in some works, such as Ferraz-Mello et al. (2003) or Ferraz-Mello et al. (2008).

3.3 Expansion of the additional tidal potential

Let the perturber reside in the point \,\mbox{{\boldmathr}}^{\,*}\, relative to the centre of a deformable near-spherical primary. In an exterior point  $r$ , the tidally deformed body generates the additional tidal potential calculated by Kaula (1964) 101010 A partial sum of this series, with $\,l,\,|q|,\,|j|\,\leq\,2\,$ and $\,p=h=0\,$, was developed by Darwin (1879). In modern notation, his derivation is discussed by Ferraz-Mello et al. (2008). Mind that in Ibid. the convention on the meaning of the notations  $r$  and \,\mbox{{\boldmathr}}^{\,*}\, is opposite to ours.

[TABLE]

where $\,G\,=\,6.674\,\times\,10^{-11}$ m3 kg*-1* s${}^{-2}\,$ is Newton’s gravity constant. As ever, the orbital elements with and without asterisk pertain to \,\mbox{{\boldmathr}}^{*}\, and  $r$ , correspondingly, while

[TABLE]

In the expressions (134 - 135), we assume that, generally, the perturber located at \,\mbox{{\boldmathr}}^{\,*}\, does not coincide with the secondary residing in  $r$ . In the special case, when they are the same body, we must first carry out the differentiation over  $r$  and only then set \,\mbox{{\boldmathr}}^{\,*}\,=\,\mbox{{\boldmathr}}\,.

Both the dynamical Love numbers $\,k_{l}\,$ and the phase lags $\,\epsilon_{l}\,$ are functions of the tidal Fourier modes $\,\omega_{lmpq}\,=\,(l-2p)\,\dot{\omega}^{\,*}\,+\,(l-2p+q)\,n^{\,*}\,+\,m\,(\dot{\Omega}^{\,*}-\,\dot{\theta}^{\,*}\,)\,$. After the secondary and the fictitious perturber are set to be the same body, and \mbox{{\boldmathr}}^{\,*} is set equal to $r$, the modes become 111111 While in Section 2.6 we were using the osculating mean motion (defined in a standard way on the line between equations 71 and 72), here and hereafter we are using the anomalistic mean motion defined as in Table 1. We assume that the two are close, and therefore we interchangeably use the same notation n for both. The legitimacy of this is discussed in Efroimsky & Makarov (2014, Appendix B).

[TABLE]

Interestingly, Kaula himself never addressed the Fourier modes in his works, probably (mis-)assuming that both the dynamical Love numbers $\,k_{l}\,$ and the phase lags $\,\epsilon_{l}\,$ are frequency-independent. The later development of geophysics demonstrated that the forms of the frequency-dependencies of $\,k_{l}\,$ and $\,\epsilon_{lmpq}\,$ play an important role in many situations. Hence the necessity to introduce the Fourier modes $\,\omega_{lmpq}\,$ (Efroimsky 2012).

For a reader-friendly introduction to the Kaula theory, see Efroimsky & Makarov (2013). It can be understood from equation (15) in that paper, that the tidal modes’ absolute values,

[TABLE]

are the physical forcing frequencies excited in the tidally deformed body.

4 Tidal evolution of the semimajor axis

4.1 The general formula

In the Lagrange-type planetary equation for the semimajor axis rate (116), we should insert formula (132), and should perform differentiation over the mean motion. We further average the result over the mean anomaly and over the argument of pericentre as in Kaula (1964). This work, carried out in Appendix C, leads to:

[TABLE]

where we employed a shortened notation for the quality functions of the primary:

[TABLE]

the Fourier tidal modes excited in the primary being

[TABLE]

Likewise, for the quality functions of the secondary we introduced the notation

[TABLE]

the Fourier tidal modes excited in the secondary being

[TABLE]

Here $\,\Omega\,$, $\,i\,$, $\,\omega\,$ are the Euler angles of the orbit on the primary’s equator, while $\,\Omega^{\,\prime}\,$, $\,i^{\,\prime}\,$, $\,\omega^{\,\prime}\,$ are those on the secondary’s. The rotation rates of the primary and secondary are $\,\dot{\theta}\,$ and $\,\dot{\theta}^{\,\prime}\,$.

Our expression (138) differs from its counterpart in Kaula (1964) by the factor of $\,(M+M^{\,\prime})/M\,$. The reason for this is explained above, in Section 3.2.

Finally, we would mention that our expression (138) behaves well when $\,M^{\,\prime}\rightarrow 0\,$ or $\,M\rightarrow 0\,$, because $\,K_{l}\,=\,O(M^{\,2})\,$. This can be proven via formulae (31), (40b) and (42) from Efroimsky (2015).

4.2 The leading inputs

By the formulae derived in Appendix G, the quadrupole part of the major semiaxis’ rate is

[TABLE]

This long formula can obviously be split into two parts:

[TABLE]

where the first part is due to the tides in the primary and comprises the terms with $\,K_{2}(\omega_{lmpq})\,$. The second part is due to the tides in the secondary and comprises the terms with $\,K^{\,\prime}_{2}(\omega^{\,\prime}_{lmpq})\,$.

4.3 The case when the spin of neither partner is synchronised

If none of the partners is synchronised and both $\,i\,$ and $\,i^{\,\prime}\,$ are small, the leading terms are semidiurnal, i.e., those with $\,\{lmpq\}\,=\,\{2200\}\,$. Approximated with these terms, the major semiaxis’ rate is:

[TABLE]

To compare the inputs, write the above as

[TABLE]

$\rho\,$ and $\,\rho^{\,\prime}\,$ being the mean densities of the primary and the secondary, correspondingly.

When the role of the secondary is negligible, we are left with

[TABLE]

In the case when the spin is faster than orbiting, the Fourier mode $\,\omega_{2200}\,=\,2\,(n-\dot{\theta})\,$ is negative, and so is the phase lag 121212 Recall that the lag $\,\epsilon_{l}\,$ is an odd function of the Fourier mode $\,\omega_{lmpq}\,$. $\,\epsilon_{2200}\,\equiv\,\epsilon_{2}(\omega_{2200})\,$. Then the above expression becomes:

[TABLE]

where $\,Q_{2}\,$ is the quadrupole tidal quality factor defined through

[TABLE]

We reiterate that in expression (148) the quality function $\,k_{2}/Q_{2}\,$, mass $\,M\,$, and radius $\,R\,$ are those of the partner tides wherein are dominant (the primary). The mass $\,M^{\,\prime}\,$ is that of the secondary (the tides wherein we neglected in equations 147 - 148).

Approximation (146 - 147) contains only one tidal mode, the semidiurnal one. So this approximation looks the same, no matter what the frequency-dependence $\,K_{2}(\omega_{2mpq})\,$. For this reason, our equation (147) agrees with the corresponding formulae by both Hut (1981, eqn 9) and Emelyanov (2018, eqn 18) who used the CTL (constant time lag) model. It also coincides with equation (A1) from Lainey et al. (2012) who relied on the CPL (constant phase lag) tidal model.

4.4 The case when the primary is not synchronised,

while the secondary is

If the primary is not synchronised ($\,\dot{\theta}\neq n\,$), the part $\;\left({da}/{dt}\right)_{l=2}^{(prim)}\;$ is approximated with formulae (147 - 148). If the secondary is synchronised ($\,\dot{\theta}^{\,\prime}=n\,$), the terms with $\,K^{\,\prime}_{2}(2n-2\dot{\theta}^{\,\prime})\,$ and $\,K^{\,\prime}_{2}(n-\dot{\theta}^{\,\prime})\,$ in equation (4.2) get nullified. In their absence, we are left with

[TABLE]

where we took into account that $\,K^{\,\prime}_{2}\,$ is an odd function.

As expression (150) contains only one tidal mode, $\,n\,$, the form of this expression is independent of the shape of the frequency-dependence $\,K_{2}^{\,\prime}(\omega_{2mpq})\,$. So our answer coincides with equation (28) from Emelyanov (2018) and also with equation (9) from Hut (1981) if we set $\,\Omega=n\,$ in Hut’s equation. Our answer, however, differs from the first equation (A2) in the paper by Lainey et al. (2012), which contains an erroneous factor of $\;-\,21\,$ instead of $\;-\;57\,$. The same oversight is contained in equation (1) in Barnes et al. (2008) and in the first equation (40) in Shoji & Kurita (2014).

Together, the tides in both the primary and the secondary generate the rate $\,da/dt\,$ obtained by summing up the rates (147) and (150):

[TABLE]

The dissipation rate in a synchronised satellite (and the corresponding input in the elements’ rates) may be considerably amplified by longitudinal librations, when the satellite has a pronounced dynamical triaxiality (Frouard & Efroimsky 2017a, Efroimsky 2018).

4.5 Beyond quadrupole

Bills et al. (2005) argued that to attain a high precision in the modeling of Phobos’ tidal dynamics the knowledge of $\,k_{3}\,$ and perhaps even $\,k_{4}\,$ may be needed. Later, Taylor & Margot (2010) suggested that for very close asteroidal binaries the degrees $\,l\,$ up to $\,6\,$ may matter.

In the quadrupole ($\,l=2\,$) approximation (4.2), the smallest terms taken into account are of order $\,(R/a)^{5}\,e^{4}\,$. Now, if we choose to go beyond the quadrupole approximation and take into account the $\,l=3\,$ terms, the largest of those will be of order $\;(R/a)^{7}\,e^{0}\,i^{0}\,=\,(R/a)^{7}\,$. Such are the terms with $\,\{lmpq\}\,=\,\{3300\}\,$ and $\,\{3110\}\,$. We may neglect them insofar as

[TABLE]

a somewhat stringent condition not necessarily obeyed by all close-in binaries.

At the same time, had we kept in expression (4.2) only the terms up to $\,(R/a)^{5}e^{2}\,$, the neglect of the $\,l=3\,$ would be justified under a more relaxed condition

[TABLE]

4.6 Final caveat

In both Sections 4.3 and 4.4 we passingly dropped the terms containing $\,e^{4}\,K_{2}(2n)\,$ and $\,e^{4}\,K^{\,\prime}_{2}(2n)\,$, see the last line in equation (4.2). At first glance, this is legitimate when the eccentricity is not too high. In reality, the issue is subtle owing to the extremely sharp shapes of the frequency-dependencies of both $\,K_{2}(2n)\,$ and $\,K^{\,\prime}_{2}(2n)\,$. When the peak frequency happens to be equal or very close to $\,2n\,$, these terms may become prominent, even for modest values of $\,e\,$.

5 Tidal evolution of the eccentricity

The planetary equation for the eccentricity evolution is given in equation (117). The insertion of expressions (132), (185b) and (187b) in this equation and the removal of the short and long period oscillating terms gives us

[TABLE]

where we used the shortened notation (139) and (141).

As expected (see Section 3.2 above), our expression (154) differs from the corresponding formula in Kaula (1964, eqn 38) by a factor of $\,(M+M^{\,\prime})/M\,$.

The quadrupole part of expression (154) reads (see Appendix H):

[TABLE]

where the contribution $\;(de/dt)^{(prim)}\;$ comprises all the primary-generated terms (those with $\,K_{2}\,$), while $\;(de/dt)^{(sec)}\;$ comprises the secondary-generated terms (those with $\,K_{2}^{\,\prime}\,$).

5.1 The case when neither partner is synchronised

When the spin of neither body is synchronised, while both inclinations are small, the leading terms in the above equation are those linear in $\,e\;$:

[TABLE]

Specifically, when both partners satisfy the Constant Time Lag (CTL) model (i.e., when both $\,K_{2}\,$ and $\,K_{2}^{\,\prime}\,$ are linear in the tidal mode), the above expression becomes

[TABLE]

This agrees with equations (10) from Hut (1981) and (19) from Emelyanov (2018). Apart from the afore-mentioned factor of $\,(M+M^{\,\prime})/M\,$, this also agrees with the first line of equation (46) from Kaula (1964). (On the second line, Kaula lost the factor of 4 in the denominator.)

When both partners satisfy the Constant Phase Lag (CPL) model (so both $\,K_{2}\,$ and $\,K_{2}^{\,\prime}\,$ are constants) and both $\,\dot{\theta}\,$ and $\,\dot{\theta}^{\,\prime}\,$ exceed $\,3n/2\,$, we have:

[TABLE]

and

[TABLE]

wherefrom

[TABLE]

This is in agreement with the second expression in equation (A1) from Lainey et al. (2012), but differs from the corresponding formulae in some other works.

5.2 The case when the primary is not synchronised,

while the secondary is

If the primary is not synchronised ($\dot{\theta}\neq n$), the part $\,\left({de}/{dt}\right)_{l=2}^{(prim)}$ is still approximated by the no-asterisk terms from the expressions above. If at the same time the secondary is synchronised ($\dot{\theta}^{\,\prime}=n$), then in formula (156) the term with $\,K^{\,\prime}_{2}(2n-2\dot{\theta}^{\,\prime})\,$ must be set zero. Thence,

[TABLE]

where we took into consideration that $\,K^{\,\prime}_{2}\,$ is an odd function.

Containing only one frequency, expression (159) bears no dependence on the shape of the frequency-dependence of $\,K_{2}^{\,\prime}\,$. So it coincides with the corresponding expressions from Emelyanov (2018, eqn 29) and Hut (1981, eqn 10) both of whom employed the CTL model. It also is in agreement with equation (A2) by Lainey et al. (2012) who used the CPL model.

5.3 Beyond quadrupole

Under what condition may the degree-$3\,$ terms be ignored in the expression for $\,de/dt\;$?

In the quadrupole ($\,l=2\,$) approximation (155), the smallest terms taken into account are of order $\,(R/a)^{5}\,e^{3}\,$. Had it been our intention to include there also the $\,l=3\,$ terms, the largest of those would be the ones with $\,\{lmpq\}\,=\,3300\,$ and $\,3110\,$. Being of order $\;(R/a)^{7}\,e^{1}\,i^{0}\,=\,(R/a)^{7}\,e\,$, they may be ignored if

[TABLE]

Had we kept in expression (155) only the terms up to $\,(R/a)^{5}e\,$, the neglect of the $\,l=3\,$ would be justified in all realistic situations:

[TABLE]

6 Tidal evolution of the inclination

In the Lagrange-type equation for the inclination rate (118), we should insert expression (132), perform the differentiation over $\,\omega\,$, $\,\omega^{\,\prime}\,$, $\,\Omega\,$, $\,\Omega^{\,\prime}\,$ and $i^{\,\prime}$ and then extract the secular terms by averaging over ${\cal M}$, $\omega$ and $\omega^{\,\prime}$ to be consistent with Kaula (1964) who also removed the oscillating part. This work is carried out in Appendices D and E, the result being

[TABLE]

[TABLE]

Apart from the multiplier of $\,(M+M^{\,\prime})/M\,$ discussed previously (see Section 3.2), the first term of the above expression coincides with the first term of the expression (38) in Kaula (1964), while the second term in our formula differs from that in Kaula (1964, eqn 38) considerably. The second term in Kaula (1964, eqn 38) is the equivalent of the first with primed Keplerian elements. However, as explained in Section 3.2, the rotation matrix ${}^{\mathrm{T}}{\varmathbb{R}}\varmathbb{R}^{\,\prime}$ from the secondary’s frame to the primary’s frame breaks the symmetry between primed and unprimed Keplerian elements in $di/dt$. According to equation (118), the primed equivalent of the first line of equation (162) is multiplied by the slow oscillating $\cos(\omega-\omega^{\,\prime})$ and vanishes once averaged over the argument of pericentres $\omega$ and $\omega^{\,\prime}$ (see the non-averaged expression in Appendix J). Moreover, in our derivation of the inclination rate, the second line of equation (162) is an inertial force due to the non-Galilean nature of the coprecessing frame of the primary. More explicitly, this term expresses the variation of the orbital inclination with respect to the primary’s equator induced by the motion of the primary’s spin axis.

We would emphasize that the apparent lack of symmetry between the two components in the expression of $di/dt$ is due to the fact that the inclination $i$ is reckoned from the primary’s equatorial plane. Unlike $a$ and $e$ which have the same definition in both body frames and whose rates are symmetric in primed and unprimed variables, here the symmetry is recovered by writing the time derivatives of the orbital inclination with respect to the secondary’s equatorial plane, namely,

[TABLE]

[TABLE]

In addition, it shall be reminded that neither the expression (162) for $\,di/dt\,$ nor the expression (163) for $\,di^{\,\prime}/dt\,$ allows, by itself, to deduce the motion of the orbital plane with respect to the inertial frame. For example, if the inclination rate $\,di/dt\,$ happens to be zero, the orbit is fastened to the primary’s equator, and its orientation in the inertial frame is then governed by its longitude of the ascending node $\,\Omega\,$ and the primary’s Euler angles $\,(\psi,\varepsilon)\,$. Reciprocally, when $\,di/dt\,$ is non-zero, the orbit can still remain at rest in the inertial frame, in which case the apparent inclination evolution is only due the motion of the primary’s spin axis. This would have precisely been the situation if most of the total angular momentum of the system were associated with the orbital motion. Nevertheless, as in Kaula (1964), we are not interested in the motion of the orbit plane relative to the inertial frame. This is why the orientation of the orbit is measured either in the primary’s or in the secondary’s frame. The presence of two inclinations to represent the orientation of a single orbit plane might seem odd at first sight; however, the angles $\,i\,$ and $\,i^{\,\prime}\,$ can also be interpreted as the obliquities of the bodies relative to the orbit.

In a situation where the orbital angular momentum is much lower than the angular momentum of the primary’s spin, and where the dissipation in the perturber can be neglected, our result coincides with that of Kaula (1964, eqn 38), up to the afore-mentioned factor of $\,(M+M^{\,\prime})/M\,$.

The quadrupole inputs of expression (162) reads:

[TABLE]

with $\,\varrho\,=\,\beta\,n\,a^{2}/(C\,\dot{\theta})\,$. For a very small eccentricity, a cruder approximation is acceptable

[TABLE]

Without loss of precision, $\;\sin i\;$ may be changed to $\;i\;$ in both equations (164) and (165).

In many realistic situations, the inclination is stabilised in the sense that $\,di/dt\propto-i\,$. Specifically, it can be seen from equation (165) that at small eccentricities stabilisation is taking place for a synchronous orbit ($\,\dot{\theta}=n\,$), for the  2:1  spin-orbit resonance ($\,\dot{\theta}=2n\,$), and for fast prograde rotation ($\,\dot{\theta}\gg n\,$). For other values of the angular velocity, however, it is not possible to determine the sign of $\,di/dt\,$ which, generally, depends on the rheology and on the value of the parameter $\,\varrho\,$.

7 Tidal evolution of $\,\omega\,$, $\Omega$ and $\,{\cal{M}}_{0}\,$

Tidal evolution of $\,\omega\,$, $\,\Omega$ and $\,{\cal{M}}_{0}\,$ can be described by the same tools. Be mindful, though, that the rate of these angles contains an input due to the oblateness of the primary (see, e.g., Efroimsky 2005a). Therefore, even if a total rate is measured, it will not be easy to single out the tidal input in it. Also, for $\,d\omega/dt\,$ of a very close-in planet, the relativistic correction may supersede the tidal effect, like in the case of Mercury.

8 Conclusions

In this paper, we have revisited the theory of Kaula (1964) basing our calculation on a non-canonical Hamiltonian formalism with constraint. We have written down the rates $\,da/dt\,$, $\,de/dt\,$, and $\,di/dt\,$ to order $\,e^{4}\,$, inclusively. They differ from Kaula’s expressions which contain a redundant factor of $\,M/(M+M^{\,\prime})\,$, with $\,M\,$ and $\,M^{\,\prime}\,$ being the masses of the primary and the secondary. Since Kaula was interested in the Earth-Moon system, this redundant factor was close to unity and was unimportant. This omission, however, must be corrected when Kaula’s theory is applied to a binary composed of partners of comparable masses.

We have pointed out that, while it is legitimate to simply sum the primary’s and secondary’s inputs in $\,da/dt\,$ or $\,de/dt\,$, this is not the case for $\,di/dt\,$, so our expression for the inclination rate $\,di/dt\,$ differs from that of Kaula in two regards. First, in the expression for the primary’s $\,di/dt\,$ the contribution due to the dissipation in the secondary averages out completely, provided the apsidal precession is uniform. Second, we have an additional term which emerges owing to the conservation of the angular momentum: a change in the inclination of the orbit causes a change of the primary’s plane of equator.

We have carried out our developments in the gyroscopic approximation (which implies that the spin of a body is much faster than the evolution of the spin axis’ orientation). As a by-product, our work also provides a full set of equations of motion as it reads before the this approximation is applied.

Acknowledgements

One of the authors (ME) would like to thank Valeri V. Makarov and Dimitri Veras for numerous stimulating discussions.

Conflict of interest The authors hereby state that they have no conflict of interest to declare.

Appendix

Appendix A Differential of the rotation matrix $\,\varmathbb{F}={}^{\mathrm{T}}{\varmathbb{R}}\varmathbb{R}^{\prime}\,$

Let $\varmathbb{R}_{1}(\varphi)$ and $\varmathbb{R}_{3}(\varphi)$ be the rotation matrices of angle $\varphi$ around the first axis $\bm{i}={}^{\mathrm{T}}{(1,0,0)}$ and around the third one $\bm{k}={}^{\mathrm{T}}{(0,0,1)}$, respectively:

[TABLE]

The derivatives of these matrices read

[TABLE]

A direct calculation gives

[TABLE]

We introduce the hat operator which associate to any vector $\bm{v}={}^{\mathrm{T}}{(x,y,z)}\in\mathbb{R}^{3}$ the skew-symmetrix matrix $\hat{\bm{v}}$ defined as

[TABLE]

This skew-symmetric matrix is such that for any two vectors $\bm{a},\bm{b}\in\mathbb{R}^{3}$, their vector product $\bm{a}\times\bm{b}$ reads $\hat{\bm{a}}\bm{b}$. From Eq (168), one notices that

[TABLE]

Let us recall the definition of the rotation matrix $\varmathbb{F}$:

[TABLE]

Applying the rule (170), we get

[TABLE]

To simplify the result, we introduce the vectors $\bm{i}^{\prime}$, $\bm{K}$, $\bm{K}^{\prime}$ defined as

[TABLE]

The associated skew-symmetric matrices $\hat{\bm{i}}^{\prime}$, $\hat{\bm{K}}$ and $\hat{\bm{K}}^{\prime}$ are given by

[TABLE]

A direct calculation shows that $d\varmathbb{X}$, defined as $\varmathbb{R}_{1}(-i)\varmathbb{R}_{3}(-\bar{\Omega})\;(d\varmathbb{F}\,{}^{\mathrm{T}}{\varmathbb{F}})\;\varmathbb{R}_{3}(\bar{\Omega})\varmathbb{R}_{1}(i)$, is equal to

[TABLE]

Appendix B Equations of motion in terms of the Keplerian elements

Let $\,\varmathbb{M}_{\mathrm{C}/\mathrm{C^{\,\prime}}}\,$ be the Jacobian $\,\partial(\bm{p},\,\bm{r})/\partial(\bm{p}^{\,\prime},\,\bm{r}^{\,\prime})\,$ describing the transition between the two cartesian coordinate systems and given by

[TABLE]

As shown by equations (37) and (40) in the main text, the Hamiltonian equations for $\,(\bm{p},\,\bm{r})\,$ are

[TABLE]

We shall precise that in our problem ${\cal H}$ bears not dependence on $\,\bm{p}^{\,\prime}\,$, wherefore $\,\partial{\cal H}/\partial\bm{p}^{\,\prime}=\bm{0}\,$. Now, let $\,\varmathbb{M}_{\mathrm{K}/\mathrm{C}}\,$ and $\,\varmathbb{M}_{\mathrm{K}^{\,\prime}/\mathrm{C}^{\,\prime}}\,$ be the Jacobian matrices $\,\partial(\bm{Y},\bm{y})/\partial(\bm{p},\bm{r})\,$ and $\,\partial(\bm{Y}^{\,\prime},\bm{y}^{\,\prime})/\partial(\bm{p}^{\,\prime},\bm{r}^{\,\prime})\,$, respectively, where $\,\bm{Y}=(a,e,i)\,$ and $\,\bm{y}=({\cal M},\omega,\bar{\Omega})\,$. These matrices describe transitions to Keplerian (K) variables from the Cartesian (C) ones, hence the notation. Applying the chain rule, we have

[TABLE]

Then, we define the Poisson matrices $\varmathbb{B}$ and $\varmathbb{B}^{\,\prime}$ as

[TABLE]

and introduce the total Jacobian matrix

[TABLE]

Combined, equations (178-180) render us

[TABLE]

Appendix C Differentiation of $\,{\cal{R}}\,$ with respect to the mean motion

First, we differentiate $\,U\,$ with respect to the mean motion $\,{\cal{M}}\,$:

[TABLE]

After we set  $r$  and \,\mbox{{\boldmathr}}^{*}\, equal to one another (and drop the now-redundant asterisks), the  secular part  of the above derivative will become

[TABLE]

Similarly, the secular part of the derivative of the potential created by the tidally deformed secondary, at the place where the primary resides, will read:

[TABLE]

While in equation (183) $\,i\,$ is standing for the secondary’s inclination on the planetary equator, in equation (184) $\,i^{\,\prime}\,$ denotes the inclination of the primary’s apparent orbit on the secondary’s equator. Likewise, while $\,k_{l}\,\sin\epsilon_{lmpq}\,$ is the quality function of the primary, $\,k_{l}^{\,\prime}\,\sin\epsilon^{\,\prime}_{lmpq}\,$ is that of the secondary.

According to formula (132), the sum of the primary’s and secondary’s inputs in the derivative of the disturbing function over $\,{\cal{M}}\,$ will be

[TABLE]

Appendix D Differentiation of $\,{\cal{R}}\,$ with respect to the argument of the pericentre

Differentiation of $\,U\,$ over $\,\omega\,$ renders us

[TABLE]

For \,\mbox{{\boldmathr}}\,=\,\mbox{{\boldmathr}}^{*}\,, the  secular part  of the derivative reduces to

[TABLE]

Equivalently, the secular part of the derivative of the potential created by the tidally deformed secondary and acting on the primary is:

[TABLE]

where $\,i^{\,\prime}\,$, $\,\omega^{\,\prime}\,$, and $\,\Omega^{\,\prime}\,$ denote the inclination, the argument of the pericentre, and the longitude of the node of the planet’s apparent orbit as seen from the perturber.

Combining the last two equations with expression (132) for $\,{\cal{R}}\,$ as a function of $\,U\,$ and $\,U^{\,\prime}$, we obtain:

[TABLE]

Appendix E Differentiation of $\,{\cal{R}}\,$ with respect to the longitude of the node

Differentiation of $\,U\,$ over $\,\Omega\,$ gives us

[TABLE]

For \,\mbox{{\boldmathr}}\,=\,\mbox{{\boldmathr}}^{*}\,, the  secular part  of the above expression becomes

[TABLE]

Combined with equation (132), the above formula yields:

[TABLE]

In this situation,

[TABLE]

Similarly, differentiation over the longitude of the node $\,\Omega^{\,\prime}\,$ reckoned from the secondary’s equator gives

[TABLE]

Appendix F Differentiation of $\,{\cal{R}}\,$ with respect to the inclination

The derivative of $\,U\,$ with respect to $\,i\,$ is

[TABLE]

For \,\mbox{{\boldmathr}}\,=\,\mbox{{\boldmathr}}^{*}\,, the  secular part  of this expression takes the form of

[TABLE]

It should be noted that in the differentiation of $U$ with respect to $i$, the effect of the primary’s oblateness $J_{2}$ does not average out as it was the case in the differentiation over $\,\cal M\,$, $\,\omega\,$ or $\,\Omega\,$. Combining this with formula (132), we obtain:

[TABLE]

which gives

[TABLE]

Similarly, the secular part of the derivative of the potential created by the tidally deformed secondary and acting on the primary is:

[TABLE]

Appendix G Details of the calculation of $\,da/dt\,$

Writing $\,da/dt\,$ in the leading order over the inclination requires the knowledge of the squares of the inclination functions $\;F_{201}^{2}\,=\,\frac{\textstyle 1}{\textstyle 4}\,+\,O(i^{2})\;$ and $\;F_{220}^{2}\,=\,9\,+\,O(i^{2})\;$, the other $\,F^{2}_{lmp}(i)\,$ being of order $\,O(i^{2})\,$ or higher. So we shall work with the sets of integers $\;(lmpq)\,=\,(201q)\;$ and $\;(lmpq)\,=\,(220q)\;$. The corresponding eccentricity functions are:

[TABLE]

Also mind that for $\;(lmpq)\,=\,(2010)\;$ the expression $\;(2-2p+q)\;$ is zero — and so is the input $\,(da/dt)_{2010}\,$. Below is an inventory of the relevant inputs:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where we made use of notation (139) and (141).

Appendix H Details of the calculation of $\,de/dt\,$

In notation (139 - 142), our expression (154) becomes

[TABLE]

its quadrupole part being

[TABLE]

To write expression (207) in the leading order over the inclination, we shall need the squares of the two relevant $\,F_{lmp}(i)\,$ functions: $\;F_{201}^{2}\,=\,\frac{\textstyle 1}{\textstyle 4}\,+\,O(i^{2})\;$ and $\;F_{220}^{2}\,=\,9\,+\,O(i^{2})\;$, all the other $\,F^{2}_{lmp}(i)\,$ being of order $\,O(i^{2})\,$ or higher. This way, we shall be interested in the following sets of integers: $\;(lmpq)\,=\,(201q)\;$ and $\;(lmpq)\,=\,(220q)\;$. The relevant eccentricity functions are given by equations (196) above.

For $\;(lmpq)\,=\,(2010)\;$, both the expressions $\;(2-2p+q)\;$ and $\;(2-2p)\;$ vanish — and so does the $\,(de/dt)_{2010}\;$ input, up to higher-order terms in the inclinations:

[TABLE]

Thence, of the sets $\;(lmpq)\,=\,(201q)\;$, only those with $\,q\,=\,-2,\,-1,\,1,\,2\,$ are important:

[TABLE]

[TABLE]

Of the sets $\,(lmpq)\,=\,(220q)\,$, we shall be interested in the ones with $\,q\,=\,-1,\,0,\,1,\,2\;:$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where notation (139 - 142) was employed.

Appendix I Details of the calculation of $\,di/dt\,$

According to equation (118), the evolution rate of the inclination involves derivatives of the perturbing function with respect to $\,\omega\,$, $\,\omega^{\,\prime}\,$, $\,\Omega\,$, $\,\Omega^{\,\prime}\,$ and $\,i^{\,\prime}\,$. Nevertheless, in the secular expression (162), only differentiations over $\,\omega\,$ and $\,\Omega\,$ remain. These are given in Appendices D and E.

To write expression (162) in the leading order over the inclination, we should keep in mind that in this expression the squared inclination functions $\,F^{\,2}_{lmp}(i)\,$ are accompanied by a factor of either $\,\alpha_{lmp}\,=\,\left[(l-2p)\,\cos i-m\right]/\sin i\;$ or $\;\beta_{lmp}\,=\,\left[m\,\cos i-(l-2p)\right]/\sin i\,$. The functions $\,F^{\,2}_{201}(i)\,$ and $\,F^{\,2}_{220}(i)\,$ are both of order $\,O(i^{\,0})\,$, but for $\,(lmp)=(201)\,$, the two factors $\,\alpha_{201}\,$ and $\,\beta_{201}\,$ vanish. In the case $\,(lmp)=(220)\,$, we have

[TABLE]

Here we also have to consider the functions $\,F^{\,2}_{lmp}(i)\,$ of order $\,O(i^{\,2})\,$, namely $\,\displaystyle F^{\,2}_{210}\,=\,\frac{9}{4}\,\sin i^{\,2}\,+\,O(i^{\,4})\,$ and $\,\displaystyle F^{\,2}_{211}\,=\,\frac{9}{4}\,\sin i^{\,2}\,+\,O(i^{\,4})\,$, all the other $\,F^{2}_{lmp}(i)\,$ being of order $\,O(i^{\,4})\,$. The corresponding factors are $\,\alpha_{210}\,=\,1/\sin i\,+\,O(i)\,$, $\;\alpha_{211}\,=\,-\,1/\sin i\,+\,O(i)\,$, $\;\beta_{210}\,=\,-1/\sin i\,+\,O(i)\,$, $\;\beta_{211}\,=\,1/\sin i\,+\,O(i)\,$.

Therefore, we shall be interested in the following sets of integers: $\,(lmpq)\,=\,(220q)\,$, $\,(lmpq)\,=\,(210q)\,$ and $\,(lmpq)\,=\,(211q)$, the corresponding eccentricity functions being given by equations (196) above. Below we provide the resulting contributions. Deriving these, we used $\,1/\sqrt{1-e^{2}}\,=\,1\,+\,e^{2}/2\,+\,3e^{4}/8\,+\,O(e^{6})\,$ and then, in each contribution, truncated this expansion as necessary to keep the overall answer precise up to $\,e^{4}\,$, inclusively.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Without loss of precision, in all these formulae $\;\sin i\;$ may be changed to $\;i\;$.

Appendix J Long-period oscillating terms in the inclination rate

To derive the inclination rate $\,di/dt\,$ non-averaged over the long period oscillating terms, it is necessary to include the contribution due to the $\,J_{2}\,$ in the potential energy: 131313 This contribution comprises the $\,\{lmhj\}\,=\,\{20hj\}\,$ terms from the expansion for the potential energy (Frouard & Efroimsky 2017a, eqn 115).

[TABLE]

To calculate $\,di/dt\,$, we insert the above in the formula (132) for $\,\mathcal{R}\,$, then plug the result in the orbital equation (118), and finally carry out an averaging over the mean anomaly $\mathcal{M}$. This entails:

[TABLE]

with $j=2h-2$ and $r=2(h-p)+q$.