Viscoelastic tides: models for use in Celestial Mechanics

Abstract

This paper contains equations for the motion of linear viscoelastic bodies interacting under gravity. The equations are fully three dimensional and allow for the integration of the spin, the orbit, and the deformation of each body. The goal is to present good models for the tidal forces that take into account the possibly different rheology of each body. The equations are obtained within a finite dimension Lagrangian framework with dissipation function. The main contribution is a procedure to associate to each spring–dashpot model, which defines the rheology of a body, a potential and a dissipation function for the body deformation variables. The theory is applied to the Earth (solid part plus oceans) and a comparison between model and observation of the following quantities is made: norm of the Love numbers, rate of tidal energy dissipation, Chandler period, and Earth–Moon distance increase.

Appendices

Appendix 1: The kinetic energy when the angular momentum of tidal waves is relevant

The following analysis of a particular set of solutions to Eq. (12) shows that the Lagrangian function (11) sometimes misses important physical features of the dynamics.

Equation (12) have an interesting set of “tidal wave” (or “quadrupole wave”) solutions:

$$\begin{aligned} {\varvec{\varOmega }}=0,\qquad \mathbf {B}={\hat{a}}\mathbf {R}(\omega _\circ t){\mathbf {S}}\mathbf {R}^T(\omega _\circ t), \end{aligned}$$

(83)

where ${\hat{a}}>0$ is the wave amplitude,

$$\begin{aligned} \omega _\circ =\frac{1}{2}\sqrt{\frac{\gamma }{\mu }} \end{aligned}$$

(84)

is the wave angular velocity,

$$\begin{aligned} \mathbf {R}(\omega _\circ t)= \left( \begin{matrix} \cos (\omega _\circ t) &{}\quad -\sin (\omega _\circ t) &{}\quad 0\\ \sin (\omega _\circ t) &{}\quad \cos (\omega _\circ t) &{}\quad 0\\ 0 &{} 0 &{} 1 \end{matrix} \right) \qquad \text {and}\qquad {\mathbf {S}}=\left( \begin{matrix} 1 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad -1 &{}\quad 0\\ 0 &{}\quad 0 &{}\quad 0 \end{matrix} \right) . \end{aligned}$$

(85)

Since ${\varvec{\varOmega }}=0$, the body frame $\mathrm {K}$ is an inertial frame (if $\mathbf {Y}(0)=\mathbb {I}$ then $\mathrm {K}=\kappa$ for all t) and the angular momentum ${\mathrm{I}_\circ }({\varvec{\varOmega }}+\mathbf {B}{\varvec{\varOmega }}+{\varvec{\varOmega }}\mathbf {B})\in \mathrm {K}$ is null in both $\mathrm {K}$ and $\kappa$. This seems to be not physically correct since waves do usually have an associated momentum, as explained by the following interesting remarks taking from Peskin (2010). “The phenomenon of wave momentum is remarkable in several respects. First, it is not clear a priori that waves ought to have associated momentum. Waves are commonly divided into two types: transverse and longitudinal waves. In the transverse case, since nothing is moving in the direction of propagation, how can there be associated momentum in that direction? Even in the longitudinal case, since the wave motion is typically oscillatory, one would think that the average momentum density would be zero. How can there be net momentum in the direction of the wave?”. The same type of question applies to the tidal wave above, which is transverse to the radial direction. The answer is essentially the same as that given in Peskin (2010) for water waves, which are also transverse waves: “The motion of fluid particles in water waves is circular, and one might think that the net momentum in the direction of propagation would be zero. This reasoning is incorrect, however, because of the correlation between the height of the water and the direction of horizontal motion. As any swimmer knows, the water is moving forward (i.e., in the direction of the wave) at the crest of the wave, and backward in the trough. This asymmetry is the fundamental source of net momentum in the direction of wave propagation$\ldots$”

The arguments in the previous paragraph shows that a term must be added to the Lagrangian function (12) in order to a tidal wave solution to have an associated angular momentum even when $\varOmega =0$ (or even when there is no mass net rotation). The term we add is motivated by an “$\varepsilon$-expansion” of the pseudo-rigid body Lagrangian function in Ragazzo and Ruiz (2015) (it is the only term of order $\varepsilon ^{5/2}$ in that expansion). The new Lagrangian function is

$$\begin{aligned} \mathscr {L}= \frac{{\mathrm{I}_\circ }}{4}\left( |{\varvec{\varOmega }}|^2+2\ \mathrm{Tr}\ ({\varvec{\varOmega }}{\varvec{\varOmega }}^T\mathbf {B})\right) + \frac{{\mathrm{I}_\circ }}{4} \left( \mu |\dot{\mathbf {B}}|^2-\gamma |\mathbf {B}|^2\right) +{\mathrm{I}_\circ }\beta \mu \ \mathrm{Tr}\ (\mathbf {B}{\dot{\mathbf {B}}}{\varvec{\varOmega }}), \end{aligned}$$

(86)

where $\beta$ is a dimensionless constant. The Lagrangian function (86) is invariant under the action $\text {SO(3)}\times \text {SO(3)}\rightarrow \text {SO(3)}$ given by $(U,Y)\rightarrow UY$. As in the case of the rigid body the angular momentum ${\varvec{\ell }}$ defined in Eq. (10) is a conserved quantity. Again a computation gives that ${\varvec{\ell }}=\mathbf {Y}\,\mathbf {L}\,\mathbf {Y}^T$ but in this case the angular momentum $\mathbf {L}$ in the body frame $\mathrm {K}$ is given by $\mathbf {L}={\mathrm{I}_\circ }({\varvec{\varOmega }}+\mathbf {B}{\varvec{\varOmega }}+{\varvec{\varOmega }}\mathbf {B})+{\mathrm{I}_\circ }\beta \mu [{\dot{\mathbf {B}}},\mathbf {B}].$ In this expression, the first term ${\mathrm{I}_\circ }({\varvec{\varOmega }}+\mathbf {B}{\varvec{\varOmega }}+{\varvec{\varOmega }}\mathbf {B}):\mathrm {K}\rightarrow \mathrm {K}$ represents the angular momentum of the body relative to the inertial frame $\kappa$ as if it had no motion relative to $\mathrm {K}$. The second term ${\mathrm{I}_\circ }\beta \mu [{\dot{\mathbf {B}}},\mathbf {B}]$ represents the angular momentum of the body relative to the frame $\mathrm {K}$. The value of $\beta$ is directly connected to the choice of the body frame $\mathrm {K}$. In particular, unless $\beta \mu [{\dot{\mathbf {B}}},\mathbf {B}]=0$, the frame $\mathrm {K}$ is not a Tisserand frame (see Munk and MacDonald 1961).

As before, the equation for ${\varvec{\varOmega }}$ derived from the Lagrangian function (86) is

$$\begin{aligned} {\dot{\mathbf {L}}}+[{\varvec{\varOmega }},\mathbf {L}]=0 \qquad \text {where:}\qquad \mathbf {L}={\mathrm{I}_\circ }({\varvec{\varOmega }}+\mathbf {B}{\varvec{\varOmega }}+{\varvec{\varOmega }}\mathbf {B})+{\mathrm{I}_\circ }\beta \mu [{\dot{\mathbf {B}}},\mathbf {B}]. \end{aligned}$$

(87)

The equations of motion for $\mathbf {B}$ can be obtained in the following way. The set of matrices $\mathbf {B}$ can be considered as the subset of all $3\times 3$ matrices that satisfy the constraints $f_{km}(\mathbf {B})=B_{km}-B_{mk}=0$ and $f_0(\mathbf {B})=\ \mathrm{Tr}\ \mathbf {B}=0$. Let $\chi _{mk}$ and $\chi _0$ denote the Lagrange multipliers associated to these constraints. The equations of motion are obtained from the extended Lagrangian function

$$\begin{aligned} {\widehat{\mathscr {L}}}= \mathscr {L}-\chi _0 f_0-\sum _{m,k}\chi _{mk}f_{mk} \end{aligned}$$

in the usual way:

$$\begin{aligned} \frac{d}{dt}\left( \frac{\partial {\widehat{\mathscr {L}}}}{\partial \dot{\mathbf {B}}} \right) -\frac{\partial {\widehat{\mathscr {L}}}}{\partial \mathbf {B}}=0. \end{aligned}$$

Using the expressions for ${\widehat{\mathscr {L}}}$ we obtain

$$\begin{aligned} \frac{{\mathrm{I}_\circ }}{2}\left( \mu {\ddot{\mathbf {B}}}-2\beta \mu (\mathbf {B}\dot{{\varvec{\varOmega }}}+[\dot{\mathbf {B}},{\varvec{\varOmega }}]) +\gamma \mathbf {B}{\varvec{\varOmega }}^2 \right) +\chi _0\mathbb {I} +\varvec{\chi }-\varvec{\chi }^T =0. \end{aligned}$$

In order to eliminate the Lagrange multipliers $\chi _{ij}$ we add this equation to its transpose and divide by two to obtain

$$\begin{aligned} \frac{{\mathrm{I}_\circ }}{2}\left( \mu {\ddot{\mathbf {B}}}-\mu \beta ([\mathbf {B},\dot{{\varvec{\varOmega }}}]+2[\dot{\mathbf {B}},{\varvec{\varOmega }}]) +\gamma \mathbf {B}+ {\varvec{\varOmega }}^2 \right) +\chi _0\mathbb {I} =0. \end{aligned}$$

In order to determine $\chi _0$ we take the trace of this expression to get

$$\begin{aligned} \chi _0=-\frac{{\mathrm{I}_\circ }}{2}\frac{1}{3}\ \mathrm{Tr}\ {\varvec{\varOmega }}^2 . \end{aligned}$$

So, the differential equation for $\mathbf {B}$ is

$$\begin{aligned} \mu {\ddot{\mathbf {B}}}-\beta \mu [\mathbf {B},\dot{{\varvec{\varOmega }}}]-2\beta \mu [\dot{\mathbf {B}},{\varvec{\varOmega }}] +\gamma \mathbf {B}=- {\varvec{\varOmega }}^2+\frac{1}{3}(\ \mathrm{Tr}\ {\varvec{\varOmega }}^2)\mathbb {I} . \end{aligned}$$

(88)

Equations (87) and (88) have the same tidal wave solution as that in Eq. (83). The angular momentum associated to this wave is

$$\begin{aligned} {\varvec{\ell }}={\mathrm{I}_\circ }\beta \mu [{\dot{\mathbf {B}}},\mathbf {B}]={\mathrm{I}_\circ }\beta \mu 4\omega _\circ {\hat{a}}^2 \left( \begin{matrix} 0 &{} -1 &{} 0 \\ 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 \end{matrix} \right) , \end{aligned}$$

which implies that the angular momentum vector have components $\ell _1=0$, $\ell _2=0$, and $\ell _3={\mathrm{I}_\circ }\beta \mu 4\omega _\circ {\hat{a}}^2$. The energy associated to this wave is

$$\begin{aligned} E=\frac{{\mathrm{I}_\circ }}{4} \left( \mu \ \mathrm{Tr}\ {\dot{\mathbf {B}}}^2+\gamma \ \mathrm{Tr}\ \mathbf {B}^2\right) = \frac{{\mathrm{I}_\circ }}{4} \left( \mu {\hat{a}}^28\omega _\circ ^2+\gamma {\hat{a}}^2 2\right) ={\hat{a}}^2\gamma {\mathrm{I}_\circ }. \end{aligned}$$

As in Sect. 2, consider a mass m of homogeneous inviscid liquid under self-gravity. At rest the liquid has a spherical shape and moment of inertia ${\mathrm{I}_\circ }=0.4\,ma^2$, where a is the radius of equilibrium. In this case $\mu =\mu _f=1$ and $\gamma =\gamma _f=2{\mathrm{I}_\circ }G\left( \frac{5}{2}\frac{{\mathrm{I}_\circ }}{m}\right) ^{-5/2}$ as given in Eqs. (16) and (15), respectively. Suppose that the liquid is rotating uniformly with constant angular velocity $\varOmega =\varOmega _3\mathrm{e}_3$. The equation for the small free oscillations of $\mathbf {B}$ about the equilibrium shape is

$$\begin{aligned} {\ddot{\mathbf {B}}}-2\beta [\dot{\mathbf {B}},{\varvec{\varOmega }}] +\gamma _f \mathbf {B}=0. \end{aligned}$$

The decomposition of matrix $\mathbf {B}\rightarrow (b_0,b_1,b_2)$ given in Eqs. (39) and (40) can be used to rewrite this equation as:

$$\begin{aligned} {\ddot{b}}_j-2j\beta \varOmega _3i{\dot{b}}_j+\gamma _f b_j=0,\quad j=0,1,2. \end{aligned}$$

(89)

The frequency of small oscillations $\sigma$ obtained from this equation is

$$\begin{aligned} \sigma =\pm \sqrt{\gamma _f}+\beta j\varOmega _3+\mathscr {O}(\varOmega _3^2). \end{aligned}$$

Comparing this expression to that in Eq. (87) in Bryan (1889) we obtain that for a mass of perfect liquid

$$\begin{aligned} \beta =\beta _f=\frac{1}{2}. \end{aligned}$$

(90)

For several types of waves (electromagnetic waves, sound waves, water waves, and certain kinds of traveling waves on strings under tension) the following simple relation holds (Peskin 2010)

$$\begin{aligned} \text {wave momentum}\cdot \text {wave phase velocity}=\text {wave energy}. \end{aligned}$$

For a wave moving along a circular path, like the tidal wave above, this relation could be changed as

$$\begin{aligned} \underbrace{\text {momentum}\cdot R}_{angular\ momentum} \cdot \overbrace{\frac{1}{R}\cdot \text {phase velocity}}^{angular\ velocity}=\text {energy} \end{aligned}$$

where $R>0$ is some reference radius. If we assumed that the previous relation would hold for the tidal wave solution then, using the expressions for $\ell _3$ and E above, we would obtain $\beta =1$. Therefore, it seems that $\beta \ge 0$ can vary depending on the rheological nature of the body.

To finish, consider the tide response equations [compare to Eq. (37)]:

$$\begin{aligned} \begin{aligned} \mu {\ddot{\mathbf {B}}}-2\beta \mu [\dot{\mathbf {B}},{\varvec{\varOmega }}] +\gamma \mathbf {B}+{\varvec{\varLambda }}&={\mathbf {A}}(t)\\ {\ddot{{\varvec{\varLambda }}}}+c_1{\dot{{\varvec{\varLambda }}}}+c_2{\varvec{\varLambda }}&=c_3{{\ddot{\mathbf {B}}}}+c_4{\dot{\mathbf {B}}}, \end{aligned} \end{aligned}$$

(91)

where ${\dot{{\varvec{\varOmega }}}}=0$ and $\varOmega =\varOmega _3\mathrm{e}_3$. If we apply the decomposition $\mathbf {B}\rightarrow (b_0,b_1,b_2)$, ${\mathbf {A}}\rightarrow (a_1,a_2,a_3)$, and ${\varvec{\varLambda }}\rightarrow (\lambda _0,\lambda _1,\lambda _2)$ given in Eqs. (39) and (40) we rewrite Eq. (91) as [compare to Eq. (49)]:

$$\begin{aligned} \begin{aligned} \mu {\ddot{b}}_j-2j\beta \mu \varOmega _3i{\dot{b}}_j+\gamma b_j+\lambda _j&=a_j(t)\\ {\ddot{\lambda }}_j+c_1{\dot{\lambda }}_j+c_2\lambda&=c_3{\ddot{b}}_j+c_4{\dot{b}}_j\qquad j=0,1,2. \end{aligned} \end{aligned}$$

(92)

As before $a_j(t)$ can be Fourier expanded. For a term of the form $a_j(t)={\hat{a}}\mathrm{e}^{i(\omega t+\varphi )}$, there corresponds a solution of the form $b_j(t)={\hat{b}}\mathrm{e}^{i(\omega t+\varphi )}$, $\lambda _j(t)={\hat{\lambda }}\mathrm{e}^{i(\omega t+\varphi )}$, such that ${\hat{a}}$, ${\hat{b}}$, and ${\hat{\lambda }}$ are related as

$$\begin{aligned} \begin{aligned} \{\mu \omega (-\omega +2j\beta \varOmega _3)+\gamma \} {\hat{b}}+{\hat{\lambda }}&={\hat{a}}\\ (-\omega ^2 +c_1\omega i+c_2)\hat{\lambda }&=(-c_3\omega ^2 +c_4\omega i) {\hat{b}}\qquad j=0,1,2. \end{aligned} \end{aligned}$$

For the diurnal and the semi-diurnal frequencies $\omega =j\varOmega _3$, $j=1,2$, the first of these equations becomes

$$\begin{aligned} \{\mu (2\beta -1)(j\varOmega _3)^2+\gamma \} {\hat{b}}+{\hat{\lambda }} ={\hat{a}}, \end{aligned}$$

that is an equation equal to that obtained in Sect. 5.1 for the same frequencies except for the change $\mu \rightarrow \mu (1-2\beta )$. This is very convenient in the fit of the parameters $c_1,c_2,c_3,c_4$ as we did in Sect. 6.2. Essentially, the same procedure used in that section can be used in this case and since $\mu (1-2\beta )$ can be either positive or negative this relaxes the positivity condition we had on the parameter $\mu$.

Appendix 2: The planetary dissipation of tidal energy and the average work done by the tidal force: relations to the dissipation function and to the Lagrangian function, respectively

In this appendix we present a classical formula of Zschau and Platzman for the time-average planetary dissipation rate of tide energy W. We apply this formula to our model and show that it gives the time average of the dissipation function of the model, which shows that our choice of dissipation function is consistent with results obtained from continuum mechanics. The clear and elegant deduction of the formula of Zschau and Platzman in Platzman (1984) lead us to an analogous, apparently new, formula for the time-average of the work done by the primary tidal force in deforming the planet. This new formula, deduced in Appendix “The average work done by the tidal force and the Lagrangian function” section, gives the “time-average planetary tidal balance of potential and kinetic energy”⁹ J. This new formula is related to the Lagrangian function of our model in the same way as the Zschau and Platzman formula is related to the dissipation function.

1.1 Appendix 2.1: A formula of Zschau and Platzman for the planetary dissipation function and the Rayleigh dissipation function

Following Platzman (1984), let $\overline{\varPsi }$ be the primary astronomical potential (due to the satellite), $\varPsi ^\prime$ be the secondary potential due to the tidal redistribution of the mass of the planet, and $\varPsi =\overline{\varPsi }+\varPsi ^\prime$ be the complete tide potential. Let

$$\begin{aligned} \varPsi =\sum _n\left[ \overline{\varPsi }_n\left( \frac{r}{a}\right) ^n+ \varPsi _n^\prime \left( \frac{a}{r}\right) ^{n+1}\right] \end{aligned}$$

(93)

be the solid harmonic decomposition of $\varPsi$, where: $\overline{\varPsi }_n$ and $\varPsi _n^\prime$ are spherical surface harmonics of degree n, $(r,\theta ,\phi )$ are spherical coordinates in the rotating reference frame $\mathrm {K}$, and a is the radius of a spherical surface $S_a$ that contains the planet inside and does not contain the satellite. In Platzman (1984, Eq. (6)), it is shown that the time-average of the planetary dissipation of tidal energy is given by

$$\begin{aligned} W=\frac{1}{4\pi G a}\sum _n(2n+1)\int _{S_{a}}\langle \overline{\varPsi }_n \partial \varPsi _n^\prime /\partial t\rangle dS, \end{aligned}$$

(94)

where dS is the surface element and $\langle \cdot \rangle$ denotes the average over a period tide. Since $\overline{\varPsi }$ is dominated by $\overline{\varPsi }_2$, all terms in this formula but those of second-degree can be neglected

$$\begin{aligned} W=\frac{5}{4\pi G a}\int _{S_{a}}\langle \overline{\varPsi }_2 \partial \varPsi _2^\prime /\partial t\rangle dS. \end{aligned}$$

(95)

The energy dissipation rate W is related to the Rayleigh dissipation function $\mathscr {D}$ in Eq. (23) as follows.

In our model, the primary potential has only the the second-degree harmonic

$$\begin{aligned} \overline{\varPsi }(x)=-\frac{1}{2}(x\cdot {\mathbf {A}}x)= \frac{|x|^2}{a^2}\overline{\varPsi }_2({\hat{x}}), \end{aligned}$$

(96)

where $x=|x|{\hat{x}}$. The secondary potential is

$$\begin{aligned} \varPsi ^\prime (x)=-\frac{3}{2}\frac{G{\mathrm{I}_\circ }}{|x|^5}(x\cdot \mathbf {B}x) =\frac{a^3}{|x|^3}\varPsi _2^\prime ({\hat{x}}). \end{aligned}$$

(97)

Using that $dS=a^2\sin \theta d\theta d\phi =a^2d{\hat{x}}$ on $S_a$, where $(\theta , \phi )$ are polar coordinates on the unit sphere, Eq. (95) can be written as

$$\begin{aligned} W= & {} \frac{5}{4\pi G a}\int _{S_{a}}\langle \overline{\varPsi }_2 \partial \varPsi _2^\prime /\partial t\rangle dS\\= & {} \frac{15}{16\pi }{\mathrm{I}_\circ }\left\langle \int _{S_1}({\hat{x}}\cdot {\mathbf {A}}{\hat{x}}) ({\hat{x}}\cdot {\dot{\mathbf {B}}}{\hat{x}}) d{\hat{x}} \right\rangle \\= & {} \frac{15}{16\pi }{\mathrm{I}_\circ }\left\langle \sum _{ijkl}A_{ij}{\dot{B}}_{kl}\underbrace{ \int _{S_1}{\hat{x}}_i{\hat{x}}_j{\hat{x}}_k{\hat{x}}_l d{\hat{x}}}_{=\varGamma _{ijkl}} \right\rangle . \end{aligned}$$

Now, the same argument as that given at the end of Sect. 3 implies that $\varGamma _{ijkl}$ is a linear combination of the three tensors $\delta _{ij}\delta _{kl}$, $\delta _{ik}\delta _{jl}$, and $\delta _{il}\delta _{jk}$. An easy computation using polar coordinates shows that

$$\begin{aligned} \varGamma _{ijkl}=\frac{4\pi }{15}(\delta _{ij}\delta _{kl}+ \delta _{ik}\delta _{jl}+\delta _{il}\delta _{jk}) \end{aligned}$$

and

$$\begin{aligned} W=\frac{15}{16\pi }{\mathrm{I}_\circ }\left\langle \sum _{ijkl}A_{ij}{\dot{B}}_{kl}\varGamma _{ijkl} \right\rangle =\frac{{\mathrm{I}_\circ }}{2}\langle \ \mathrm{Tr}\ {\mathbf {A}}{\dot{\mathbf {B}}}\rangle = \frac{{\mathrm{I}_\circ }}{2}\frac{1}{T}\int _0^T\ \mathrm{Tr}\ {\mathbf {A}}{\dot{\mathbf {B}}} dt. \end{aligned}$$

(98)

In principle the time average must be taken over a period tide. It happens that the tide is not exactly periodic but only almost periodic usually with a dominance of the semi-diurnal mode, so the average is still well defined as $T\rightarrow \infty$. For the computations, T must be taken as a large multiple of the semi-diurnal period.

Finally, multiplying both sides of the first equation in (36) by ${\dot{\mathbf {B}}}$, taking the trace and the averaging of both sides, and using that the average of a time derivative is zero $\langle \frac{df}{dt}\rangle =0$, we obtain

$$\begin{aligned} \langle \ \mathrm{Tr}\ {\varvec{\varLambda }}{\dot{\mathbf {B}}}\rangle =\langle \ \mathrm{Tr}\ {\varvec{\varLambda }}_1{\dot{\mathbf {B}}}\rangle + \langle \ \mathrm{Tr}\ {\varvec{\varLambda }}_2 {\dot{\mathbf {B}}}\rangle = \langle \ \mathrm{Tr}\ {\mathbf {A}}{\dot{\mathbf {B}}}\rangle . \end{aligned}$$

Multiplying both sides of the second and third equations in (36) by ${\varvec{\varLambda }}_1$ and ${\varvec{\varLambda }}_2$, respectively, taking the trace and the averaging of both sides, we obtain

$$\begin{aligned} \langle \ \mathrm{Tr}\ {\mathbf {A}}{\dot{\mathbf {B}}}\rangle = \langle \ \mathrm{Tr}\ {\dot{\mathbf {B}}}{\varvec{\varLambda }}_1\rangle + \langle \ \mathrm{Tr}\ {\dot{\mathbf {B}}}{\varvec{\varLambda }}_2\rangle = \eta _1^{-1}\langle |{\varvec{\varLambda }}_1|^2\rangle + \eta _2^{-1}\langle |{\varvec{\varLambda }}_2|^2\rangle = \eta _1\langle |\dot{{\tilde{\mathbf {B}}}}_1|^2\rangle + \eta _2\langle |{\dot{{\tilde{\mathbf {B}}}}}_2|^2\rangle , \end{aligned}$$

where we used that ${\varvec{\varLambda }}_1=\eta _1\dot{{\tilde{\mathbf {B}}}}_1$ and ${\varvec{\varLambda }}_2=\eta _2\dot{{\tilde{\mathbf {B}}}}_2$, similarly to $\eta _1{\dot{{\tilde{x}}}}_1=\lambda _1$ and $\eta _2{\dot{{\tilde{x}}}}_2=\lambda _2$ in Eq. (20). This last equation and Eqs. (98) and (23) imply that

$$\begin{aligned} W= \frac{{\mathrm{I}_\circ }}{2}\langle \ \mathrm{Tr}\ {\mathbf {A}}{\dot{\mathbf {B}}}\rangle = \frac{1}{T}\int _0^T(2\mathscr {D}) dt. \end{aligned}$$

(99)

This equation shows that the planetary dissipation of tidal energy as given in Eq. (95) coincides with the average power dissipated by the planet in our model.

1.2 Appendix 2.2: The average work done by the tidal force and the Lagrangian function

There are two fundamental functions in our theory of tide response: the Lagrangian $\mathscr {L}_{\mathbf {B}}$ and the dissipation $\mathscr {D}$ functions both given in Eq. (23). After we showed the relation between the formula of Zschau and Platzman in Eq. (95) and the dissipation function $\mathscr {D}$, we wandered whether there would exist a similar relation for $\mathscr {L}_{\mathbf {B}}$. Indeed there is, $\mathscr {L}_{\mathbf {B}}$ is equal to the time-average of the work done by the tidal force to deform the planet or equivalently to the “time-average planetary tidal balance of energy”(see footnote at the first page).

Let ${\tilde{J}}$ be the work done by the primary tidal force to deform the planet

$$\begin{aligned} {\tilde{J}}=-\int _{B_a}\rho (d\cdot \nabla \overline{\varPsi })dV \ , \end{aligned}$$

where: $B_a$ is a solid ball that contains the planet and does not contain the satellite, d is the displacement vector, $\rho$ is the density of the deformed state, and $-\rho \nabla \overline{\varPsi }$ is the primary tidal force [see Lambeck (1980), p. 7 for details]. It is convenient to suppose that the density is a smooth function that vanishes outside the planet, so the nontrivial part of the integral is restricted to the planet interior. The result we obtain below is the same as that obtained supposing that $\rho$ varies discontinuously at some interfaces [see Platzman (1984) for the treatment of discontinuities]. Partial integration, the use of $\nabla ^2\varPsi ^\prime =-4\pi G\nabla (\rho d)$ [see Lambeck (1980), p. 7], and that $\rho =0$ over the spherical surface $S_a$ leads to

$$\begin{aligned} {\tilde{J}}=-\frac{1}{4\pi G}\int _{B_a}\overline{\varPsi }\nabla ^2\varPsi ^\prime dV. \end{aligned}$$

Partial integration twice and the use of $\nabla ^2\overline{\varPsi }=0$ in $B_a$ implies that

$$\begin{aligned} {\tilde{J}}=-\frac{1}{4\pi G}\left\{ \int _{S_a}\overline{\varPsi }(\nabla \varPsi ^\prime \cdot \mathrm{e}_r)dS -\int _{S_a}\varPsi ^\prime (\nabla \overline{\varPsi }\cdot \mathrm{e}_r)dS\right\} {.} \end{aligned}$$

Then, by using the solid harmonic decomposition of $\varPsi ^\prime$ and $\overline{\varPsi }$ in Eq. (93) the integral becomes

$$\begin{aligned} {\tilde{J}} =\frac{1}{4\pi G a}\sum _n(2n+1)\int _{S_{a}}\overline{\varPsi }_n \varPsi _n^\prime dS. \end{aligned}$$

Finally, taking the time average of this equation we obtain

$$\begin{aligned} J=\langle {\tilde{J}}\rangle =\frac{1}{4\pi G a}\sum _n(2n+1)\int _{S_{a}}\langle \overline{\varPsi }_n \varPsi _n^\prime \rangle dS, \end{aligned}$$

(100)

that is similar to Eq. (94) of Zschau and Platzman. Since $\overline{\varPsi }$ is dominated by $\overline{\varPsi }_2$, all terms but those of second-degree can be neglected in the above formula, so

$$\begin{aligned} J=\frac{5}{4\pi G a}\int _{S_{a}}\langle \overline{\varPsi }_2 \varPsi _2^\prime \rangle dS, \end{aligned}$$

(101)

that is similar to Eq. (95).

The same reasoning that lead to Eq. (98) shows that the average planetary tidal balance of energy can be written as

$$\begin{aligned} J =\frac{{\mathrm{I}_\circ }}{2}\langle \ \mathrm{Tr}\ {\mathbf {A}}\mathbf {B}\rangle = \frac{{\mathrm{I}_\circ }}{2}\frac{1}{T}\int _0^T\ \mathrm{Tr}\ {\mathbf {A}} \mathbf {B}dt. \end{aligned}$$

(102)

Finally, multiplying both sides of the equation of motion (36) by $\frac{{\mathrm{I}_\circ }}{2} \mathbf {B}$, using the constraints $\mathbf {B}=\mathbf {B}_1+{\tilde{\mathbf {B}}}_1$ and $\mathbf {B}=\mathbf {B}_2+{\tilde{\mathbf {B}}}_2$, using the relations ${\varvec{\varLambda }}_1=\alpha _1\mathbf {B}=\eta _1\dot{{\tilde{\mathbf {B}}}}_1$ and ${\varvec{\varLambda }}_2=\alpha _2\mathbf {B}=\eta _2\dot{{\tilde{\mathbf {B}}}}_2$, which are similar to those in Eq. (20), and using Eq. (102) we obtain

$$\begin{aligned} J=\frac{1}{T}\int _0^T \frac{{\mathrm{I}_\circ }}{2}\Bigl \{-\mu |{\dot{\mathbf {B}}}|^2+\gamma |\mathbf {B}|^2 +\alpha _1|\mathbf {B}_1|^2+\alpha _2 |\mathbf {B}_2|^2\Bigr \} dt= \frac{1}{T}\int _0^T(-2\mathscr {L}_{\mathbf {B}})dt, \end{aligned}$$

(103)

that is similar to Eq. (99). Therefore, J is related to the Lagrangian function $\mathscr {L}_{\mathbf {B}}$ given in Eq. (23) in the same way as the rate of dissipation of tidal energy W is related to the dissipation function $\mathscr {D}$ given in Eq. (23).