Celest Mech Dyn Astr
DOI 10.1007/s10569-016-9682-3
ORIGINAL ARTICLE
Dynamic Elastic Tides
Jack Wisdom1 · Jennifer Meyer1,2
Received: 17 November 2015 / Revised: 26 February 2016 / Accepted: 7 March 2016
© Springer Science+Business Media Dordrecht 2016
Abstract This is an exploration of dynamic tides on elastic bodies. The body is thought of
as a dynamical system described by its modes of oscillation. The dynamics of these modes
are governed by differential equations that depend on the rheology. The modes are damped
by dissipation. Tidal friction occurs as exterior bodies excite the modes and the modes act
back on the tide raising body. The whole process is governed by a closed set of differential
equations. Standard results from tidal theory are recovered in a two-timescale approximation
to the solution of these differential equations.
Keywords Tides · Tidal evolution · Tidal friction · Lags · Orbital evolution · Wobble
damping · Dynamical astronomy
1 Introduction
Solar System studies of tides have, for the most part, been based directly or indirectly on the
works of Darwin (1908). In these works the tide raising potential is Fourier expanded, and
the response of the planet to each Fourier component is separately specified. The exterior
potential of the tidally distorted planet takes the form of the tide raising potential, but it
is phase delayed and given an amplitude proportional to the potential Love number, which
encapsulates the structure of the planet. The large number of parameters is made tractable by
choosing the phase delays to be proportional by a factor that has some particular frequency
dependence. This frequency dependence has been either chosen for convenience or motivated
by the behavior of materials at nontidal frequencies. For instance, if the tidal phase lags
are taken to be proportional to frequency then the various tidal components form a single
tidal bulge that is delayed in time from the tide raising potential (Mignard 1981; Hut 1981;
B Jack Wisdom
wisdom@mit.edu
1 Massachusetts Institute of Technology, Cambridge, MA 02139, USA
2 Present Address: Clovis Community College, Fresno, CA 93730, USA
123
J. Wisdom, J. Meyer
Eggleton et al. 1998). If the phase lags are taken to be independent of the frequency then we
get the tidal models of Kaula (1966) and Goldreich (1963, 1966). There is scant knowledge
of the actual frequency dependence of tides on planetary bodies. Williams (2008) determined
the phase lags for the Moon at two frequencies. He reports a Q of approximately 30 for a
forcing frequency of one month and a Q of about 35 for one year. Note that Williams actually
measures the phase lag δ, and converts the phase lag to Q using the relation 1/Q = tan δ
(Williams 2011, personal communication).
The method of Darwin (and those who follow him) is peculiar in that the process of raising
tides is not described by a simple set of differential equations. If we want to calculate the
tidal effect on an exterior body, we first let the body raise a set of tidal bulges (for each
Fourier component), we impose a phase delay in the response by a specified amount for each
component, then we calculate the effect of the phase delayed bulges on the exterior body. The
end result may be a set of differential equations that model tidal friction, but the derivation
of those equations is not through the solution of a closed set of differential equations.
We take a different approach, inspired by, but distinct from, the theory of dynamic tides
on stars (Zahn 1975). We think of the planet as a dynamical system described by its modes
of oscillation. The dynamics of the modes are determined by the internal distribution of
density and the rheology. Dissipation occurs entirely in the modal dynamics, which may
be more complicated than a simple second order differential equation modelling a damped
harmonic oscillator, but is nevertheless an ordinary differential equation. Time varying tidal
potentials excite the modes and the modes respond according to their own dynamics. In turn,
the exterior potential of the modally distorted planet affects the tide raising body. The whole
process is described by a closed set of differential equations. Multiple time scale analysis of
these differential equations allow us to recover Darwin-like approximate expressions for tidal
evolution and compare to standard results. Our framework is logically clean. The dynamics
of the modes determines all the phase lags in the approximate expressions, but we need not
go to these approximate expressions—we can directly integrate the full equations. Our goal
is to develop a framework into which it will be possible to incorporate different rheologies
in a logically clear and modular fashion. This paper describes the framework and works out
the details for one particular model of dissipation in the modes. For simplicity we assume
the body is homogeneous with uniform material properties.
The plan of the paper is as follows. In the next section we develop the expressions for the
distortion of a body subject to an external tidal potential. We then develop the expressions for
the tidal potential due to an exterior perturbing body. We then derive the free elastic modes
of our body. This is followed by a decomposition of the tidal distortion in terms of the free
elastic modes. We treat the modes as a dynamical system and next give the Lagrangian for the
motion of the free modes. The Lagrangian for the motion of the planet with tidal distortion is
developed in the following two sections. First, we develop expressions for the kinetic energy
of a planet with modal distortion, then we develop the potential energy. The full Lagrangian is
then formed and the equations ofmotion are developed. As a first application, we calculate the
effect of elasticity on the frequency of the Chandler wobble. This is followed by a section on
dissipation: we incorporate dissipation by modifying the differential equations governing the
dynamics of the modes. This is the point at which one would incorporate different rheologies.
We then calculate the rate of damping of the Chandler wobble. We come next to tidal friction.
We do this in stages. We first calculate the rate of change of semimajor axis and rotation rate
for a body in a circular orbit about a extended body with zero obliquity. We then extend these
results to an eccentric orbit. We consider both the case where the rotation of the body is in
spin-orbit resonance and where it is not. Then finally we extend the results to the general
case with eccentricity and obliquity. This is followed by our concluding remarks.
123
Dynamic Elastic Tides
2 Static tidal distortion
In this section we develop expressions for the distortion of a body subject to an external tidal
potential. Our principal source for this section is Love (1944). We assume the body is elastic
and incompressible.
The external disturbing potential is
UT (x, y, z) =
∑
lm
wml X¯
m
l (x/R, y/R, z/R), (1)
where wml are the coefficients of solid spherical harmonics X¯
m
l [see Eq. (303)], and R is the
radius of the body. The external disturbing potential leads to a distortion of the body: this
distortion contributes to the disturbing potential. The total potential can be written (Jefferys
1976) V = ∑lm V ml where
V ml (x, y, z) = wml X¯ml (x/R, y/R, z/R) +
3g
2l + 1
m
l X¯
m
l (x/R, y/R, z/R). (2)
The second term arises from the potential of the distortion of the body. Given wml we must
solve for ml .
The displacement of a material particle in the body is u(t, x, y, z), where the coordi-
nates (x, y, z) refer to the position in the body of material particle before displacement. The
equations to be solved are
− ∇ p + μ∇2u + ρ∇V = 0. (3)
Incompressibility implies
∇ · u = 0. (4)
Taking the divergence of the first equation and using the second equation, we see that ∇2 p =
0. We then put
p(x, y, z) = ρV (x, y, z) + p˜(x, y, z), (5)
where
p˜(x, y, z) =
∑
lm
pml X¯
m
l (x/R, y/R, z/R), (6)
and where X¯ml is a solid spherical harmonic. Let u1 be
u1(x, y, z) =
∑
lm
(
Aml r
2∇ X¯ml (x/R, y/R, y/R) + Bml xX¯ml (x/R, y/R, y/R)
)
, (7)
where X¯ml is, again, a solid harmonic. Using Eqs. (313–315), we find that u1 of the form of
Eq. (7) satisfies Eqs. (3) and (4) if the following two equations are satisfied:
pml /μ = (4l + 2)Aml + 2Bml (8)
0 = 2l Aml + (l + 3)Bml , (9)
respectively. For l = 2, these have the solution
Am2 =
5pm2
42μ
(10)
Bm2 = −
4pm2
42μ
. (11)
123
J. Wisdom, J. Meyer
To this solution we may add arbitrary solutions of the equation
∇2u = 0, (12)
with
∇ · u = 0. (13)
It will be enough to add solutions of the form
u2(x, y, z) =
∑
lm
f ml ∇ X¯ml (x/R, y/R, z/R). (14)
The boundary conditions at the surface (see Love) require
0 = (l Aml + Bml )R2 pml + l f ml − Rml (15)
0 = (2l Aml + Bml )R2 pml + 2(l − 1) f ml (16)
ρwml = (μ(2l Aml + (l + 2)Bml ) − 1)pml + gρ
(
1 − 3
2l + 1
)
ml . (17)
The solutions for l = 2 are
pm2 = −(21/2)(μ/gR)wm2 /Δ (18)
f m2 = 2(R/g)wm2 /Δ (19)
m2 = (5/2)(wm2 /g)/Δ, (20)
where
Δ = 1 + 19μ/(2ρgR). (21)
The displacement Love number is h2 = (5/2)/Δ. The total displacement u1 + u2, for this
m, is
u(x, y, z) = u mT (x, y, z)h2wm2 /(gR), (22)
defining the tidal shape function u mT . For m = 0, the tidal shape function u 0T (0, 0, R) has
the components (0, 0, R), and the displacement at the surface along the z axis, the axis of
maximum tidal distortion, is h2w02/g.
3 Tidal disturbing potential
In this section we develop expressions for the tidal potential due to an external perturber. The
tidal disturbing potential is the second harmonic contribution to the gravitational potential
energy per unit mass
U2(x, R′) = − Gm
(R′)3
(
3
2
(x · Rˆ′)2 − 1
2
(x · x)
)
, (23)
where x has components (x, y, z), m is the mass of the perturbing, tide-raising body, and
the R′ is the vector from center of mass of the body to the tide raising body. The distance
between the bodies is R′, and the direction to the disturbing body is Rˆ′. Let α, β, and γ
be the direction cosines of Rˆ′ = αxˆ + β yˆ + γ zˆ. Then, let N be the rotation that takes
the direction Rˆ′ to the zˆ direction: zˆ = N Rˆ′. Let N (θ, ψ) = Mx (θ)Mz(ψ), where Mi are
123
Dynamic Elastic Tides
rotations about the indicated axes. Then, θ = acos(γ ) and ψ = atan(α, β).1 Using the
property (Nx) · (Ny) = x · y, we see
U2(x, R′) = U2(Nx, z) (24)
=
∑
m
wm2 (α, β, γ )X¯
m
2 (x/R, y/R, z/R). (25)
We find wm2 (α, β, γ ) = W X¯m2 (α, β, γ ), where W = −(Gm/R′)(R/R′)2. The correspond-
ing tidal shape function can be obtained with an appropriate rotation of u 0T . We have
(N (α, β, γ ))−1u 0T (N (α, β, γ )x) =
∑
m
X¯m2 (α, β, γ )u
m
T (x). (26)
4 Elastic free modes
Next we develop expressions for the free modes of oscillation. Our principal source for this
section is Lamb (1882). In his terminology modes corresponding to tidal distortions are
vibrations in the second class. We will assume the material is incompressible. Our notation
differs slightly from his.
The displacement is described as a sum over modes. The displacement is
u(t, x, y, z) =
∑
lmn
RAmlnφ
m
ln(x, y, z) cos(ω
n
l t + mln), (27)
where Amln and 
m
ln are the amplitude and phase of each mode. The modes of interest in the
tidal problem are l = 2, m runs from −l to l, and n = 1, 2, . . .. The amplitude Amln and the
function φmln are dimensionless.
Let
φmln(x, y, z)/R (28)
= (Υln + ψml−1(κlnr))∇ X¯ml −
l
l + 1
κ2lnr
2l+3
(2l + 1)(2l + 3)ψl+1(κlnr)∇
(
X¯ml
r2l+1
)
, (29)
where
ψn(θ) = (−1)n(2n + 1)!!
(
1
θ
d
dθ
)n ( sin θ
θ
)
, (30)
where κln and Υln are explained below. The functions ψn are related to the spherical Bessel
functions.
The modal frequencies are determined by the condition
aldl − blcl = 0, (31)
where, in the incompressible case,
al = (κln R)
2
2l + 1 − 2(l − 1) (32)
bl = 1 (33)
1 We use the two argument arctangent with the property ψ = atan(sin(ψ), cos(ψ)), rather than the one
argument arctangent with the property ψ = atan(tan ψ), in order to retain quadrant information.
123
J. Wisdom, J. Meyer
cl = 2(l − 1)ψl−1(κln R) − (κln R)
2
2l + 1 ψl(κln R) (34)
dl = l
l + 1
(
ψl(κln R) + 2(l + 2)
κln R
ψ ′l (κln R)
)
(35)
where ψ ′(x) = dψ(x)/dx . Note that κln occurs only in the combination κln R, determined
by Eq. (31). The frequencies ωln of the modes satisfy
κ2ln = ω2lnρ/μ, (36)
where ρ is the density and μ is the rigidity. Numerically, we find the lowest frequency
tidal mode has κ21R/π = 0.8484938956, the next lowest frequency mode has κ22 R/π =
1.7421226796, and then κ23R/π = 2.8257142846. Finally,
Υln = (κln R)
2ψl(κln R) − 2(l − 1)(2l + 1)ψl−1(κln R)
(κln R)2 − 2(l − 1)(2l + 1) . (37)
5 Representation of tidal distortion by elastic modes
Define the overlap integral
〈u1, u2〉 = 1
V
∫
V
(u1 · u2)dV (38)
where V = (4/3)π R3, the volume of the body. The elastic modes have zero overlap. Let
(β2n)
2 = 〈φm2n, φm2n
〉
. (39)
Note that β2n is independent of m. We find: β21 = 0.5325432017, β22 = 0.2145631038,
and β23 = 0.0826504258. Define the l = 2 normalized elastic modes:
uMnm = R φm2n
/
β2n . (40)
The tidal shape function u 0T can be expanded in terms of the normalized elastic modes with
m = 0. Let
u 0T (x, y, z) =
∑
n
gnuMn0(x, y, z), (41)
then
gn =
〈
u 0T , un
〉
. (42)
We find g1 = 0.5608256130, g2 = −0.0381757369, and g3 = 0.0039974227. The tidal
shape function and its representation in terms of the elastic modes is shown in Fig. 1.
If the tidal distortion is rotated, then the representing modes rotate accordingly. For a
disturbing body with direction cosines (α, β, γ ) the tidal distortion is
u(x, y, z) = h2W
gR
∑
m
X¯m2 (α, β, γ )u
m
T (x, y, z) (43)
= h2W
gR
∑
mn
gn X¯
m
2 (α, β, γ )u
M
nm(x, y, z). (44)
123
Dynamic Elastic Tides
Fig. 1 The solid line shows the z
component of the tidal shape
function u 0T . The dashed line
shows the representation of the
tidal shape function using the
n = 1 elastic mode. The dotted
line shows the representation
using the n = 1 and n = 2 modes.
The three mode representation is
indistinguishable from the solid
line on this scale
z/R
(u
0 T
(0
,0
,z
))
z
/
R
1.00.50.0-0.5-1.0
1.0
0.5
0.0
-0.5
-1.0
6 Lagrangian for elastic modes
We treat each mode as a degree of freedom. We assume here that the configuration of the
body is given as a sum of modal distortions (with l = 2). We write the displacement as
u(t, x, y, z) =
∑
mn
Πnm(t)uMnm(x, y, z) (45)
where Πnm(t) is the dimensionless modal coordinate at time t .
The kinetic energy in each mode is
Tnm(t,Πnm, Π˙nm) = 1
2
MR2Π˙2nm . (46)
with M = ρV , the mass of the body. Recall that the modes are normalized. We know the
frequency of each mode, so we can write the elastic potential energy for each mode:
Vnm(t,Πnm, Π˙nm) = 1
2
MR2ω22nΠ
2
nm (47)
The Lagrangian for the free elastic modes is then
L(t,Π, Π˙) = 1
2
MR2
∑
mn
(
Π˙2nm − ω22nΠ2nm
)
. (48)
This gives the equation of motion
MR2(Π¨nm + ω22nΠnm) = 0, (49)
confirming a free oscillation with frequency ω2n . The modes will be forced when the full
potential energy is developed.
123
J. Wisdom, J. Meyer
7 Kinetic energy
The next task is to develop the kinetic energy of the rotating body, taking into account the
fact that the shape is changing. We assume the undistorted starting configuration of the body
is a sphere of radius R and densityρ, with moment of inertia (2/5)MR2. The configuration
of the body at time t is obtained by distorting and rotating this reference body. We give the
body a distortion uA that gives the principal moments of inertia and a time-dependent modal
tidal distortion uT (t), followed by a time-dependent rotation M(t) in space. The position of
each constituent is
xα(t) = M(t)(x0 + uA + uT (t)). (50)
The velocity of the constituent is
x˙α(t) = M˙(t)(x0 + uA + uT (t)) + M(t)(u˙T (t)) (51)
= M˙(t)M(t)−1M(t)(x0 + uA + uT (t)) + M(t)(u˙T (t)) (52)
= ω(t) × (M(t)(x0 + uA + uT (t))) + M(t)(u˙T (t)) (53)
= M(t)(ω′(t) × (x0 + uA + uT (t)) + u˙T (t)) (54)
The square of the velocity is
(x˙α(t))2 = (ω′(t) × (x0 + uA)) · (ω′(t) × (x0 + uA)) (55)
+ 2(ω′(t) × (x0 + uA)) · (ω′(t) × uT (t) + u˙T (t)) (56)
+ (ω′(t) × uT (t) + u˙T (t)) · (ω′(t) × uT (t) + u˙T (t)) (57)
The distortion uA that gives the principal moments is volume preserving; it can be repre-
sented as a gradient
uA(x, y, z) = aR2∇ X¯02(x/R, y/R, z/R) + bR2∇ X¯22(x/R, y/R, z/R), (58)
where
a = (B + A)/2 − C
A + B + C (59)
b = (
√
3/2)(B − A)
A + B + C . (60)
The sum of the principal moments is 3I where I = (2/5)MR2, the moment of inertia of a
homogeneous sphere. In these expressions we are ignoring second order contributions to the
moments in a and b; we are assuming a and b are small. Note that we can write
R∇ X¯02 = (−X¯11,−X¯−11 , 2X¯01). (61)
and
R∇ X¯22 = (
√
3 X¯11,−
√
3 X¯−11 , 0). (62)
For convenience, introduce
u0A(x, y, z) = (x, y, z) (63)
u1A(x, y, z) = R2∇ X¯02(x/R, y/R, z/R) (64)
u2A(x, y, z) = R2∇ X¯22(x/R, y/R, z/R). (65)
123
Dynamic Elastic Tides
The first term in the kinetic energy integral is, by design,
1
2
∫
V
ρ((ω′(t) × (x0 + uA)) · (ω′(t) × (x0 + uA)))dV (66)
= 1
2
(
A(ωa)2 + B(ωb)2 + C(ωc)2
)
, (67)
where the components of ω′ are (ωa, ωb, ωc).
Introduce the integrals
I (1)knmi j =
1
V
∫
V
(ei × ukA) · (e j × uMnm)dV, (68)
I (2)knmi =
1
V
∫
V
(ei × ukA) · uMnmdV, (69)
I (3)nmn′m′i j =
1
V
∫
V
(ei × uMnm) · (e j × uMn′m′)dV, (70)
I (4)nmn′m′i =
1
V
∫
V
(ei × uMnm) · uMn′m′dV . (71)
In terms of these integrals, the complete kinetic energy is
T (t; q,Π;ω′, Π˙) = 1
2
(
A(ωa)2 + B(ωb)2 + C(ωc)2
)
(72)
+ MR2
∑
knmi j
ak I
(1)
knmi j (ω
′)i (ω′) jΠnm (73)
+ MR2
∑
knmi
ak I
(2)
knmi (ω
′)i Π˙nm (74)
+ 1
2
MR2
∑
nmn′m′i j
I (3)nmn′m′i j (ω
′)i (ω′) jΠnmΠn′m′ (75)
+ MR2
∑
nmn′m′i
I (4)nmn′m′i (ω
′)iΠnmΠ˙n′m′ (76)
+ 1
2
MR2
∑
nm
(
Π˙nm
)2
, (77)
where a0 = 1, a1 = a, and a2 = b, and q are the coordinates that specify the orientation of
the body in space. For k = 0, the nonzero coefficients are:
I (1)0,n,0,0,0 = I (1)0,n,0,1,1 = Cn (78)
I (1)0,n,0,2,2 = −2Cn (79)
I (1)0,n,1,0,2 = I (1)0,n,2,0,0 = I (1)0,n,−1,1,2 = I (1)0,n,−2,0,1 = −
√
3Cn (80)
I (1)0,n,2,1,1 =
√
3Cn, (81)
with C1 = 0.1747793752, C2 = −0.0501544874, and C3 = 0.0138166088, and with
I (1)0nmi j = I (1)0nmji . These have a simple representation I (1)0nmi j = −Cn∂i∂ j X¯m2 (ωa, ωb, ωc).
123
J. Wisdom, J. Meyer
8 Potential energy
The gravitational potential energy is the integral of Eq. (23) over the mass of the body:
VG(t; q,Π;ω′, Π˙) =
∫
V
ρU2(x, R′)dV (82)
It is convenient to introduce a number of integrals, and then write the potential energy in
terms of them. Let
J (1)knm =
1
V
∫
V
ukA · u MnmdV, (83)
J (2)nmn′m′ =
1
V
∫
V
u Mnm · u Mn′m′dV = δnn′δmm′ (84)
J (3)knmi j =
1
V
∫
V
(eˆi · ukA)(eˆ j · u Mnm)dV, (85)
J (4)nmn′m′i j =
1
V
∫
V
(eˆi · u Mnm)(eˆ j · u Mn′m′)dV . (86)
In terms of these the gravitational potential energy is
VG(t; q,Π;ω′, Π˙) (87)
= − Gm
(R′)3
(
(1 − 3(α′)2)A + (1 − 3(β ′)2)B + (1 − 3(γ ′)2)C
2
)
(88)
+ GmM
R′
(
R
R′
)2 ∑
knm
J (1)knmakΠnm (89)
+ 1
2
GmM
R′
(
R
R′
)2 ∑
nmn′m′
J (2)nmn′m′ΠnmΠn′m′ (90)
− 3GmM
R′
(
R
R′
)2 ∑
i jknm
J (3)knmi jα
′
iα
′
j akΠnm (91)
− 3
2
GmM
R′
(
R
R′
)2 ∑
i jnmn′m′
J (4)nmn′m′i jα
′
iα
′
jΠnmΠn′m′ , (92)
where α′0 = α′, α′1 = β ′, α′2 = γ ′, and (α′, β ′, γ ′) are the direction cosines of the perturbing
body with respect to the body axes:
(α′, β ′, γ ′) = (M(q))−1 (α, β, γ ), (93)
and (α, β, γ ) are the direction cosines of the perturbing body with respect to a spatially fixed
rectangular basis, and M(q) is the rotation that carries the body in its reference orientation
(with principal axes aligned with the spatial axes) to the actual orientation specified by the
coordinates q .
The coefficients J (1)0nm are all zero, but there are some nonzero k = 0 terms among the
J (3)knmi j . Specifically,
J (3)0,n,0,0,0 = J (3)0,n,0,1,1 = −Cn (94)
J (3)0,n,0,2,2 = 2Cn (95)
123
Dynamic Elastic Tides
J (3)0,n,2,0,0 = J (3)0,n,−2,0,1 = J (3)0,n,1,0,2 = J (3)0,n,−1,1,2 =
√
3Cn (96)
J (3)0,n,2,1,1 = −
√
3Cn (97)
with Cn as before, and J
(3)
i j0nm = J (3)j i0nm . These have the simple representation J (3)0nmi j =
Cn∂i∂ j X¯m2 (α
′, β ′, γ ′).
Recall that the elastic potential energy, following Eq. (47), is
VE (t; q,Π;ω′, Π˙) = 1
2
MR2
∑
nm
ω22nΠ
2
nm . (98)
The total potential energy is V = VG + VE .
9 Lagrangian and equations of motion
The Lagrangian for the system is the difference of the kinetic energy T and the potential
energy V . Note that the kinetic energy does not depend on q and the potential energy does
not depend on ω′. Further, the potential energy depends on q only through (α′, β ′, γ ′).
It will be convenient to use Euler-like angles to specify the orientation of the body with
respect to its reference orientationwith the principal axes (aˆ, bˆ, cˆ) alignedwith the rectangular
spatial axes (xˆ, yˆ, zˆ). (These are the principal axes of the body without tidal distortion.) We
choose the rotation that carries the body to its actual orientation as
M(θ, φ, ψ) = Mz(φ)Mx (θ)My(ψ). (99)
Note that these are not the usual Euler angles.
The equations of motion for the angular velocities (ωa, ωa, ωc) and the coordinates
(θ, φ, ψ) are Poincaré’s equations (Poincaré 1901). The equations of motion for the veloci-
ties Π˙ and the coordinates Π are the Lagrange equations. In the derivation of the Poincaré
equations, we use the vector field basis corresponding to infinitesimal rotations about the
body axes (aˆ, bˆ, cˆ). In the chosen coordinates, the basis vector fields are
ea = cosψ ∂
∂θ
− sin ψ
cos θ
∂
∂φ
+ sin ψ sin θ
cos θ
∂
∂ψ
(100)
eb = ∂
∂ψ
(101)
ec = sin ψ ∂
∂θ
+ cosψ
cos θ
∂
∂φ
− cosψ sin θ
cos θ
∂
∂ψ
(102)
The commutators of these basis fields satisfy
[ei , e j ]( f ) =
∑
i jk
cki j ek( f ), (103)
with structure constants
cki j = i jk, (104)
where i jk is 1 for (i, j, k) equal to a cyclic permutation of (0, 1, 2), −1 for a cyclic permu-
tation of (2, 1, 0), and 0 otherwise. Concretely,
[ea, eb] = ec (105)
[eb, ec] = ea (106)
123
J. Wisdom, J. Meyer
[ec, ea] = eb. (107)
The Lagrangian depends on the coordinates (θ, φ, ψ) only through (α′, β ′, γ ′). In the
Poincaré equations we need ei (Lˆ), where Lˆ is the Lagrangian written in terms of the angular
velocities, and the vector field takes the derivative of the coordinate slot of the Lagrangian.
We have
ea(Lˆ) = γ ′ ∂ Lˆ
∂β ′
− β ′ ∂ Lˆ
∂γ ′
(108)
eb(Lˆ) = α′ ∂ Lˆ
∂γ ′
− γ ′ ∂ Lˆ
∂α′
(109)
ec(Lˆ) = β ′ ∂ Lˆ
∂α′
− α′ ∂ Lˆ
∂β ′
. (110)
The Poincaré equations are, in this case,
d
dt
(
∂ Lˆ
∂ωa
)
= γ ′ ∂ Lˆ
∂β ′
− β ′ ∂ Lˆ
∂γ ′
+ ωc ∂ Lˆ
∂ωb
− ωb ∂ Lˆ
∂ωc
(111)
d
dt
(
∂ Lˆ
∂ωb
)
= α′ ∂ Lˆ
∂γ ′
− γ ′ ∂ Lˆ
∂α′
+ ωa ∂ Lˆ
∂ωc
− ωc ∂ Lˆ
∂ωa
(112)
d
dt
(
∂ Lˆ
∂ωc
)
= β ′ ∂ Lˆ
∂α′
− α′ ∂ Lˆ
∂β ′
+ ωb ∂ Lˆ
∂ωa
− ωa ∂ Lˆ
∂ωb
(113)
For our Lagrangian, including the I (1)knmi j terms in the kinetic energy with k = 0, but ignoring
the other I (i) contributions, we find
∂ Lˆ
∂ωa
= Aωa − 2MR2
∑
nm
Cn∂0 X¯
m
2 (ω
a, ωb, ωc)Πnm (114)
∂ Lˆ
∂ωb
= Bωb − 2MR2
∑
nm
Cn∂1 X¯
m
2 (ω
a, ωb, ωc)Πnm (115)
∂ Lˆ
∂ωc
= Cωc − 2MR2
∑
nm
Cn∂2 X¯
m
2 (ω
a, ωb, ωc)Πnm (116)
Similarly, keeping the k = 0 terms in J (3)knmi j , but ignoring the other J (i) terms, we find
∂ Lˆ
∂α′
= − Gm
(R′)3
3Aα′ + MR2
∑
nm
FnΠnm∂0 X¯m2 (α′, β ′, γ ′) (117)
∂ Lˆ
∂β ′
= − Gm
(R′)3
3Bβ ′ + MR2
∑
nm
FnΠnm∂1 X¯m2 (α′, β ′, γ ′) (118)
∂ Lˆ
∂γ ′
= − Gm
(R′)3
3Cγ ′ + MR2
∑
nm
FnΠnm∂2 X¯m2 (α′, β ′, γ ′), (119)
where
Fn = Gm
(R′)3 (
6Cn) ≈ Gm
(R′)3
(
Rω22nh2gn
g
)
, (120)
123
Dynamic Elastic Tides
where the last form of Fn is valid for large rigidity. The Poincaré equations are constructed
from these components.
Keeping the same terms as before, the Lagrange equations governing the motion of Πnm
are
Π¨nm + ω22nΠnm = −2Cn X¯m2 (ωa, ωb, ωc) + Fn X¯m2 (α′, β ′, γ ′)). (121)
Note that if we ignore the terms involving Π Poincaré’s equations become
A
dωa
dt
= (B − C)ωbωc − 3Gm
(R′)3
(B − C)β ′γ ′ (122)
B
dωb
dt
= (C − A)ωaωc − 3Gm
(R′)3
(C − A)α′γ ′ (123)
C
dωc
dt
= (A − B)ωaωb − 3Gm
(R′)3
(A − B)α′β ′, (124)
which are just Euler’s equations for the motion of a rigid body subject to a gravity-gradient
torque.
10 Chandler wobble
Elasticity affects the period of the Eulerian wobble of the Earth (Chandler 1891; Newcomb
1882). As an illustration of the formalism we have developed we will calculate the elastic
correction to the frequency of the Eulerian wobble. We assume the body is rotating very
nearly with its spin axis aligned with the cˆ principal axis, and that the body is axisymmetric
A = B. We assume there is no external perturbing body. For simplicity we take into account
only the gravest elastic modes (those with n = 1). We assume here that the elastic modes are
in equilibrium, so, from Eq. (121), we have
Π1m = −2C1 X¯m2 (ωa, ωb, ωc)/ω221. (125)
We will assume ωa and ωb are much smaller than ωc; we will systematically ignore terms
that are second and higher order in ωa and ωb. With this assumption
Π1,0 = −2C1 (ω
c)2
ω221
(126)
Π1,1 = −2
√
3C1
ωaωc
ω221
(127)
Π1,−1 = −2
√
3C1
ωbωc
ω221
(128)
Π1,2 = 0 (129)
Π1,−2 = 0. (130)
We have, in this case,
∂ Lˆ
∂ωa
= Aωa + 2C1MR2Π1,0ωa − 2
√
3C1MR
2Π1,1ω
c (131)
∂ Lˆ
∂ωb
= Bωb + 2C1MR2Π1,0ωb − 2
√
3C1MR
2Π1,−1ωc (132)
123
J. Wisdom, J. Meyer
∂ Lˆ
∂ωc
= Cωc − 4C1MR2Π1,0ωc. (133)
The Poincaré equations, with these assumptions, show that ωc is constant, so Π1,0 is also
constant. We define
A′ = A + 2C1MR2Π10 = A − 4C21MR2
(
ωc
ω21
)2
(134)
B ′ = B + 2C1MR2Π10 = B − 4C21MR2
(
ωc
ω21
)2
(135)
C ′ = C − 4C1MR2Π10 = C + 8C21MR2
(
ωc
ω21
)2
, (136)
and assume that (A′, B ′, C ′) are the observed principal moments of inertia. The Poincaré
equations are then
A′ω˙a
(
1 + 12C21
MR2
A′
(
ωc
ω21
)2)
= ωbωc(B ′ − C ′) + 12C21
(
ωc
ω21
)2
MR2ωbωc
(137)
B ′ω˙b
(
1 + 12C21
MR2
B
(
ωc
ω21
)2)
= ωaωc(C ′ − A′) − 12C21
(
ωc
ω21
)2
MR2ωaωc
(138)
Taking the time derivative of the first equation and using the second, we find
ω¨a = −ω2Cωa, (139)
where the Chandler frequency is
ωC = ωc
(
C ′ − A′
A′
− 12C21
MR2
A
(
ωc
ω21
)2)
, (140)
and we have ignored a quantity of second order. We see that the Eulerian frequency, ωc(C ′ −
A′)/A′, is reduced by elasticity. As so far developed, our theory does not apply to a body such
as theEarthwhere the body has significant radial variation in density.Nevertheless,we can see
what effective rigidity a homogeneous Earth must have to get the observed Chandler period of
434 days. Using a density of 5500kg/m3 and a radius of 6,371,000m, C ′/(MR2) = 0.3307,
and A′/(MR2) = 0.3296, we find μ = 1.8 × 1011 N/m2, which is actually comparable to
modern estimates of the rigidity of the Earth at a depth of about 1000km (Stacey 1992).
11 Dissipation
The Q of an oscillator is defined as
1
Q
= 1
2π E
∮
d E
dt
dt, (141)
where d E/dt is the rate at which work is done on the oscillator, E is the energy stored in the
oscillator, and the integral is over one cycle of the oscillation. Consider the damped driven
oscillator
m(x¨ + dx˙ + ω20x) = A cos(ωt). (142)
123
Dynamic Elastic Tides
The forced response is
x(t) = B cos(ωt − δ), (143)
where
B2 =
(
A
m
)2 1
(ω20 − ω2)2 + (ωd)2
. (144)
Let
δ0 = arctan(ωd, ω20 − ω2), (145)
then for ω < ω0 the phase shift δ = δ0 is positive, meaning that the response of the oscillator
lags the forcing, and B is positive. For ω > ω0, the phase shift δ0 is greater than π/2. We
could take δ = δ0 for all ω and always keep B positive. Alternatively, we can bring δ into
the range of −π/2 to π/2 (the range of the one argument arctan) by subtracting π from δ0
making δ = δ0 − π negative. Then we flip the sign of B to be negative, so that the solution
remains valid. So for ω < ω0, both δ and B are positive, and for ω > ω0 both δ and B are
negative. This reflects more clearly that the response is out of phase with the forcing.
The energy dissipated per cycle is
ΔE = π AB sin δ, (146)
which is positive. For now, we take the peak energy stored in the oscillator to be
E = mω20 B2/2, (147)
and return to this definition later. The Q is then given by
1
Q
= ((ω
2
0 − ω2)2 + (ωd)2)1/2
ω20
sin δ (148)
≈ ((ω
2
0 − ω2)2)1/2
ω20
sin δ (149)
≈ sin δ, (150)
where the first approximation is for small dissipation (small d), and the second approximation
additionally assumes small forcing frequency (ω  ω0). In the same limit
tan δ ≈ ωd/ω20. (151)
So for large Q
1
Q
≈ ωd
ω20
. (152)
Note that in this model, with these definitions, and in these limits the Q is inversely propor-
tional to ω.
The definition of the energy stored is somewhat problematic. In Eq. (147) we took the
energy stored to be the maximum potential energy of the spring, the potential energy of the
spring at its maximum extension. Goldreich (1963) took a different definition: he defined
E to be the integral over the first quarter cycle of the rate at which work is done on the
oscillator:
E =
∫ T/4
0
Fvdt (153)
= −
∫ T/4
0
(A cos(ωt)ωB sin(ωt − δ))dt (154)
123
J. Wisdom, J. Meyer
= ABω
(
cos(2ωt − δ)
4ω
− t sin δ
2
)∣∣∣∣
T/4
0
(155)
= AB
2
(
cos(δ) − π
2
sin δ
)
(156)
Ignoring the second term term in E for small δ, we find that
1
Q
= tan δ (157)
More exactly, there is a second term, which Goldreich ignores as an approximation. Peale
(2011, personal communication) argues that the second term should be excluded in general
because the energy stored in the oscillator should be the integral of the conservative rate of
work and the second term corresponds to dissipative work. Thus Eq. (157) should be regarded
as exact. We may compare this to the result we derived above, Eq. (149), in which 1/Q is
approximately sin δ and then only for small d and ω  ω0. For this model, to distinguish
these is splitting hairs because the phase lags are not known well enough to distinguish
between tan δ and sin δ. But in other models of dissipation, different definitions of the energy
stored lead to more qualitatively different results.
We incorporate tidal dissipation into our modal model by adding an ad hoc dissipative
term to the modal equations, Eq. (121),
Π¨nm + dnΠ˙nm + ω22nΠnm = −2Cn X¯m2 (ωa, ωb, ωc) + Fn X¯m2 (α′, β ′, γ ′). (158)
We assume dn is independent of m. We have
1
Qn
≈ ωdn
ω22n
, (159)
where ω is the forcing frequency, assumed to be much less than the modal frequency.
The fact that this model gives phase lags proportional to frequency makes it equivalent to
popularly used sets of tidal evolution equations (Mignard 1981; Hut 1981). The disadvantage
is that real materials and planets do not behave this way. The frequency dependence of
the phase lags for the Moon are different (Williams 2008). In a subsequent paper we will
incorporate different rheologies into our model.
12 Wobble damping
Here we consider the decay of the Eulerian wobble due to dissipation in the elastic modes.
We again assume that the spin axis is nearly aligned with the cˆ axis. The damped modal
equations (keeping only n = 1 modes) are
Π¨1,0 + d1Π˙1,0 + ω221Π1,0 = −2C1(ωc)2 (160)
Π¨1,1 + d1Π˙1,1 + ω221Π1,1 = −2
√
3C1ω
cωa (161)
Π¨1,−1 + d1Π˙1,−1 + ω221Π1,−1 = −2
√
3C1ω
cωb. (162)
As ωc is constant, we can assume
Π1,0 = −2C1(ωc/ω21)2, (163)
123
Dynamic Elastic Tides
is constant. To construct the Poincaré equations we need
∂ Lˆ
∂ωa
= Aωa + 2C1MR2Π1,0ωa − 2
√
3C1MR
2Π1,1ω
c (164)
∂ Lˆ
∂ωb
= Bωb + 2C1MR2Π1,0ωb − 2
√
3C1MR
2Π1,−1ωc (165)
ignoring second order quantities. Define (A′, B ′, C ′) as before, and assume these are the
observed values of the principal moments. Then the left hand sides of the Poincaré equations
are
d
dt
∂ Lˆ
∂ωa
= A′ω˙a − 2√3C1MR2ωcΠ˙1,1 (166)
d
dt
∂ Lˆ
∂ωa
= B ′ω˙b − 2√3C1MR2ωcΠ˙1,−1 (167)
and the right hand sides are
ωc
∂ Lˆ
∂ωb
− ωb ∂ Lˆ
∂ωc
= ωcωb(B ′ − C ′) − 2√3C1MR2(ωc)2Π1,−1 (168)
ωa
∂ Lˆ
∂ωc
− ωc ∂ Lˆ
∂ωa
= ωcωa(C ′ − A′) + 2√3C1MR2(ωc)2Π1,1. (169)
Collecting the equations, the wobble damping is governed by the constant coefficient linear
differential equations:
0 = Π¨1,1 + d1Π˙1,1 + ω221Π1,1 + ξωcωa (170)
0 = Π¨1,−1 + d1Π˙1,−1 + ω221Π1,−1 + ξωcωb (171)
0 = A′ω˙a − ξωcMR2Π˙1,1 − ωcωb(B ′ − C ′) + ξωcMR2ωcΠ1,−1 (172)
0 = B ′ω˙b − ξωcMR2Π˙1,−1 − ωcωa(C ′ − A′) − ξωcMR2ωcΠ1,1 (173)
where ξ = 2√3C1. We find the solution decays exponentially, proportional to e−t/τ , with
1
τ
= ξ2
(
ωc
ω21
)4 (MR2
A′
) (ωC
ωc
)
d1. (174)
Using Eq. (152), we define the effective Q for an oscillation at the Chandler frequency to
satisfy
1
QC
= ωCd1
ω221
, (175)
then
1
τ
= ξ2 (ω
c)3
(ω21)2
(
MR2
A′
) (
1
QC
)
. (176)
The wobble damping timescale is approximately
τ ≈ 19.38 A
MR2
μQC
ρR2(ωc)3
. (177)
This agrees with Peale (1973), though the method of calculation is very different.
123
J. Wisdom, J. Meyer
13 Tidal friction
Consider a perturbing body in orbit about a dissipative elastic body. We assume the spin
axis of the body is perpendicular to the orbit plane (γ ′ = 0) and that there is no wobble
(ωa = ωb = 0). We will use the same coordinates (θ, φ, ψ) to specify the orientation as
before. Let n be the mean motion of the perturbing body, with true longitude λ. We assume
that A′ = B ′. We will consider only the gravest elastic modes (those with n = 1). With these
assumptions the equations of motion are
C ′ω˙c = √3MR2F1
(
Π1,22α
′β ′ − Π1,−2
(
(α′)2 − (β ′)2)) (178)
ω˙a = ω˙b = 0 (179)
with
Π¨1,2 + d1Π˙1,2 + ω221Π1,2 = F1
√
3
2
(
(α′)2 − (β ′)2) (180)
Π¨1,−2 + d1Π˙1,−2 + ω221Π1,−2 = F1
√
3
2
2α′β ′ (181)
Π1,1 = Π1,−1 = 0, (182)
and with
Π1,0 ≈ −(2C1(ωc)2 + F1/2)/ω221. (183)
The average polar moment is
C ′ = C + 8C21MR2
(
ωc
ω21
)2
+ 2C1MR2 F1
ω221
. (184)
For this geometry α′ = cos(λ − φ) and β ′ = sin(λ − φ), with ωc = φ˙. Define φ′ = φ − λ,
then the equations of motion are
C ′ω˙c = √3MR2F1
(−Π1,2 sin(2φ′) − Π1,−2 cos(2φ′)
)
(185)
Π¨1,2 + d1Π˙1,2 + ω221Π1,2 = F1
√
3
2
cos(2φ′) (186)
Π¨1,−2 + d1Π˙1,−2 + ω221Π1,−2 = −F1
√
3
2
sin(2φ′) (187)
Let’s assume that the rotation is not near synchronous, and that φ′ moves approximately
uniformly
φ′ ≈ (ωc − n)t. (188)
(We are taking the orbital eccentricity to be zero.) Then the Π equations are periodically
forced damped harmonic oscillators, with solutions (assuming ω21 	 |ωc − n|)
Π1,2 =
√
3
2
F1
ω221
cos(2(ωc − n)t − δ2ωc−2n) (189)
Π1,−2 = −
√
3
2
F1
ω221
sin(2(ωc − n)t − δ2ωc−2n) (190)
where
tan δ2ωc−2n ≈ 2(ω
c − n)d1
ω221
= 1
QT
, (191)
123
Dynamic Elastic Tides
defining QT as the effective Q at the frequency 2(ωc−n). Substituting these into the equation
for ω˙c we find
ω˙c = −3
2
MR2
C ′
(F1)2
ω221
sin δ2ωc−2n . (192)
This gives the rate of deceleration of rotation for a perturber of mass m in a circular orbit. It
can apply to either the case of the deceleration of a planet by a satellite or of the rotation of
a satellite by the primary. This expression is easily generalized using Eqs. (144–148).
To reduce this expression to the usual expression we have to make an additional assump-
tion. For one of the factors of F1, we use the approximate form for large rigidity in Eq. (120).
We also use k2 = (3/5)h2. We find then
ω˙c = −3
2
k2
QT
(
R
R′
)3 (Gm
a3
) (
m
M
)
MR2
C ′
χ, (193)
where
χ = 10C1g1 ≈ 0.98. (194)
In each case m is the tidal perturber, and M is the perturbed solid body. Except for the χ
factor Eq. (193) gives the expression in Goldreich and Peale (1966). Note that the expression
given by Peale (1973) is too big by a factor of 2, as mentioned by Dobrovolskis (2007). Also
note if we include more modes then the factor of χ is closer to 1:
∑∞
i=1 10Ci gi = 1.
This formula can be applied to get the tidal decay of rotation for a non-synchronous
satellite, but also for the tidal evolution of semimajor axis. Assuming the angular momentum
of the system is conserved, the angular momentum lost by the rotating body goes into the
orbit. We find, using Eq. (193), that the rate of tidal evolution of the semimajor axis a of a
circular orbit is
1
a
da
dt
= 3n k2
QT
R5
a5
m
M
, (195)
which agrees with Peale (1986).
14 Tidal friction with eccentricity
Consider a perturbing body in an elliptic orbit about a dissipative elastic body. To start, we
will assume the spin axis of the body is perpendicular to the orbit plane (γ ′ = 0) and that
there is no wobble (ωa = ωb = 0). Let n be the mean motion of the perturbing body, a
the semimajor axis, e the orbital eccentricity,  the longitude of pericenter, and λ the true
longitude. Following, Eqs. (185–187), the equations of motion are
C ′ω˙c = √3MR2F1
(−Π1,2 sin(2(φ − λ)) − Π1,−2 cos(2(φ − λ))
)
(196)
Π¨1,2 + d1Π˙1,2 + ω221Π1,2 = F1
√
3
2
cos(2(φ − λ)) (197)
Π¨1,−2 + d1Π˙1,−2 + ω221Π1,−2 = −F1
√
3
2
sin(2(φ − λ)), (198)
Π¨1,0 + d1Π˙1,0 + ω221Π1,0 = −(2C1(ωc)2 + F1/2), (199)
where both F1 and λ vary nonuniformly because of the Keplerian orbital motion. Recall that
F1 ∝ (a/R′(t))3 (see Eq. 120).
123
J. Wisdom, J. Meyer
To solve these equations we Fourier expand the periodic forcing. We use
(
a
R′(t)
)3
cos(k f (t)) =
∞∑
j=−∞
X−3,kj (e) cos( jnt), (200)
and (
a
R′(t)
)3
sin(k f (t)) =
∞∑
j=−∞
X−3,kj (e) sin( jnt), (201)
where f is the true anomaly and X−3,kj (e) are Hansen functions.
Substituting these into the Π equations of motion gives
Π¨1,2 + d1Π˙1,2 + ω221Π1,2 = G1
√
3
2
∞∑
j=−∞
X−3,2j (e) cos(2φ − 2 − jnt) (202)
Π¨1,−2 + d1Π˙1,−2 + ω221Π1,−2 = −G1
√
3
2
∞∑
j=−∞
X−3,2j (e) sin(2φ − 2 − jnt) (203)
Π¨1,0 + d1Π˙1,0 + ω221Π1,0 = −2C1(ωc)2 − G1
1
2
∞∑
j=−∞
X−3,0j (e) cos( jnt), (204)
where φ = φ0 + ωct (with φ0 giving the orientation of the planet at time t = 0) and
G1
(
a
R′(t)
)3
= F1. (205)
These linear equations have the approximate solutions
Π1,2 =
√
3
2
G1
ω221
∞∑
j=−∞
X−3,2j (e) cos(2φ − 2 − jnt − δ2ωc− jn) (206)
Π1,−2 = −
√
3
2
G1
ω221
∞∑
j=−∞
X−3,2j (e) sin(2φ − 2 − jnt − δ2ωc− jn) (207)
Π1,0 = −2C1 (w
c)2
ω221
− 1
2
G1
ω221
∞∑
j=−∞
X−3,0j (e) cos( jnt − δ jn) (208)
where
tan δ2ωc− jn ≈ (2ω
c − jn)d1
ω221
, (209)
tan δ jn ≈ jnd1
ω221
. (210)
In writing these expressions we have used the approximation that |2ωc − jn| << ω21. These
expressions are easily generalized.
Substituting these expressions into the equation for ω˙c and averaging over the orbital
period, we find
ω˙c = −3
2
MR2
C ′
(G1)2
ω221
∞∑
j=−∞
(
X−3,2j (e)
)2
sin δ2ωc− jn, (211)
123
Dynamic Elastic Tides
where we have ignored the small time variation in C ′ which exists because F1 (see Eq. 184)
is now time dependent. Note that X−3,20 (e) = 0. This expression looks more familiar if we
assume large rigidity:
ω˙c = −3
2
k2
(
R
R′
)3 (Gm
a3
) ( m
M
) MR2
C ′
∞∑
j=−∞
(
X−3,2j (e)
)2
sin δ2ωc− jn, (212)
where we have ignored the factor of χ .
Next we compute the total energy dissipation. For the simple driven oscillator the rate of
dissipation is given by Eq. (146). Generalizing this to our case, where there are many forcing
terms, and integrating over an orbit period, we find the average rate of energy dissipation to
be
d E
dt
= 3
4
M R2G21
ω221
∞∑
j=−∞
(
X−3,2j (e)
)2
(2ωc − jn) sin δ2ωc− jn (213)
+1
4
M R2G21
ω221
∞∑
j=−∞
(
X−3,0j (e)
)2
( jn) sin δ jn . (214)
For the special case of synchronous rotation we set ωc = n.
To derive the usual expression for tidal heating in synchronous rotation at low eccentricity,
we select the lowest order terms in eccentricity:
d E
dt
= 21
2
M R2G21
ω221
e2 sin δn, (215)
where δn = −δ−n = δ2ωc−n = −δ2ωc−3n . To get the familiar formula, we make a large μ
approximation (as before) and find
d E
dt
= 21
2
k2
Q
Gm2nR5
a6
e2, (216)
where we have again ignored a factor of χ = 0.98 and approximated sin δn by 1/Q. This
agrees with Peale and Cassen (1978) and Peale et al. (1979).
For the case of phase lags proportional to frequency we use δ2ωc− jn = (2 − j)δn and
δ jn = jδn . Making the large μ approximation again, we derive
d E
dt
= 21
2
k2
Q
Gm2nR5
a6
e2η(e), (217)
where the tidal heating enhancement factor is
η(e) = 1 + 18e2 + 3329
28
e4 + 55551
112
e6 + · · · . (218)
This agrees with Peale and Cassen (1978) and Wisdom (2008). Wisdom (2008) gave an
expression valid at arbitrary eccentricity:
ζ(e) = e2η(e) = 2
7
f0(e)
β15
− 4
7
f1(e)
β12
+ 2
7
f2(e)
β9
, (219)
123
J. Wisdom, J. Meyer
where
f0(e) = 1 + 31
2
e2 + 255
8
e4 + 185
16
e6 + 25
64
e8, (220)
f1(e) = 1 + 15
2
e2 + 45
8
e4 + 5
16
e6, (221)
f2(e) = 1 + 3e2 + 3
8
e4, (222)
with β = (1 − e2)1/2.
14.1 Orbital evolution not in spin-orbit resonance
Next, we can use conservation of energy and angular momentum with the above results to
compute the rates of change of semimajor axis and eccentricity for the case in which there
is no spin-orbit resonance. Based on conservation of energy,
da
dt
= − 2a
2
GMm
(
Cωcω˙c + d E
dt
)
. (223)
Substituting in Eqs. (211) and (214) yields
1
a
da
dt
= κ
2
∞∑
j=−∞
j
(
3(X−3,2j (e))
2 sin δ2ωc− jn − (X−3,0j (e))2 sin δ jn
)
, (224)
where
κ = a
GMm
MR2
G21
ω221
≈ k2nR
5
a5
m
M
, (225)
where the approximation is for large rigidity. Conservation of angular momentum gives
1
e
de
dt
= κ
4
√
1 − e2
e2
⎛
⎝−
√
1 − e2
∞∑
j=−∞
(X−3,0j (e))
2 j sin δ jn (226)
+ 3
∞∑
j=−∞
(
j
√
1 − e2 − 2
)
(X−3,2j (e))
2 sin δ2ωc− jn
⎞
⎠ . (227)
Note that the leading term in (1/e) de/dt as a polynomial in e is proportional to a constant—
the 1/e2 factor is cancelled by a leading factor of e2 in the subsequent expression.
If we keep only the lowest order terms in eccentricity, the expressions reduce to
1
a
da
dt
= 3κ sin δ2ωc−2n (228)
and
1
e
de
dt
= −3κ
4
(
3
2
sin δn + 1
4
sin δ2ωc−n + sin δ2ωc−2n − 49
4
sin δ2ωc−3n
)
, (229)
which confirms Equation 7 of Goldreich (1963). (Goldreich’s 3 is actually δn , which is
associated with frequency n, not 32n, as stated.)
123
Dynamic Elastic Tides
For Mignard tides, where phase lags are proportional to frequency,
1
a
da
dt
= 6 κ
Q
[(
ωc
n
− 1
)
+ e2
(
27
2
ωc
n
− 23
)
+ e4
(
573
8
ωc
n
− 180
)
+ e6
(
3961
16
ωc
n
− 6765
8
)
+ · · ·
]
(230)
1
e
de
dt
= 33
2
κ
Q
[(
ωc
n
− 18
11
)
+ e2
(
13
2
ωc
n
− 369
22
)
+ e4
(
181
8
ωc
n
− 3645
44
)
+ · · ·
]
. (231)
14.2 Orbital evolution in spin-orbit resonance
Similar calculations yield the rates of change of semimajor axis and eccentricity in spin-orbit
resonance, where ωc = kn and therefore
ω˙c = −3
2
k
n
a
da
dt
. (232)
For the case of a tide-raising perturber orbiting around an extended body,
1
a
da
dt
= 2a
GMm
(
3k2Γm − 1
)−1 d E
dt
, (233)
with Γm = C/(ma2). Substituting in Eq. (214) yields
1
a
da
dt
= κ
2
(
3k2Γm − 1
)−1
Υ, (234)
where
Υ =
∞∑
j=−∞
(
3(X−3,2j (e))
2(2k − j) sin δ2ωc− jn + (X−3,0j (e))2 j sin δ jn
)
. (235)
Conservation of angular momentum gives
1
e
de
dt
= −κ
4
√
1 − e2
e2
(√
1 − e2 − 3kΓm
1 − 3k2Γm
)
Υ (236)
Note that the first parenthesized factor in Eq. (236) is approximately one for small eccentricity
in synchronous rotation, in both the limit that Γm = C/(ma2) is large and small compared
to one. Both these limits occur in the Solar System: for the satellites of the giant planets
Γm 	 1, whereas for the Moon Γm  1.
For the case of an extended body orbiting around a perturbing central mass,
1
a
da
dt
= κ
2
(
3k2ΓM − 1
)−1
Υ, (237)
with ΓM = C/(Ma2), and
1
e
de
dt
= −κ
4
√
1 − e2
e2
(√
1 − e2 − 3kΓM
1 − 3k2ΓM
)
Υ (238)
Note that ΓM = C/(Ma2) is typically small compared to one.
123
J. Wisdom, J. Meyer
For Mignard tides we can expand Υ as a power series in e
Υ = 12(k − 1)2 + e2(90k2 − 324k + 276) + e4(315k2 − 1719k + 2160) (239)
+1
2
e6(1575k2 − 11883k + 20295) + · · · (240)
Notice that for synchronous rotation the leading term in Υ is order e2, but for other spin-orbit
commensurabilities the leading term is a constant. Thus a satellite would not be expected to
be in nonsynchronous spin-orbit resonance at low eccentricity. The eccentricity would decay
rapidly to zero in finite time and then the resonance lock would be lost.
If we specialize to synchronous rotation (k=1) and keep only the lowest order terms in
eccentricity, the expressions reduce to
1
a
da
dt
= 3κ
4
e2 (3Γm − 1)−1
(
3 sin δn + 1
2
sin δ2ωc−n − 49
2
sin δ2ωc−3n
)
(241)
for a perturber with mass m orbiting around a synchronous extended body with mass M and
1
a
da
dt
= 3κ
4
e2 (3ΓM − 1)−1
(
3 sin δn + 1
2
sin δ2ωc−n − 49
2
sin δ2ωc−3n
)
(242)
for an extended body with mass M orbiting around a central body with mass m. For Mignard
tides, where the phase lags are proportional to frequency, the common factor in parentheses
is 28δn . The rate of change of eccentricity is the same in both cases:
1
e
de
dt
= −3κ
8
(
3 sin δn + 1
2
sin δ2ωc−n − 49
2
sin δ2ωc−3n
)
, (243)
which again confirms Equation 8 of Goldreich (1963), if his 23 is actually the lag for
frequency n, not 32n. For Mignard tides, where the phase lags are proportional to frequency,
1
e
de
dt
= −21
2
κ
Q
. (244)
15 Tidal friction with eccentricity and obliquity
Next we consider a perturbing body in an elliptic orbit about a dissipative body with obliquity.
We continue to assume that there is no wobble (ωa = ωb = 0).
The equations of motion are
C ′ω˙c = √3MR2F1
(
Π1,1β
′γ ′ + Π1,22α′β ′ − Π1,−1α′γ ′ + Π1,−2((β ′)2 − (α′)2)
)
(245)
with
Π¨1,0 + d1Π˙1,0 + ω221Π1,0 = −2C1(ωc)2 + F1
2(γ ′)2 − (α′)2 − (β ′)2
2
(246)
Π¨1,1 + d1Π˙1,1 + ω221Π1,1 = F1
√
3α′γ ′ (247)
Π¨1,2 + d1Π˙1,2 + ω221Π1,2 = F1
√
3((α′)2 − (β ′)2)/2 (248)
Π¨1,−1 + d1Π˙1,−1 + ω221Π1,−1 = F1
√
3β ′γ ′ (249)
Π¨1,−2 + d1Π˙1,−2 + ω221Π1,−2 = F1
√
3α′β ′ (250)
We have again ignored a small periodic time dependence in C ′.
123
Dynamic Elastic Tides
The direction cosines are, in this approximation, o · s where
o = Rz(Ω)Rx (i)Rz(λ)(1, 0, 0), (251)
with λ = f + ω, the sum of the true anomaly and the argument of pericenter, and
s = Rz(Ω)Rx (i)Rz(Λ)Rx (I )Rz(φ)s0, (252)
and s0 is either (1, 0, 0), (0, 1, 0), or (0, 0, 1), for α′, β ′, and γ ′, respectively. Here I is the
obliquity of the body to the orbit, and Λ is the longitude of the equator minus the longitude
of the ascending node of the orbit, and φ = ωct . We find
α′ = g0 cos(λ′ − φ) + g1 cos(λ′ + φ) (253)
β ′ = g0 sin(λ′ − φ) − g1 sin(λ′ + φ) (254)
γ ′ = g2 sin(λ′), (255)
with λ′ = λ − Λ, and
g0 = (1 + cos I )/2 (256)
g1 = (1 − cos I )/2 (257)
g2 = − sin I. (258)
In terms of these, we can write down terms that enter the equations of motion
α′β ′ = g
2
0
2
sin(2λ′ − 2φ) − g
2
1
2
sin(2λ′ + 2φ) − g0g1 sin(2φ) (259)
α′γ ′ = g0g2
2
sin(2λ′ − φ) + g1g2
2
sin(2λ′ + φ)
+ g2(g0 − g1)
2
sin(φ) (260)
β ′γ ′ = −g0g2
2
cos(2λ′ − φ) + g1g2
2
cos(2λ′ + φ)
+ g2(g0 − g1)
2
cos(φ) (261)
(α′)2 − (β ′)2
2
= g
2
0
2
cos(2λ′ − 2φ) + g
2
1
2
cos(2λ′ + 2φ)
+ g0g1 cos(2φ) (262)
2(γ ′)2 − (α′)2 − (β ′)2
2
= g
2
2 − g20 − g21
2
− g
2
2 + 2g0g1
2
cos(2λ′). (263)
Proceeding as before, we Fourier expand the right-hand sides using the Hansen functions,
and find the solutions of the Π equations under the assumption that the orbital period is much
longer than the modal periods:
Π1,1 =
√
3
2
G1
ω221
∞∑
j=−∞
(
− X−3,2j (e)g0g2 sin(φ − jnt − 2(ω − Λ) − δωc− jn)
+ X−3,2j (e)g1g2 sin(φ + jnt + 2(ω − Λ) − δωc+ jn)
+ 1
2
X−3,0j (e)g2(g0 − g1) sin(φ + jnt − δωc+ jn)
+ 1
2
X−3,0j (e)g2(g0 − g1) sin(φ − jnt − δωc− jn)
)
(264)
123
J. Wisdom, J. Meyer
Π1,2 =
√
3
2
G1
ω221
∞∑
j=−∞
(
X−3,2j (e)g
2
0 cos(2φ − jnt − 2(ω − Λ) − δ2ωc− jn)
+X−3,2j (e)g21 cos(2φ + jnt + 2(ω − Λ) − δ2ωc+ jn)
+ 1
2
X−3,0j (e)2g0g1 cos(2φ + jnt − δ2ωc+ jn)
+ 1
2
X−3,0j (e)2g0g1 cos(2φ − jnt − δ2ωc− jn)
)
(265)
Π1,−1 =
√
3
2
G1
ω221
∞∑
j=−∞
(
−X−3,2j (e)g0g2 cos(φ − jnt − 2(ω − Λ) − δωc− jn)
+ X−3,2j (e)g1g2 cos(φ + jnt + 2(ω − Λ) − δωc+ jn)
+ 1
2
X−3,0j (e)g2(g0 − g1) cos(φ + jnt − δωc+ jn)
+ 1
2
X−3,0j (e)g2(g0 − g1) cos(φ − jnt − δωc− jn)
)
(266)
Π1,−2 =
√
3
2
G1
ω221
∞∑
j=−∞
(
−X−3,2j (e)g20 sin(2φ − jnt − 2(ω − Λ) − δ2ωc− jn)
− X−3,2j (e)g21 sin(2φ + jnt + 2(ω − Λ) − δ2ωc+ jn)
− 1
2
X−3,0j (e)2g0g1 sin(2φ + jnt − δ2ωc+ jn)
− 1
2
X−3,0j (e)2g0g1 sin(2φ − jnt − δ2ωc− jn)
)
(267)
Π1,0 = −2C1 (ω
c)2
ω221
+ 1
2
G1
ω221
∞∑
j=−∞
(
X−3,0j (e)(g
2
2 − g20 − g21) cos( jnt − δ jn)
− X−3,2j (e)(g22 + 2g0g1) cos( jnt + 2(ω − Λ) − δ jn)
)
(268)
Generalizing Eq. (146), we find the rate of energy dissipation to be
d E
dt
= M R
2G21
ω221
∞∑
j=−∞
[
F0(X
−3,2
j (e))
2 + F1(X−3,0j (e))2 (269)
+ F2X−3,2j (e)X−3,0j (e) cos(2(ω − Λ)) (270)
+ F3X−3,2j (e)X−3,2− j (e) cos(4(ω − Λ))
]
(271)
with
F0 = 3
4
[
(g0g2)
2(ωc − jn) sin δωc− jn + (g1g2)2(ωc + jn) sin δωc+ jn
+ (g0g0)2(2ωc − jn) sin δ2ωc− jn + (g1g1)2(2ωc + jn) sin δ2ωc+ jn
]
+1
8
(g22 + 2g0g1)2 jnδ jn (272)
123
Dynamic Elastic Tides
F1 = 3
8
[
(g2(g0 − g1))2((ωc − jn) sin δωc− jn + (ωc + jn) sin δωc+ jn)
+ (2g0g1)2((2ωc − jn) sin δ2ωc− jn + (2ωc + jn) sin δ2ωc+ jn)
]
+ 1
4
(g22 − g20 − g21)2 jnδ jn (273)
F2 = 3
2
[
g22(g1 − g0)(g0(ωc − jn) sin δωc− jn − g1(ωc + jn) sin δωc+ jn)
+ 2g30g1(2ωc − jn) sin δ2ωc− jn + 2g0g31(2ωc + jn) sin δ2ωc+ jn
]
− 1
2
(g22 − g20 − g21)(g22 + 2g0g1) jn sin δ jn (274)
F3 = 3
64
g42
[−4((ωc − jn) sin δωc− jn + (ωc + jn) sin δωc+ jn)
+ ((2ωc − jn) sin δ2ωc− jn + (2ωc + jn) sin δ2ωc+ jn)
]
+ 9
32
g42 jn sin δ jn . (275)
For the special case of synchronous rotation we set ωc = n. To derive the familiar case, we
set δkωc− jn = (k− j)δ, and then approximate sin δ by1/Q.Making the largeμ approximation
as before, we derive
d E
dt
= 21
2
k2
Q
Gm2nR5
a6
ζ(e, I ), (276)
where
ζ(e, I ) = 2
7
f¯0(e) − 4
7
f¯1(e) cos I + 1
7
f¯2(e)(1 + (cos I )2) (277)
+ 3
28
e2 f¯3(e)(sin I )
2 cos(2(ω − Λ)). (278)
with
f¯0(e) = 1 + 23e2 + 180e4 + 6765
8
e6 + · · · (279)
f¯1(e) = 1 + 27
2
e2 + 573
8
e4 + 3961
16
e6 + · · · (280)
f¯2(e) = 1 + 15
2
e2 + 105
4
e4 + 525
8
e6 + · · · (281)
f¯3(e) = 1 + 14
3
e2 + 105
8
e4 + 231
8
e6 + · · · . (282)
This agrees with Wisdom [2008], where these series are given a finite representation:
f¯0(e) = f0(e)
β15
, (283)
f¯1(e) = f1(e)
β12
, (284)
f¯2(e) = f2(e)
β9
, (285)
f¯3(e) = f3(e)
β13
, (286)
123
J. Wisdom, J. Meyer
where
f3(e) = 1 − 11
6
e2 + 2
3
e4 + 1
6
e6. (287)
16 Conclusion
In this paper we have presented a new framework for understanding and calculating the
effects of tidal friction on solid bodies. Our approach is to consider the body as a dynamical
system described by its modes of oscillation. The dynamics of the modes are determined by
the internal distribution of matter and the rheology. Dissipation occurs in the dynamics of the
modes. Herewe consider only a simplemodel for dissipation: we assume themodal dynamics
is governed by a damped harmonic oscillator-like equation of motion. We derive closed
differential equations that govern the tidal process. This is to be contrasted with the usual
Darwinian approach in which tides are first raised, then given a phase lag before considering
their effects on other bodies. We are able to rederive the Darwin result as an approximate
solution to our differential equations. Note that we need not make this approximation—we
can directly integrate the differential equations governing the tidal process. We rederive many
of the standard results of tidal theory directly in our new framework. We plan to extend these
results to different rheologies in a future work.
Appendix: Surface and solid spherical harmonics
We will use the terms “solid spherical harmonic” and “surface spherical harmonic.” The
normalized surface spherical harmonics are
Cml (θ, φ) = Pml (cos θ) cos(mφ) (288)
Sml (θ, φ) = Pml (cos θ) sin(mφ), (289)
where the normalized associated Legendre polynomials satisfy
Pml (x) = Nml Plm(x), (290)
with
Nml =
[
(2 − δm,0)(2l + 1) (l − m)!
(l + m)!
]1/2
, (291)
and
Plm(x) = 1
2l l! (1 − x
2)m/2
dl+m
dxl+m
[
(x2 − 1)l
]
. (292)
A few of the Plm are
P20(cos θ) = 3
2
(cos θ)2 − 1
2
(293)
P21(cos θ) = 3 sin θ cos θ (294)
P22(cos θ) = 3(sin θ)2. (295)
For convenience we introduce
Xml (θ, φ) = Cml (θ, φ) m ≥ 0 (296)
= S−ml (θ, φ) m < 0 (297)
123
Dynamic Elastic Tides
The Xml (C
m
l and S
m
l ) are orthonormal in that
δll ′δmm′ = 14π
∫ π
0
∫ 2π
0
Xml (θ, φ)X
m′
l ′ (θ, φ) sin θdθdφ (298)
The solid spherical harmonics are
X˜ml
( x
R
,
y
R
,
z
R
)
=
( r
R
)l
Xml (θ, φ), (299)
where r2 = x2 + y2 + z2, cos θ = z/r , φ = atan(y, x). The solid harmonics satisfy the
orthogonality relations
δll ′δmm′
2l + 3 (300)
= 1
4π R3
∫ a
0
∫ π
0
∫ 2π
0
X˜ml
( x
R
,
y
R
,
z
R
)
X˜m
′
l ′
( x
R
,
y
R
,
z
R
)
r2 sin θdrdθdφ (301)
= 1
4πa3
∫
V
X˜ml X˜
m′
l ′ dV, (302)
where the V is the sphere of radius R, and the last line introduces an abbreviated notation.Note
that whenever we write a X˜ml without arguments, they are assumed to be (x/R, y/R, z/R).
It is convenient to introduce
X¯ml =
1√
2l + 1 X˜
m
l . (303)
The first few X¯ml are
X¯00(x, y, z) = 1 (304)
X¯01(x, y, z) = z (305)
X¯11(x, y, z) = x (306)
X¯−11 (x, y, z) = y (307)
X¯02(x, y, z) = (2z2 − x2 − y2)/2 (308)
X¯12(x, y, z) =
√
3xz (309)
X¯22(x, y, z) =
√
3(x2 − y2)/2 (310)
X¯−12 (x, y, z) =
√
3yz (311)
X¯−22 (x, y, z) =
√
3xy. (312)
Let X˜ml be a solid spherical harmonic. It is a homogeneous function of degree l in the
coordinates (x/R, y/R, z/R). Euler’s theorem tells us that
(x · ∇) X˜ml (x/R, y/R, z/R) = l X˜ml (x/R, y/R, z/R). (313)
A little calculation shows that
∇2(xX˜ml (x/R, y/R, z/R)) = 2∇ X˜ml (x/R, y/R, z/R), (314)
and
∇2(rk X˜ml (x/R, y/R, z/R)) = k(k + 2l + 1)rk−2 X˜ml (x/R, y/R, z/R). (315)
123
J. Wisdom, J. Meyer
References
Chandler, S.C.: On the variation of latitude. Astron. J. 11, 59–61 (1891)
Darwin, G.: Scientific Papers, vol. 2. Cambridge University Press, Cambridge (1908)
Dobrovolskis, A.: Spin states and climates of eccentric exoplanets. Icarus 192, 1–23 (2007)
Eggleton, P., Kiseleva, L., Hut, P.: The equilibrium tide model for tidal friction. Astrophys. J. 499, 853–870
(1998)
Goldreich, P.: On the eccentricity of satellite orbits in the solar system. MNRAS 126, 257–268 (1963)
Goldreich, P.: History of the lunar orbit. Rev. Geophys. Space Phys. 4, 411–439 (1966)
Goldreich, P., Peale, S.J.: Spin-orbit coupling in the solar system. Astron. J. 71, 425–438 (1966)
Hut, P.: Tidal evolution in close binary systems. Astron. Astrophys. 99, 126–140 (1981)
Jefferys, H.: The Earth, 6th edn. Cambridge University Press, Cambridge (1976)
Kaula, W.M.: Tidal dissipation by solid friction and the resulting orbital evolution. Rev. Geophys. 2, 661–685
(1966)
Lamb, H.: On the vibrations of an elastic sphere. Proc. Lond. Math Soc. 13, 189–212 (1882)
Love, A.E.H.: A Treatise on the Mathematical Theory of Elasticity, 4th edn, pp. 257–259. Dover, New York,
NY (1944)
Mignard, F.: The lunar orbit revisited, III. Moon Planets 24, 189–207 (1981)
Newcomb, S.: On the dynamics of the earth’s rotation, with respect to the periodic variations of latitude.
MNRAS 52, 336–341 (1882)
Peale, S.J.: Rotation of solid bodies in the solar system. Rev. Geophys. Space Phys. 11, 767–793 (1973)
Peale, S.J., Cassen, P.: Contribution of tidal dissipation to lunar thermal history. Icarus 36, 245–269 (1978)
Peale, S.J., Cassen, P., Reynolds, R.: Melting of Io by tidal dissipation. Science 203, 892–894 (1979)
Peale, S.J.: Resonances. In: Burns, J.A., Matthews, M.S. (eds.) Satellites, pp. 159–223. University of Arizona
Press, Tucson (1986)
Poincaré, H.: Sur une forme novelle des équations de la mécanique. Comptes rendus de l’Académie des
Sciences 132, 369–371 (1901)
Stacey, F.D.: Physics of the Earth, 3rd edn. Brookfield Press, Brisbane (1992)
Williams, J.G.: Lunar core and mantle. What does LLR see? Proceedings of the 16th International Workshop
on Laser Ranging, pp. 101–120 (2008)
Wisdom, J.: Tidal heating at arbitrary eccentricity and obliquity. Icarus 193, 637–640 (2008)
Zahn, J.-P.: The dynamical tide in close binaries. Astron. Astrophys. 41, 329–344 (1975)
123