Spin-orbit resonances: selection by dissipation

F. Bonetto^{a} , G. Gallavotti^{b} and G. Gentile^{c}

a Mathematics Department, Georgia Institute of Technology, Atlanta, GA 30332 b INFN & Fisica, Università di Roma “La Sapienza”, I-00185 c Dipartimento di Matematica, Università di Roma Tre, Roma, I-00146

Abstract. We consider a class of quasi-integrable Hamiltonian systems with “one and a half degrees of freedom” and study how friction forces, analytically depending on the perturbation parameter, can stabilize some particular periodic orbits. As a general scenario most periodic orbits can, in absence of friction, be continued under perturbation but disappear when the system becomes dissipative: given a friction value only a finite, system-dependent, set among them is left. Applications to Celestial Mechanics models of spin-orbit type are discussed in connection with resonance-locking between revolution and rotation periods of satellites.

1. Introduction sec.1

Many planetary motions are approximately periodic: for example the revolution of planets around the Sun, of satellites around planets or the rotation of planets and satellites around their axis. In some cases the periods are linked by a very simple rational relation. A well known example is provided by the large satellites of the major planets: in almost all cases rotation and revolution periods are equal. A deviation from this rule is Mercury for which the rotation period is two thirds of the revolution period (i.e. for Mercury there are 3 days in 2 years). The latter very interesting phenomenon is not easy to explain from a Hamiltonian point of view: indeed in a Hamiltonian system close to an integrable one essentially all periods corresponding to unperturbed orbits can occur even in the perturbed systems, independently of the period.

To simplify one normally assumes that the energy contained in the revolution is much larger than that contained in the rotation, leading to assume that the satellite moves on a fixed Keplerian orbit. Several simple mathematical models describing such a situation are called spin-orbit models, see for instance [10] and [5]. The models are periodically forced one degree of freedom Hamiltonian systems, i.e. “one and a half degrees of freedom” systems.

We consider several quasi–integrable Hamiltonian systems with one and a half degrees of freedom; in absence of friction the system will present a variety of motions among which, if the perturbation is not too large, isolated periodic orbits remnants of “broken-down” resonant tori.^{1} In general for every given period we expect an even number of such orbits, half of which are elliptic and half hyperbolic. Introducing a small friction term in the equations of motion, some elliptic periodic orbits of the considered models become asymptotically stable and acquire a non-empty basin of attraction. It is an important fact that for a given positive friction only a few of them remain stable or just existing: it seems likely that, in interesting cases (like the ones we consider here), the union of the basins of attraction of the stable orbits is the full phase space up to a set of zero volume. This seems confirmed by numerical simulations in similar cases; see [4].

1 That is periodic orbits that can be followed continuously as the perturbation parameter ε grows from 0 to (small)

positive values.

One can imagine that the system is subjected to an initial friction that “eventually” becomes negligible (but after a time scale larger than the characteristic periods of the system) because of a change in the state of the system (e.g. from fluid to solid): then in the long time limit the system will be found near one of the above mentioned stable periodic orbits.

We shall concentrate on three particular cases. The first is the pendulum with periodically driven point of support. This system was numerically studied in [4]: here it will be considered for small values of the parameters, in order to study it analytically by perturbation theory. The analytical results that we find are indeed inspired by the quoted numerical analysis.

The second model is a special case of a gyroscopic motion, periodically forced. And the third case is the above mentioned spin–orbit model, more closely related to Celestial Mechanics.

The latter model is classical and very simple: it assumes the equations of motion that would be obeyed by an asymmetric rigid body constrained to rotate around an axis which moves on a Keplerian elliptic orbit orthogonally to the orbit plane. In Appendix A4 we discuss several theoretical questions pertaining the approximation by studying the equations that would describe the motion of a rigid body with the center of mass moving on a Keplerian orbit but free to change the orientation of the spin and symmetry axes. The point of our analysis has been to check in which cases the equatios of motion are analytic in a proper system of coordinates so that neglecting higher order correctons might be justfied by a suitable perturbation analysis (which we do not discuss, see [3]).

In the first two cases, considered mainly for illustration purposes, we assume a mathematically simple friction model imagined to be due to an external background. Then we extend the results to friction forces which might be more reasonable for the spin-orbit model: friction is imagined as due to tidal dissipation on a fluid planet whose orbital motion occurs on a fixed orbit.

Further numerical study, like the one performed in [4], of the basins of attraction of the periodic orbits would be highly desirable in order to detect which periodic orbits are really the natural candidates for attracting the majority of confined motions when at least two of them coexist.

The mechanism of capture into resonance arising in systems differing from integrable ones by a small perturbation has been studied by several authors starting with the theory of capture into the 3:2 resonance of Mercury [10], [11]. The general mechanism is discussed and summarized in the review article [15]. In the latter paper friction is considered either periodic or just not depending on time. Here we regard the friction as not periodic in time and, ideally, abruptly changing order of magnitude from a small value to a neglegible value: a situation that we consider possible in the formation of a planet. At the beginning, when the planet can be considered in a fluid state the dissipation (due to tidal effects on an ellipsoidic rotating fluid) is sensible (though small), while after a suitably large (but astronomically not so large) time it becomes negligible: we call this time the “solidification time”.

Therefore the problem is: which are the possible stable orbits at the solidification time? We know from [10], [11], [15] that the system will settle into a periodic orbit with a probabilistic pattern. However the orbit on which, randomly, the system will settle has to be among the ones that are stable at the solidification epoch. Once the orbits are known a well developed theory [15] will even allow us to estimate the probability that a randomly chosen initial datum will be asymptotic to a given periodic orbit (among the existing ones). We study here possible criteria to determine periodic orbits and to select the (very few) ones that can compete in the random selection of the “final” periodic orbit.

In the Celestial Mechanics cases our model is oversimplified, and more realistic pictures could be devised [12]; nevertheless, because of its simplicity, it is suitable for analytical investigations (as opposed to just numerical ones) on the relevance of friction in the early stages of evolution of heavenly bodies and for the selection of structurally stable periodic motions. An important role is plaid by the assumption on the size of the friction and on the actual existence of stable periodic motions (which has to be checked). Moreover we interpret the extension studied in Appendix A5, and supplementing the general analysis in [15], as checking that, at least at low orders in perturbation theory, only few qualitative properties of friction determine which orbit survives.

2. Statement of the analytical results sec.2

p.2.1 2.1. The models. Consider the equation

θ̈ + εG(θ, t) + γ(θ̇ − μ) = 0, (2.1) 2.1

where θ ∈ 𝕋 = 2π/ℤ, the function G(θ, t) is 2π-periodic and analytic in each variable and ε, γ ∈ ℝ, with γ > 0: we shall call ε the perturbation parameter and γ the friction constant; the parameter μ will be either 0 or 1. The case μ = 0 will be called the background friction model, while μ = 1 will be called the tidal friction model.

For γ = 0 the equation (2.1) is a typical equation that arises in several (Celestial) Mechanics problems; see [17], [14] and [13]. In this case we can derive (2.1) as the Hamilton equations of the system described by the Hamiltonian

H(θ, Θ, t, T) =^{def}ωΘ + \frac{1}{2}Θ^{2} + T + εg(θ, t), (2.2) 2.2

where (θ, Θ) ∈ 𝕋 × ℝ and (t, T) ∈ 𝕋 × ℝ are conjugated variables, ω ∈ ℝ is a parameter, and ∂_{θ}g(θ, t) = G(θ, t).

For ε = 0 the system (2.2) admits the one parameter family of solutions X(t) = (θ(t), t, Θ(t), T (t)) = (θ_{0} + ωt, t, 0, 0), with θ_{0} ∈ [0, 2π). If ω = p/q ∈ ℚ, with p, q relatively prime, each such solution is periodic and it is convenient to use coordinates in which it appears particularly simple. For this purpose the standard procedure is to define a (canonical) linear change of variables

( \begin{matrix}α \\ β \end{matrix} ) = C ( \begin{matrix}θ \\ t \end{matrix} ) ,  ( \begin{matrix}A \\ B \end{matrix} ) = (C^{T} )^{−1} ( \begin{matrix}Θ \\ T \end{matrix} ) ,  (2.3) 2.3

where

C = ( \begin{matrix}m & −n \\ −q & p \end{matrix} ) ,  C^{−1} = ( \begin{matrix}p & n \\ q & m \end{matrix} ) ,  (2.4) 2.4

with (m, n) ∈ ℤ^{2} such that mp − nq = det C = 1.

Then, by setting ω = (ω, 1), one has ω_{0} ≡ Cω = (1/q, 0), and the equations of motion become, if for instance μ = 0,

\{ \begin{matrix}α̇ = 1/q + m (mA − qB) , \\ β̇ = −q (mA − qB) , \end{matrix}  \{ \begin{matrix}Ȧ = −ε∂_{α}f(α, β) − pγ (mA − qB) − p^{2}γ/q, \\ Ḃ = −ε∂_{β}f(α, β) − nγ (mA − qB) − npγ/q, \end{matrix}  (2.5) 2.6

where f(α, β) = g(θ(α, β), t(α, β)). For γ = 0 the corresponding Hamiltonian is

H = \frac{1}{2}m^{2}A^{2} + \frac{1}{2}q^{2}B^{2} − mqAB + \frac{A}{q} + εf(α, β), (2.6) 2.7

and, for ε = 0, the one parameter family of periodic solutions X(t) considered above is transformed into X_{0}(t) = CX(t) = (α_{0} + t/q, β_{0}, 0, 0), with qα_{0} + mβ_{0} = 0: this makes clear the main feature of the resonance (i.e. only one angle is really rotating while the motion of the other will be entirely controlled by the perturbation).

p.2.2 2.2. Persistence of periodic solutions. If γ ≠ 0 there can be no periodic solution close to X_{0}(t) for the system (2.5). Indeed due to the presence of a non-vanishing friction the energy of the system decreases after every period (in fact, for instance in the case μ = 0 for the Hamiltonian H defined in (2.2), equations (2.1) imply Ḣ = −γθ̇^{2}). However, as usual in the study of forced system, we shall call periodic a solution such that the (α, β) variables (or equivalently the (t, θ) variables) are periodic in time. And instead of (2.5) or the corresponding one for the case μ = 1 we consider the (equivalent) system

\{ \begin{matrix}α̈ = −mε (m∂_{α}f − q∂_{β}f) − γα̇ − nγ + μmγ, \\ β̈ = qε (m∂_{α}f − q∂_{β}f) − γβ̇ + pγ − μqγ, \end{matrix} (2.7) 2.8

where μ = 0, 1, and look for the existence of a periodic solution merging as ε, γ → 0 with the motion t → (α_{0} + t/q, β_{0}) (obtained by looking at the angle variables of X_{0}(t)) when ε, γ ≠ 0.

It will be convenient to write

f(α, β) = ∑_{ν∈ℤ} e^{iνα}f_{ν}(β), (2.8) 2.9

where, for all ν ∈ ℤ, the coefficient f_{ν}(β) is a 2π-periodic function of β.

The form of (2.7) shows that if γ > 0 and ε = 0 there is no periodic solution close to the unperturbed one. We nevertheless expect a periodic solution to exist as ε, γ vary within a set of parameters values containing a cone of the form C_{−}ε < γ < C_{+}ε. For this reason we shall fix γ = Cε, with C a parameter to be varied. We could also consider, more generally, that γ is an analytic function of ε (divisible by ε), and the discussion to which the rest of the paper is devoted could be extended to cover such a case; however we prefer do not overwhelm the analysis with unessential technical intricacies that wider generality would inevitably generate, see Appendix A5. From a physical viewpoint one should imagine that, in concrete examples, the friction parameter γ is fixed to some value, then we could write it as γ = Cε which, for given ε, fixes C to some numerical value and one should then check that the value of C fulfills the conditions that we shall find.

We shall look for a solution which is analytic in ε for ε small enough. This means that we shall write

α(t) = α_{0} + ω_{0}t + a(α_{0} + ω_{0}t, β_{0}; ε),  β(t) = β_{0} + b(α_{0} + ω_{0}t, β_{0}; ε),  ω_{0} = \frac{1}{q},  (2.9) 2.10

where a(ψ, β; ε) and b(ψ, β; ε) will be expanded as

a(ψ, β; ε) = ∑_{k=1}^{∞} ε^{k} ∑_{ν∈ℤ} e^{iνψ}a_{ν}^{(k)}(β),  b(ψ, β; ε) = ∑_{k=1}^{∞} ε^{k} ∑_{ν∈ℤ} e^{iνψ}b_{ν}^{(k)}(β).  (2.10) 2.11

If solutions of the form (2.9) exist then t ≡ qα(t) + mβ(t), so that one must have qα_{0} + mβ_{0} = 0, hence α_{0} has to be fixed as α_{0} = −mβ_{0}/q, while b(α_{0}, β_{0}; ε) will be suitably fixed and a(α_{0}, β_{0}; ε) will be so chosen to obtain qa(α_{0}, β_{0}; ε) + mb(α_{0}, β_{0}; ε) = 0. We shall prove in Section 3 the following result.

p.2.3 2.3. Theorem. Fix ω = p/q and fix C so that

min_{β∈[0,2π]} ∂_{β}f_{0}(β) < C \frac{p − μq}{q^{2}} < max_{β∈[0,2π]} ∂_{β}f_{0}(β). (2.11) 2.12

If there exists β_{0} such that

∂_{β}f_{0}(β_{0}) = \frac{p − μq}{q^{2}} C,  ∂_{β}^{2}f_{0}(β_{0}) ≠ 0,  (2.12) 2.13

then for ε small enough there is a periodic solution of the equations of motion of the form (2.9).

p.2.4 2.4. Remarks. (1) The condition (2.12) can be satisfied only if the function f_{0}(β) is not identically

constant. For instance, if g(θ, t) = − (1 + cos t) cos θ, (2.13) 2.14

one has

f(α, β) = − cos (pα + nβ) − \frac{1}{2} [cos ((p + q)α + (n + m)β) + cos ((p − q)α + (n − m)β)] , (2.14) 2.15

so that f_{0}(β) ≡ 0 except for p = q = 1 (that is ω = 1). On the other hand, if g is an analytic function, the condition that the function f_{0}(β) is not identically constant is generic.

(2) If the stationary points of the function ∂_{β}f_{0}(β) correspond only to either maxima or minima, then as the value of C increases, the two conditions (2.12) fail to be satisfied simultaneously.

(3) In the above discussion the dependence on C of the solutions is not explicitly indicated. However the proof in Section 3 implies that the solutions (2.9) are analytic in C near any C satisfying (2.12). We can therefore reformulate our theorem by stating that the solutions (2.9) are analytic in ε and γ in the intersection of a neighborhood of the origin with the cone C_{−}ε < γ < C_{+}ε, where C_{−} and C_{+} are determined by (2.11).

(4) It is important to realize that the above theorem does not provide a one parameter family of solutions: since only a finite number of values for β_{0} will, in general, be allowed by (2.12) we can say that (in general) only a finite number among the free solutions, i.e. the solutions existing when ε = 0, can be continued to ε > 0.

p.2.5 2.5. Extensions. I. The case in which ∂_{β}f_{0}(β) vanishes identically is excluded from the above analysis: such a case can be dealt with by setting γ = Cε^{2}. In the case μ = 0, for instance, we have the following result (also proved in Section 3).

p.2.6 2.6. Theorem. Let μ = 0. Fix ω = p/q such that ∂_{β}f_{0}(β) ≡ 0, and, setting D = (ω_{0}∂_{ψ})^{2}, define^{2}

\begin{aligned}F_{0}^{(2)}(β) = − & [\frac{1}{2}m^{2}∂_{α}f(β)D^{−2}∂_{α}f(β) + \\ & + \frac{1}{2}q^{2}∂_{β}f(β)D^{−2}∂_{β}f(β) − mq∂_{α}f(β)D^{−2}∂_{β}f(β)]_{0} . \end{aligned} (2.15) 2.16

Fix C such that

min_{β∈[0,2π]} F_{0}^{(2)}(β) < C \frac{p}{q^{2}} < max_{β∈[0,2π]} F_{0}^{(2)}(β). (2.16) 2.17

2 Given a function F (α, β) we denote by [F (β)]_{ν} its νth Fourier coefficient with respect to α.

If there exists β_{0} such that

F^{(2)}(β_{0}) = \frac{p}{q^{2}} C,  ∂_{β}F^{(2)}(β_{0}) ≠ 0,  (2.17) 2.18

then for ε small enough there is a periodic solution of the equations of motion of the form (2.9).

p.2.7 2.7. Remarks. (1) The condition (2.16) is a second order condition, as the analysis of Section 3.2 will show, while condition (2.11) of theorem 2.3 was a first order condition.

(2) Condition (2.15) is generic even if we restrict the analysis to trigonometric polynomials g; in particular, for g given by (2.13), the condition (2.16) is satisfied for C ≠ 0.

p.2.8 2.8. Extensions. II. We can also consider the case γ = Cε^{2} when the function f_{0}(β) is not identically constant. In such a case the following result holds: again we consider only the case μ = 0 for simplicity.

p.2.9 2.9. Theorem. Let μ = 0. Fix ω = p/q and fix β_{0} such that one has

∂_{β}f_{0}(β_{0}) = 0,  ∂_{β}^{2}f_{0}(β_{0}) ≠ 0.  (2.18) 2.19

If γ = Cε^{2} then for ε small enough there is a periodic solution of the equations of motion of the form (2.9) provided that one has |C| < C_{0} for some positive constant C_{0}.

p.2.10 2.10. Remark. More generally we can set γ = Cε^{n}, with n ≥ 2, and a result like theorem 2.9 holds. The new solution is really different from the one in theorem (2.3) as it will be seen from the fact that they differ already to first order in ε: this means that there are several families of periodic orbits which merge continuously as ε → 0 with the unperturbed orbit. In physical situations the parameters ε and γ are fixed: then for γ small enough one can have distinct periodic orbits corresponding to several values of n. However the smallness condition that we find on ε becomes more and more stringent as n increases.

p.2.11 2.11. Contents. In Section 3 we briefly discuss the proofs of the theorems stated above referring to Appendix A1 for the more technical aspects; we stress that in the present case the analysis appears much easier because of the absence of small divisors (as we are interested in periodic rather than in quasi-periodic solutions).

The periodic orbits appear in pairs of stable and unstable orbits: this is a consequence of Poincaré –Birkhoff’s theorem [1]. In Appendix A2 we briefly study the stability of the periodic solutions by performing a low orders analysis.

In Section 4, we consider a very simple form for the friction term with μ = 0 and study some elementary applications: the periodically driven pendulum and a spin-orbit model for a rigid body both in a background friction case. Finally in Section 5 we shall realize that, also by considering a tidal friction case (and possibly a more general form for the friction term), the scenario remains essentially unchanged. In particular we discuss some qualitative properties of systems relevant in Celestial Mechanics: the systems satellite-planet (like Moon-Earth) and the systems planet-star (like Mercury-Sun).

p.2.12 2.12. Conclusions. The general conclusions that we can draw about the resonance-locking in the spin-orbit problem are the following. In absence of friction each family of periodic orbits present in the unperturbed system also occurs, in general, after perturbation in the sense that at least one of the unperturbed orbits can be continued at ε > 0, independently of the value of the period if ε is small enough (how small depends on the unperturbed orbit considered, hence in particular on its period). This is a general property of Hamiltonian systems (well known: for a derivation with the methods used here see [9]).

If the friction is different from zero only a finite number of periodic orbits can be continued in ε up to a prefixed value ε_{0}. If the friction constant is large enough (but still not too large, see (2.11)) then either there is only one periodic orbit or a few of them which are accessible by continuing at ε > 0 some unperturbed ones.

We imagine that, in the history of the planetary system, friction decreased slowly: we first realize that the resonance 1:1 is stable since the beginning. As friction decreased while the planet was spinning down toward the 1:1 resonance, one after the other, other resonances^{3} became stable but in most cases (basically all except the case Mercury-Sun) the 1:1 was the only one to exist for a long time, long enough that the system could essentially fall so close to the resonance to become unaffected by the new possibilities open by the evolution into stability of the other resonances. This means that the motions were already well inside the basin of attraction of the 1:1 resonance before any other stable periodic orbit could exist: this could explain why the resonance 1:1 is almost always the dominant one. In the Mercury-Sun case as well as in all other cases considered the resonance 3:2 is the first to become stable as the friction decreases and therefore it is the most likely one, after the 1:1 resonance, to stabilize some nearby trajectory, a case that seems to have happened in the system Mercury-Sun which is locked in a resonance 3:2. We shall also find that the capture into the resonance 3:2 is more likely to occur in the case of Mercury-Sun rather than in the other cases: Mercury-Sun is the case in which such a resonance seems to appear earlier by far, essentially because of the larger value of the eccentricity. In particular this yields that Mercury-Sun is essentially the only case in which such a phenomenon could be really expected to happen, as in fact it happened.

A more precise discussion should include an estimate of the friction and of the time scales involved together with an analysis of the sizes of the basins of attraction of the different periodic orbits as functions of the friction. This goes beyond the scope of the present paper.

3. Proof of the theorems sec.3

p.3.1 3.1. Proof of theorem 2.3. The analysis follows the classical perturbation theory pattern for analytic perturbations: we first check that the problem is soluble to all orders of pertubation theory so that a power series solution can be defined up to convergence analysis and subsequently convergence is checked for small enough perturbaton parameter values. Convergence is a rather standard check once the expansion coefficients have been shown to exist and be algorithmically computable. In this section we derive the power series and the uneventful convergence check is in Appendix A1.

For γ = Cε, inserting (2.10) into the equations of motion (2.7) gives, for k = 1,

\begin{aligned}(iω_{0}ν)^{2} a_{ν}^{(1)} & = −m (imνf_{ν} − q∂_{β}f_{ν}) − C \frac{m(p − μq)}{q}δ_{ν,0}, \\ (iω_{0}ν)^{2} b_{ν}^{(1)} & = q (imνf_{ν} − q∂_{β}f_{ν}) + (p − μq)Cδ_{ν,0}, \end{aligned} (3.1) 3.1

and, for k ≥ 2,

\begin{aligned}(iω_{0}ν)^{2} a_{ν}^{(k)} & = −m [m∂_{α}f − q∂_{β}f]_{ν}^{(k−1)} − C (iω_{0}ν) a_{ν}^{(k−1)}, \\ (iω_{0}ν)^{2} b_{ν}^{(k)} & = q [m∂_{α}f − q∂_{β}f]_{ν}^{(k−1)} − C (iω_{0}ν) b_{ν}^{(k−1)}, \end{aligned} (3.2) 3.2

3 As we shall see the onset of stability can be related, or at least bounded, in terms of the size of q if the period is T = 2πq/p. The higher q, at p/q = ω roughly constant, the smaller the value of the viscosity below which the orbit exists and is stable.

where, given any function F admitting a formal series

F (ψ; ε) = ∑_{k=1}^{∞} ε^{k} ∑_{ν∈ℤ} e^{iνψ}F_{ν}^{(k)}, (3.3) 3.3

[F]_{ν}^{(k)} denotes the coefficient with Taylor label k and Fourier label ν; see footnote 2. The above equations can be solved for ν ≠ 0 provided that, for ν = 0, one has, for k = 1,

0 = mq∂_{β}f_{0}(β_{0}) − C \frac{m(p − μq)}{q},  0 = −q^{2}∂_{β}f_{0}(β_{0}) + (p − μq)C,  (3.4) 3.4

and, to order k ≥ 2,

0 = −m [m∂_{α}f − q∂_{β}f]_{0}^{(k−1)} ,  0 = q [m∂_{α}f − q∂_{β}f]_{0}^{(k−1)} .  (3.5) 3.5

The two equations (3.4) correspond to the single equation

∂_{β}f_{0}(β_{0}) = C \frac{p − μq}{q^{2}} , (3.6) 3.6

as it is immediate to check. Also the two equations (3.5) reduce to one equation,

[m∂_{α}f − q∂_{β}f]_{0}^{(k−1)} = 0, (3.7) 3.7

which can be solved by suitably fixing the sequence {b_{0}^{(k)}}_{k∈ℕ}: the sequence {a_{0}^{(k)}}_{k∈ℕ} is, finally, determined by imposing the validity of the exact relation α(0)q + β(0)m ≡ 0, see Appendix A1. More precisely [m∂_{α}f − q∂_{β}f]_{0}^{(k)} can be made equal to zero for all k by a suitable choice of the sequence {b_{0}^{(k)}}_{k∈ℕ}: this can be easily done if ∂_{β}^{2}f_{0}(β_{0}) ≠ 0. To prove such a property we can write in (3.7)

[m∂_{α}f − q∂_{β}f]_{0}^{(k−1)} = −q∂_{β}^{2}f_{0}(β_{0}) b_{0}^{(k−1)} + all the other terms , (3.8) 3.8

where ∂_{β}^{2}f_{0}(β_{0}) ≠ 0 by hypothesis (see the second relation in (2.12)), and use the tree expansion envisaged in [9] and briefly recalled in Appendix A1, to which we defer for details: here we confine ourselves that the final bound on the radius of convergence gives

ε_{0} = min \{\frac{ω_{0}}{|C|}, \frac{ω_{0}^{2}}{B_{1}} , \frac{|∂_{β}^{2}f_{0}(β_{0})|}{B′_{1}}\} ,  ω_{0} = \frac{1}{q}  (3.9) 3.9

for some positive constants B_{1} and B′_{1}.

Note that (3.6) admits a solution for C ≠ 0 only if the function f_{0}(β) is not identically vanishing, and C is chosen as in (2.11). Therefore theorem 2.3 is proved.

p.3.2 3.2. Proof of theorem 2.6. Suppose now that one has ∂_{β}f_{0}(β) ≡ 0 and γ = Cε^{2}. In such a case equations (3.1) and (3.2) have to be modified into

(iω_{0}ν)^{2} a_{ν}^{(1)} = −m (imνf_{ν} − q∂_{β}f_{ν}) ,  (iω_{0}ν)^{2} b_{ν}^{(1)} = q (imνf_{ν} − q∂_{β}f_{ν}) ,  (3.10) 3.10

for k = 1, and

\begin{aligned}(iω_{0}ν)^{2} a_{ν}^{(k)} & = −m [m∂_{α}f − q∂_{β}f]_{ν}^{(k−1)} − C (iω_{0}ν) a_{ν}^{(k−2)} − C \frac{mp}{q}δ_{ν,0}δ_{k,2}, \\ (iω_{0}ν)^{2} b_{ν}^{(k)} & = q [m∂_{α}f − q∂_{β}f]_{ν}^{(k−1)} − C (iω_{0}ν) b_{ν}^{(k−2)} + pCδ_{ν,0}δ_{k,2}, \end{aligned} (3.11) 3.11

for k ≥ 2.

The equations (3.10) are trivially solvable as ∂_{β}f_{0}(β) ≡ 0 by hypothesis, while the equations (3.10) and (3.11) can be solved for ν ≠ 0 provided that, for ν = 0, one has, for k = 2,

0 = −m [m∂_{α}f − q∂_{β}f]_{0}^{(1)} − C \frac{mp}{q},  0 = q [m∂_{α}f − q∂_{β}f]_{0}^{(1)} + pC.  (3.12) 3.12

while for k ≥ 3 equations (3.5) are satisfied.

The identity [∂_{α}f]_{0}^{(1)} ≡ 0 still holds, so that (3.12) gives

[∂_{β}f]_{0}^{(1)} = \frac{p}{q^{2}} C. (3.13) 3.13

By developing in (3.13) [∂_{β}f]^{(1)} = ∂_{βα}f a^{(1)} + ∂_{ββ}f b^{(1)}, (3.14) 3.14

and by using the expressions of a^{(1)} and b^{(1)} obtained by solving the first order equations (3.10), i.e.

a^{(1)} = −m D^{−2} [m∂_{α}f − q∂_{β}f] ,  b^{(1)} = q D^{−2} [m∂_{α}f − q∂_{β}f] ,  (3.15) 3.15

where D is the operator D = (ω_{0}∂_{ψ}), one finds that

\begin{aligned}[∂_{β}f]^{(2)} & = −m^{2}∂_{αβ}fD^{−2}∂_{α}f − q^{2}∂_{ββ}fD^{−2}∂_{β}f + mq (∂_{βα}fD^{−2}∂_{β}f + ∂_{ββ}fD^{−2}∂_{α}f) \\ & = −∂_{β} (\frac{1}{2}m^{2}∂_{α}fD^{−2}∂_{α}f + \frac{1}{2}q^{2}∂_{β}fD^{−2}∂_{β}f − mq∂_{α}fD^{−2}∂_{β}f) , \end{aligned} (3.16) 3.16

and the r.h.s. is a gradient with respect to β of a function F^{(2)}(α, β): and we check that

\begin{aligned}F_{0}^{(2)}(β) & = − [\frac{1}{2}m^{2}∂_{α}fD^{−2}∂_{α}f + \frac{1}{2}q^{2}∂_{β}fD^{−2}∂_{β}f − mq∂_{α}fD^{−2}∂_{β}f]_{0} \\ & = − ∑_{ν_{1}+ν_{2}=0} (\frac{1}{2}m^{2} (iν_{1}) f_{ν_{1}}(β) (iω_{0}ν_{2})^{−2} (iν_{2}) f_{ν_{2}}(β) \\ & +\frac{1}{2}q^{2}∂_{β}f_{ν_{1}}(β) (iω_{0}ν_{2})^{−2} ∂_{β}f_{ν_{2}}(β) − mq (iν_{1}) f_{ν_{1}}(β) (iω_{0}ν_{2})^{−2} ∂_{β}f_{ν_{2}}(β)) , \end{aligned} (3.17) 3.17

as in (2.15). So that [∂_{β}f]_{0}^{(2)} = ∂_{β}F_{0}^{(2)}(β). (3.18) 3.18

If we choose C as in (2.16) and fix β_{0} as in (2.17), then (3.12) can be solved, while all the equations (3.5) with k ≥ 3 can be solved through a suitable choice of the sequence {b_{0}^{(k)}}_{k∈ℕ}. It remains to check the convergence of the series in ε: we refer to Appendix A1 for details, see also [9]. The bound, derived in Appendix A1, on the radius of convergence reads as

ε_{0} = min \{\sqrt{\frac{ω_{0}}{|C|}}, \frac{ω_{0}^{2}}{B_{2}} , \frac{|∂_{β}^{2}f_{0}(β_{0})|}{B′_{1}}\} (3.19) 3.19

for some positive constants B_{2} and B′_{2}.

p.3.3 3.3. Proof of theorem 2.9. We can combine the results of the two above theorems. The condition on β_{0} follows from a first order analysis, by taking into account that the equations to be used in such a case are (3.10): hence the value of β_{0} is the same as in absence of friction. Then equations (3.12) give a condition on b_{0}^{(1)} which reads as

b_{0}^{(1)} = − (∂_{β}^{2}f_{0}(β_{0}))^{−1} G_{0}^{(2)},
G_{0}^{(2)} = \frac{p}{q^{2}} C + other terms ,

(3.20) 3.20

so that the condition on the radius of convergence becomes

ε_{0} = min \{\frac{q^{2}}{p|C|}, \frac{ω_{0}^{2}}{B_{3}} , \frac{|∂_{β}^{2}f_{0}(β_{0})|}{B′_{3}}\} (3.21) 3.21

for some positive constants B_{3} and B′_{3}.

p.3.4 3.4. Remark. (1) The constants B_{j}, B′_{j} can be explicitly computed, see Appendix A1 for details. (2) Note that with respect to [9] there are no small divisors, and all propagators 1/(ω_{0}ν) can be bounded by 1/ω_{0} (for ν ≠ 0).

4. Simple examples with background friction sec.4

p.4.1 4.1. Pendulum with friction and forcing. Consider the case (2.13) with p = q = 1 (so that n = 1 and m = 2 in (2.4)), and choose ε > 0 and C > 0. Then from (2.14) we obtain

∂_{β}f_{0}(β) = −\frac{1}{2}∂_{β} cos β = \frac{1}{2} sin β, (4.1) 4.1

so that the first of (2.12) gives β_{0} = arcsin 2C, (4.2) 4.2

provided 2C < 1. This means that the periodic solution X(t) has α_{0} = −2β_{0}, hence θ_{0} = −β_{0} = − arcsin 2C. Note that (4.2) corresponds to two solutions β_{1} ∈ (0, π/2) and β_{2} = π − β_{1}.

If c =^{def} cos arcsin 2C a (tedious) computation (see (A2.14)) gives for the Lyapunov multipliers (also known as Floquet multipliers) describing the stability of the orbit the values

λ_{±} = 1 ± \sqrt{2επ^{2}c} + (−Cπ + π^{2}c) ε + O(ε^{3/2}) and
c > 0  for β_{0} = β_{1}, λ_{±} = 1 ± a_{1}\sqrt{ε} + O(ε),  a_{1} ≡ \sqrt{2π^{2}c} > 0,
c < 0  for β_{0} = β_{1}, λ_{±} = 1 ± ia_{2}\sqrt{ε} + O(ε),  a_{2} ≡ \sqrt{−2π^{2}c} > 0

(4.3) 4.5

Therefore we can conclude that, to first order,one of the two periodic solutions is stable and the oher unstable, in agreement with the cited Poincaré–Birkhoff’s theorem.

It would be interesting to study (at least numerically) the basins of attractions of the stable periodic solutions: the numerical analysis in [4] suggests that, for values of the parameters not too large (as it is certainly the case in the perturbative regime), the union of the basins of attractions of all periodic orbits fills the whole phase space. Then a comparative study of the basins should determine which periodic solutions attract most of trajectories.

We have also investigated numerically what happens when the friction decreases in time, slowly with respect to the characteristic periods of the system, e.g. with an exponential decay law γ = Cεe^{−κt}, with κ small: the periodic orbits which are stable in correspondence of the value of the initial friction continue to exist and all trajectories appear to be attracted by such orbits if κ is small enough (with respect to the time needed by the attractor to be reached, i.e. a time of the order of the inverse of the Lyapunov coefficients of the stable orbits thus singled out).

p.4.2 4.2. Periodically forced gyroscope with background friction. Consider the equation (2.1) with

g(θ, t) = ∑_{j∈ℕ} a_{j} cos(2θ − jt), (4.4) 4.6

which can represent the precession of a gyroscope subjected to a periodic torque.

For instance the gyroscope could be moving on an ellipse of eccentricity e and it could be subjected to a gravitational attraction from a mass located in the ellipse focus. For concreteness one can consider the g(θ, t) arising in a spin–orbit planetary model: this is interesting because it gives us the possibility of introducing a well known model (see equation (2) in [10]) which may have some relevance in Astronomy and that we shall study in Section 5 under the presence of tidal friction. The model has been used (without friction) in [5] to study the stability of librations in the cases of a few celestial systems. It will illustrate the mechanism of resonance selection which depends in a delicate way on the relative size of the coefficients a_{j} (which depend on the eccentricity e of the orbit of the planet, see (5.1)). There are 10 such coefficients

The construction of the model is reproduced in Appendices A3 and A4, see (A4.23), and gives for the function g(θ, t) the expression (5.1) in next Section, from which the coefficients can be obtained once the value of the eccentricity is known. The coefficients have been computed, in our numerical tests, via an algebraic manipulator from the data of a few celestial bodies by assigning appropriate values to the eccentricity e.

We shall use data that arise in the following six cases: Moon-Earth, Mercury-Sun, Io-Jupiter, Enceladus-Saturn, Dione-saturn and Rhea-Saturn. We stress, however, that the following has no pretention of being a study in Celestial Mechanics because the background friction does not seem to be a sensible model for the capture into resonance of the systems considered.

For ω = p/q one can write f(α, β) = ∑_{j=−3,j≠0}^{7} a_{j} cos ((2p − jq)α + (2n − jm)β) , so that

f_{0}(β) = a_{j_{0}} cos ((2n − j_{0}m)β) ,  2p − j_{0}q = 0,  (4.5) 4.7

which gives the values listed in table 4.1.

j_{0} p q n m f_{0}(β)
1 1 2 1 3 a_{1} cos β
2 1 1 1 2 a_{2} cos 2β
3 3 2 1 1 a_{3} cos β
4 2 1 1 1 a_{4} cos 2β
5 5 2 2 1 a_{5} cos β
6 3 1 2 1 a_{6} cos 2β
7 7 2 3 1 a_{7} cos β

Table 4.1. Values of p, q, n, m and the corresponding f_{0}(β) for different resonances j_{0} : 2, i.e. p/q = j_{0}/2, with j_{0} = 1, …, 7 in (4.5).

Therefore (2.12) fixes β_{0} so that

β_{0} = \frac{1}{δ} arcsin (−\frac{p}{q^{2}} \frac{C}{a_{j_{0}}} \frac{1}{δ}) ,  δ = \{ \begin{matrix}1, & j_{0} = 1, 3, 5, 7, \\ 2, & j_{0} = 2, 4, 6, \end{matrix}  (4.6) 4.8

provided that one has

|\frac{p}{q^{2}} \frac{C}{a_{j_{0}}} \frac{1}{δ}| < 1. (4.7) 4.9

We define C_{j_{0}} to be the positive value of C for which one has equality in (4.7) and call this quantity the existence threshold for the resonance corresponding to j_{0}.

j_{0}  1  2  3  4  5  6  7
p/q  1/2  1/1  3/2  2/1  5/2  3/1  7/2
M-E q^{2}|a_{j_{0}}|δ/p 0.05488 0.99247 0.12725 0.01272 0.00116 0.00010 0.00001
M-S q^{2}|a_{j_{0}}|δ/p 0.20553 0.89475 0.43787 0.16398 0.05574 0.02024 0.00640
I-G q^{2}|a_{j_{0}}|δ/p 0.00410 0.99996 0.00957 0.00007 0.00000 0.00000 0.00000
E-S q^{2}|a_{j_{0}}|δ/p 0.00450 0.99995 0.01050 0.00009 0.00000 0.00000 0.00000
D-S q^{2}|a_{j_{0}}|δ/p 0.00220 0.99999 0.00513 0.00002 0.00000 0.00000 0.00000
R-S q^{2}|a_{j_{0}}|δ/p 0.00100 1.00000 0.00233 0.00000 0.00000 0.00000 0.00000

Table 4.2. Values of C_{j_{0}} which characterize the appearance of periodic orbits with periods p/q for the cases Moon-Earth (M-E), Mercury-Sun (M-S), Io-Jupiter (I-J), Enceladus-Saturn (E-S), Dione-Saturn (D-S) and Rhea-Saturn (R-S).

Thus we see that for C large, i.e. at least for C > C_{j_{0}} in our (not optimal) estimates, there is no periodic orbit close to an unperturbed one, while for values of C just below C_{2} only the periodic orbit with frequency ω = 1 exists. Another periodic orbit (with frequency 3/2) appears when C falls below the value C_{3}. In the case of Mercury-Sun one has C_{3} ≈ C_{2}/2, while in all other cases one has at best C_{3} ≈ C_{2}/8. This should be in agreement with the fact that Mercury-Sun is the only case in which at the end the system was captured into the resonance 3:2; however for a physical interpretation of the results we defer to next Section.

In all cases, in concrete applications, one should check that the values of ε of interest are less than the value of the radius of convergence of the perturbative series (relative to each orbit). Here we have not fixed the value of ε because the present analysis is an illustration rather than an application: hence we are free to take ε very small, as small as necessary to insure convergence of the series that define the resonant orbits. We shall discuss this matter, i.e. how to fix ε, in the context of Celestial Mechanics applications when the friction is tidal.

p.4.3 4.3. Lyapunov coefficients of the periodic orbits As far as the stability of the periodic solutions are concerned, one has α_{0} = −mβ_{0}/q + O(ε), hence θ_{0} = −β_{0}/q + O(ε), with β_{0} given by (4.6). As in deriving (4.3) a trivial computation (but still more tedious than the previous one) gives, for the cases Moon-Earth (with j_{0} = 2) and for Mercury-Sun (with j_{0} = 3) or for Mercury-Sun with j_{0} = 2, respectively

λ_{±} = 1 ± 8.85229\sqrt{ε cos 2β_{0}} + (−Cπ + 39.18154 cos 2β_{0}) ε + O(ε^{3/2}),
λ_{±} = 1 ± 14.35847\sqrt{ε cos β_{0}} + (−2Cπ + 103.08289 cos β_{0}) ε + O(ε^{3/2}),
λ_{±} = 1 ± 8.40999\sqrt{ε cos 2β_{0}} + (−Cπ + 35.36398 cos 2β_{0}) ε + O(ε^{3/2}),

(4.8) 4.10

where the stable periodic orbits are those with cos δβ_{0} < 0.

5. Spin-orbit interaction with tidal friction sec.5

p.5.1 5.1. Application to the spin-orbit model. In a system satellite-planet (or planet-Sun) there can be several types of friction: the friction between satellite layers of different composition, say one liquid and one solid (core-mantle friction), or the friction due to the tides (tidal friction). There could be also other sources of friction which we do not consider (especially those which could modify the revolution motion of the planet that, as we are implicitly using that it occurs on a fixed orbit). One can expect that such phenomena produce a friction term of the form γθ̇ r(θ, θ̇, t) + γ′, with γ, γ′ ∈ ℝ, and r a 2π–periodic function (in θ and t) analytic and positive, to be added to the forcing term G(θ, t) = ∂_{θ}g(θ, t), where, [5] and appendix A4, A5

g(θ, t) = (\frac{1}{4}e − \frac{1}{32}e^{3} + \frac{5}{768}e^{5}) cos(2θ − t) + (\frac{1}{2} − \frac{5}{4}e^{2} + \frac{13}{32}e^{4}) cos(2θ − 2t)

+ ( − \frac{7}{4}e + \frac{123}{32}e^{3} − \frac{489}{256}e^{5}) cos(2θ − 3t) + (\frac{17}{4}e^{2} − \frac{115}{12}e^{4}) cos(2θ − 4t)

+ ( − \frac{845}{96}e^{3} + \frac{32525}{1536}e^{5}) cos(2θ − 5t) (5.1) 5.1

+ (\frac{533}{32}e^{4}) cos(2θ − 6t) + ( − \frac{228347}{7680}e^{5}) cos(2θ − 7t)

+ ( − \frac{1}{96}e^{3} − \frac{11}{1536}e^{5}) cos(2θ + t) + (\frac{1}{48}e^{4}) cos(2θ + 2t) + ( − \frac{81}{2560}e^{5}) cos(2θ + 3t).

Here the parameter e is the eccentricity of the orbit, and g(θ, t), as given by (5.1), is the θ-dependent part of the power expansion in e of ω^{2}(λ_{T})ω_{T}^{−2} cos^{2}(θ − λ_{T}) (see Appendices A3 and A4 for details and notations). The perturbative parameter ε is related to the asymmetry δ_{2} = (I_{y} − I_{x})/I_{z} of the equatorial moments of inertia, which we suppose to be due to the tidal equatorial bulge of height h, so that δ_{2} = 2h/R, while C is related to the viscosity of the fluid magma forming the planet:

ε =^{def}δ_{2} = \frac{3}{2} \frac{2h}{R},  C =^{def} \frac{1}{3} \frac{Rη}{Mω} .  (5.2) 5.2

The above model can be derived from the theory of rigid motions as explained (for completeness) in Appendix A4.

Here we consider the case in which the friction is fixed as in (2.1), with μ = 1, and fix γ = Cε. The choice reflects that one expects friction to be due to tidal phenomena and to be minimized in a resonance 1:1; see Section 5.2 for further comments. Then the equations of motion become

θ̈ + εG(θ, t) + Cε (θ̇ − 1) = 0, (5.3) 5.3

with G(θ, t) = ∂_{θ}g(θ, t) given by (5.1). See table 5.1 for the values of ε and γ.

Planet γ = Cε  ε  C
Moon 1.00e − 02 3.38e − 05 (2.96e + 02)
Mercury 6.06e − 04 2.03e − 06 (2.99e + 02)
Io  1.26e − 01 7.66e − 03 (1.64e + 01)
Encelado 6.51e + 01 3.49e − 02 (1.86e + 03)
Dione 5.33e + 00 7.98e − 03 (6.68e + 02)
Rea 2.31e + 00 3.38e − 03 (6.85e + 02)

Table 5.1 Values γ, ε, C for the various cases according to the model developed in Appendix A3,A4 and in Section 5.2; they are the values adopted in this paper.

We change variables in a way suitable to study one of the seven resonances (corresponding to the seven terms in (5.1): this step is the same as the corresponding one performed in the case of the gyroscope in Section 4. The resonances will be labeled in the same way by j_{0} = 1, 2 … 7. The coefficients a_{j} are the same; see table 4.1.

By recalling (2.12) we see that one has to fix β_{0} so that

β_{0} = \frac{1}{δ} arcsin (−\frac{p − q}{q^{2}} \frac{C}{a_{j_{0}}} \frac{1}{δ}) ,  provided  |\frac{p − q}{q^{2}} \frac{C}{a_{j_{0}}} \frac{1}{δ}| < 1 (5.4) 5.4

with δ = 1 if j_{0} = 1, 3, 5, 7 and δ = 2 if j_{0} = 2, 4, 6. As a consequence the table 4.3 is replaced by table 5.2.

j_{0}  1  2  3  4  5  6  7
p/q  1/2 1/1 3/2  2/1  5/2  3/1  7/2
M-E \frac{q^{2}|a_{j_{0}} |δ}{p−q} 0.05488 ∞ 0.38176 0.02545 0.00193 0.00015 0.00001
M-S \frac{q^{2}|a_{j_{0}} |δ}{p−q} 0.20551 ∞ 1.30786 0.32796 0.09290 0.03036 0.00896
I-G \frac{q^{2}|a_{j_{0}} |δ}{p−q} 0.00410 ∞ 0.02870 0.00014 0.00000 0.00000 0.00000
E-S \frac{q^{2}|a_{j_{0}} |δ}{p−q} 0.00450 ∞ 0.03150 0.00017 0.00000 0.00000 0.00000
D-S \frac{q^{2}|a_{j_{0}} |δ}{p−q} 0.00220 ∞ 0.01540 0.00004 0.00000 0.00000 0.00000
R-S \frac{q^{2}|a_{j_{0}} |δ}{p−q} 0.00100 ∞ 0.00700 0.00001 0.00000 0.00000 0.00000

Table 5.2. Estimated values C_{j_{0}} of C below which attractive periodic orbits with periods p/q exist, for the cases Moon-Earth (M-E), Mercury-Sun (M-S), Io-Jupiter (I-J), Enceladus-Saturn (E-S), Dione-Saturn (D-S) and Rhea-Saturn (R-S). The entries 0.00000 mean that the first nonzero digit is beyond the ones written.

Evaluating the thresholds for the various cases we can estimate the various constants and we get that the radius of convergence of the series defining the periodic orbit starts being positive for C < C_{j_{0}} and rapidly tends to a limit as C → 0. If C < C_{j_{0}}/2, say, the radius of convergence has a size ε_{j_{0}} that should be compared with the actual value of ε in the case of each planet. The results of our estimates is reported in table 5.3.

j_{0}  1  3  4  5  6  7
p/q  1/2  3/2  2/1  5/2  3/1  7/2
Moon 1.23e + 00 1.28e + 01 3.05e + 01 1.94e − 01 3.63e − 01 1.99e − 03
Mercury 8.52e + 01 7.18e + 02 7.28e + 03 1.73e + 02 1.35e + 03 2.78e + 01
Io  4.02e − 04 4.22e − 03 7.50e − 04 3.57e − 07 4.94e − 08 2.03e − 11
Encelado 9.68e − 05 1.02e − 03 1.98e − 04 1.04e − 07 1.57e − 08 7.08e − 12
Dione 2.07e − 04 2.18e − 03 2.08e − 04 5.30e − 08 3.94e − 09 8.66e − 13
Rea 2.23e − 04 2.34e − 03 1.01e − 04 1.18e − 08 3.97e − 10 3.97e − 14

Table 5.3. Values ε_{j_{0}}/ε (for C = C_{j_{0}}/2). If the value is > 1 the theory can be applied to the resonance corresponding to j_{0}. The values of ε are from Table 5.1.

A rough estimate of ε_{j_{0}} easily follows from the proof in Appendix A1 and it gives a growth of the convergence radius as C → 0 proportional to

|ε| < const\sqrt{1 − \frac{|p − q|}{q^{2}} \frac{C}{|a_{j_{0}}|}},  p ≠ q.  (5.5) 5.5

Table 5.3 is based on a more detailed estimate following the proof details in Appendix A1.

The interest of table 5.3 is that it shows that aside from the cases of Io, Enceladus, Dione, Rhea the resonance 3:2 falls in the domain of convergence of our theory. In the case of the latter four heavenly bodies the interesting resonance 3:2 (as well as all others) is outside the convergence radius as estimated in Appendix A1, see (A1.15). The estimates, however, might be improvable beyond our crude results derived in Appendix A1, although it seems unlikely that the cases of the four satellites could be covered. However attempting an improvement would not be really useful because the time scale for the attraction by the 1:1 resonance turns out to be too short for all heavenly bodies except for Mercury, see below.

The main difference with respect to the case studied in Section 4 is that the periodic orbit with frequency 1 exists for all values of C; besides that we can immediately realize that no real qualitative change is produced with respect to the results obtained in the previous section.

p.5.7 5.2. A physical friction model. A detailed discussion would require assuming a more realistic model for the satellite structure and relying on some theory of evolution of the planetary system; see for instance [19] for an introduction on the subject.

Suppose that friction is essentially due to tidal effects on an entirely fluid “fast” rotating planet (i.e. with rotation speed ω at least a few times larger than the revolution speed) evolving toward a solid body with neglegible tidal effects. A model for tidal friction can be derived by assuming a tide on a fluid planet which creates a bulge lagging behind the radius vector from the center of attraction S to the planet by an angle τ. In this way the planet looses its equatorial symmetry and its smallest axis of inertia points in a direction at an angle τ with respect to the line joining the planet center to the attraction center S.

This means that in a frame fixed with the planet the axes of inertia rotate at speed λ̇_{T} and the equations of motion of the body are affected, to first order, by the presence of a torque which in the comoving frame is 2δ_{2} sin τ(θ̇ − ω_{T}), where δ_{2} = (I_{x} − I_{y})/I_{z} = (a^{2} − b^{2})/(a^{2} + b^{2}) = 2h/R, if h is the tide height and R the planet radius (because we suppose a = R − h, b = R + h), see Appendix A4. Furthermore the phase shift τ will be taken simply τ = γ_{0} (θ̇ − 1) i.e. proportional to |θ̇−ω_{T} | of the order of ω_{T} = 1 (in our dimensionless units) or smaller. Also the rotation velocity ω_{D} should enter in the dimensional analysis but we assume it to be of the order of ω_{T} because we are studying simple resonances.

The constant γ_{0} has to be a dimensionless constant formed with the following physical quantities: the viscosity η of the magmatic fluid constituting the planet, the planet radius R, the tide height h and the angular velocity ω_{T}. Simplest is to take γ_{0} = ηh/Mω_{T}. Defining ε, C as in (5.2) we obtain the equations (5.3) and the coefficients can be computed from the tidal excursion h.

The tidal excursion (as long as the planet is fluid, see [8], problem 1.4.10) is of the order of R(R/ρ)^{3}(M_{0}/M), if R is the radius of the planet, ρ is the distance between planet and central body, and M_{0} is the mass of the latter. In the early history of the solar system one can assume that the viscosity is very large: larger than that of magma, which can exceed that of water (which equals 10^{−2} C.G.S units) by a factor up to ∼ 10^{12} and smaller than that of the Earth mantle (of about ∼ 10^{22} times the water viscosity); see for instance [16] and [18].

By assuming for η a value 10^{−15} C.G.S. units and inserting the values of R, M, ω, ρ, M_{0} of the various systems considered so far (see for instance [2]), we find that the order of γ_{0} at the beginning of the life of the celestial body (in units of the inverse of the revolution period) can be quite different depending on the body and it can be essentially negligible in some cases (for instance it is around 6 10^{−4} for Mercury, and negligible for the Moon and the other satellites). This means that the characteristic time for the approach to the 1:1 resonance is quite fast for Mercury (a halving time for the spin of about 10^{3} years) and much faster for the Moon, see table 5.1. Halving the value of η does not affect the values of C_{j_{0}} but it doubles the time of the spin down: hence the actual value of η is quite important.

At the same time the value of γ = Cε has to decrease, in the case of Mercury, by a factor of about ∼ 10^{2} before the value of C reaches the treshold of stability of the first non-trivial resonance (i.e. 3:2) and of the order of 4 times more to reach the stability treshold for the (next) resonance 1:2 (compare column C in table 5.1 and columns 3 and 4 in table 5.2).

In the case of the Moon the value of γ has to decrease by a factor ∼ 10^{3} for the stability of the 3:2 resonance and 15 times more for the (next) 2:1 (again compare column C in table 5.1 and columns 3 and 4 in table 5.2).

The cooling, i.e. the moment when the viscosity can be regarded as much smaller because the body solidifies, should take quite likely the same time in both cases. Hence there is the possibility envisaged in the introduction that the planet will spin down while cooling and that its solidification, when the friction becomes comparatively rapidly neglegible, occurs when the spin is still greater than 3:2 so that the planet can be captured by the 3:2 resonance which has become stable. The cases of Mercury and of the Moon seem best suited for the appearance of the stable 3:2 resonance because it has the lowest treshold (provided the solidification time is of the order of 10^{4} revolutions of Mercury, i.e. ∼ ·10^{3} years, or ∼ 10^{2} revolutions of the Moon), see table 5.1): this seems to make Mercury the only case in which 3:2 can become a stable resonance (as the Moon would collapse on the 1:1 resonance in time much shorter because γ is quite large, a time of order of 10^{2} lunar years, i.e. ∼ 10^{2} months or ∼ 8 years).

The choice of the value of η = 10^{15} C.G.S. units affects the results but qualitatively it has the only effect of changing the time scales by a common factor. Unfortunately without a detailed model for the formation of rocky planets and the evolution of the dynamic viscosity the above seems all one can say at the moment.

Acknowledgments. We are indebted to Alessandra Celletti for useful discussions and explanations about the spin-orbit model.

Appendix A1. Tree expansion and technical details of the proofs app.A1

p.A1.1 A1.1. Trees. We prefer to follow the tree formalism introduced in [9]: other approaches based on the implicit functions theorem could also be applied, if preferred. A tree θ is defined as a partially ordered set of points, connected by lines. The lines are oriented toward the root, which is the leftmost point; the line entering the root is called the root line. If a line ℓ connects two points v_{1} and v_{2} and is oriented from v_{2} to v_{1} we say that v_{2} ≺ v_{1} and we shall write v′_{1} = v_{2} and ℓ = ℓ_{v_{2}}; we shall say also that ℓ exits from v_{2} and enters v_{1}.

root  ℓ_{0} v_{0}  v_{1} v_{2} v_{3}  v_{4}  v_{5}  v_{6} v_{7} v_{8} v_{9} v_{10} v_{11} v_{12}

Fig. A1.1. A tree θ with 9 nodes and 4 leaves: the latter ones are graphically represented as bullets. The length of the lines should be the same but it is drawn of arbitrary size.

Besides the root, there will be two kinds of points: the nodes and the leaves. The leaves can only be endpoints, i.e. they have no lines entering them, but an endpoint can be either a node or a leaf. The lines exiting from the leaves play a very different rôle with respect to the lines exiting from the nodes, as we shall see below. We shall denote by v_{0} the last (i.e. leftmost) node of the tree, and by ℓ_{0} the root line; for future convenience we shall write v′_{0} = r but r will not be considered a node. See figure A1.1 for an example of tree.

We shall denote by V (θ) the set of nodes, by L(θ) the set of leaves and by Λ(θ) the set of lines. Fixed any line ℓ_{v} ∈ Λ(θ), we shall say that the subset of θ containing ℓ_{v} as well as all nodes

w ⪯ v and all lines connecting them is a subtree of θ with root v′: of course a subtree is a tree.

Given a tree, with each node v we associate a mode label ν_{v} ∈ ℤ, and to each leaf v a leaf label κ_{v} ∈ ℕ. The quantity

k = |V (θ)| + ∑_{v∈L(θ)} κ_{v} (A1.1) A1.1

is called the order of the tree θ. With any line ℓ exiting from a node v we associate a label γ_{ℓ} assuming the symbolic values α, β and a momentum label ν_{ℓ} ∈ ℤ, which is defined as

ν_{ℓ} ≡ ν_{ℓ_{v}} = ∑_{{w∈V (θ) \atop w⪯v}} ν_{w}, (A1.2) A1.2

while with any line ℓ exiting from a leaf v we associate only the label γ_{ℓ} = β.

We can associate with each node also some labels depending on the entering lines and on the exiting one: the branching labels r_{v} and s_{v}, denoting how many lines ℓ having the label γ_{ℓ} = α

and, respectively, γ_{ℓ} = β enter v, the label η_{v}, defined as

η_{v} = \{ \begin{matrix}1, & if γ_{ℓ_{v}} = β , \\ 0, & if γ_{ℓ_{v}} = α , \end{matrix} (A1.3) A1.3

and the label j_{v} assuming the values 1, 2, 3. The label j_{v} indicates which of the three terms in (3.1) is selected.

Then with each node v we associate a node factor

\begin{aligned}F_{v} = & (−δ_{η_{v},α}m + δ_{η_{v},β}q) (δ_{j_{v},1}miν_{v} − δ_{j_{v},2}q∂_{β}) \frac{1}{r_{v}!} \frac{1}{s_{v}!}(iν_{v})^{r_{v}}∂_{β}^{s_{v}}f_{ν_{v}}(β_{0}) \\ & − Cδ_{j_{v},3} [δ_{η_{v},α} (δ_{r_{v},1}δ_{s_{v},0} − \frac{mp}{q}δ_{r_{v},0}δ_{s_{v},0}) + δ_{η_{v},β} (δ_{r_{v},0}δ_{s_{v},1} + pδ_{r_{v},0}δ_{s_{v},0})] , \end{aligned} (A1.4) A1.4

which is a tensor of rank r_{v} + s_{v} + 1, while with each leaf v we associate a leaf factor (to be defined recursively, see below)

L_{v} = b_{0}^{(κ_{v})}, (A1.5) A1.5

which is a tensor of rank 1 (i.e. a constant); to each line ℓ ≡ ℓ_{v} exiting from a node v we associate a propagator

G_{ℓ} ≡ (δ_{j_{v},1} + δ_{j_{v},2}) \frac{1}{(iω_{0}ν_{ℓ})^{2}} + δ_{j_{v},3} \frac{1}{iω_{0}ν_{ℓ}} , (A1.6) A1.6

while no divisor is associated with the lines exiting from the leaves. For consistence we can define

G_{ℓ} ≡ 1, (A1.7) A1.7

for lines exiting from leaves, so that a propagator G_{ℓ} is in fact associated with each line.

Call Θ_{k,ν,γ} the set of all trees of order k with ν_{ℓ_{0}} = ν and γ_{ℓ_{0}} = γ, if ℓ_{0} is the root line, and define the application Val: Θ_{k,ν,γ} → ℝ, as

Val(θ) = ( ∏_{v∈V (θ)} F_{v})( ∏_{v∈L(θ)} L_{v})( ∏_{ℓ∈Λ(θ)} G_{ℓ}), (A1.8) A1.8

which is called the value of the tree θ.

We can define also the reduced value of the tree θ as

Val′(θ) = F̃_{v_{0}}( ∏_{v∈V (θ)\v_{0}} F_{v})( ∏_{v∈L(θ)} L_{v})( ∏_{ℓ∈Λ(θ)\ℓ_{0}} G_{ℓ}), (A1.9) A1.9

where, as usual, ℓ_{0} denotes the root line, and F̃_{v_{0}} is defined as

F̃_{v_{0}} = δ_{j_{v_{0}} ,1} (miν_{v_{0}} − q∂_{β}) \frac{1}{r_{v_{0}}!} \frac{1}{s_{v_{0}}!}(iν_{v_{0}})^{r_{v_{0}}} ∂_{β}^{s_{v_{0}}} f_{ν_{v_{0}}} (β_{0}); (A1.10) A1.10

then, by setting Θ_{k,0,γ}^{∗} ≡ Θ_{k,0,γ} \ θ_{0}, where θ_{0} is the tree given in figure A1.2, one can define the leaf factor (A1.5) as

b_{0}^{(k)} = − [∂_{β}^{2}f_{0}(β_{0})]^{−1} ∑_{θ∈Θ_{k+1,0,γ}^{∗}} Val′(θ), (A1.11) A1.11

where the quantity ∂_{β}^{2}f_{0}(β_{0}) is nonvanishing by the second condition in (2.12). Note that in (A1.11) the value of the label γ is irrelevant as Val′(θ) does not depend on ℓ_{0}.

θ_{0} = r  v_{0}  v ν_{v_{0}} = 0 κ_{v} = k

Fig. A1.2. The tree θ_{0} with one node v_{0} and one leaf v. One has r_{v_{0}} = 0 and s_{v_{0}} = 1. The order k of the tree is given by k = 1 + κ_{v}, and the momentum is vanishing, so that ν_{v_{0}} = 0.

Then it is straightforward to prove by induction that the following tree expansion holds for the coefficients in (2.10) of the periodic solution (2.9):

a_{ν}^{(k)} = ∑_{θ∈Θ_{k,ν,α}} Val(θ), b_{ν}^{(k)} = ∑_{θ∈Θ_{k,ν,β}} Val(θ), (A1.12) A1.12

for k ∈ ℕ and ν ∈ ℤ \ {0}, while b_{0}^{(k)} is given by (A1.11) and a_{0}^{(k)} ≡ 0 has to be fixed in such a way that the identity a(α_{0}, β_{0}; ε) = −mb(α_{0}, β_{0}; ε)/q is satisfied to all orders.

p.A1.2 A1.2. Bounds. It is easy to realize that, for b_{0}^{(k)} fixed as in (A1.11) no line ℓ ∈ Λ(θ) can have ν_{ℓ} = 0, except for the lines exiting from the leaves; but for them the propagator has no divisor (see (A1.7)) and it is replaced by ∂_{β}^{2}f_{0}(β_{0}).

By using (A1.11) we can decompose iteratively the leaves into trees, so that at the end only trees without leaves appear. For such trees (A1.1) has to be replaced with k = |V (θ)|; we can write k = k_{1} + k_{2}, where k_{1} is the number of nodes v with j_{v} ∈ {1, 2} and k_{2} is the number of nodes v with j_{v} = 3.

Then, by using that

\begin{matrix}|\frac{1}{s_{v}!}∂^{s_{v}+1}f_{ν_{v}}(β_{0})| ≤ F D^{s_{v}}e^{−κ|ν_{v}|}, \\ |\frac{1}{r_{v}!}(iν_{v})^{r_{v}+1}| ≤ (r_{v} + 1) (\frac{8}{κ})^{r_{v}+1} e^{κ|ν_{v}|/8}, \\ ∑_{v∈V (θ)} (r_{v} + s_{v}) = k − 1, \end{matrix} (A1.13) A1.13

for some positive constants F, D, κ, one obtains easily

| ∏_{v∈V (θ)} F_{v}| ≤ e^{−κ|ν|/2}B_{1}^{k_{1}}|C|^{k_{2}},  | ∏_{ℓ∈Λ(θ)} G_{ℓ}| ≤ (\frac{1}{ω_{0}^{2}} )^{k_{1}} (\frac{1}{ω_{0}} )^{k_{2}} ,  (A1.14) A1.14

for a suitable constant B_{1}. If we take into account also the leaves we have no propagator to associate to the root line, but we have to take into account the factor before the sum in (A1.11), so that (3.9) follows.

A better bound can be found in the case in which f_{0} contains only one harmonic, as in the case (4.4) and as in Section 5; in such cases given the resonance j_{0} of type p : q (i.e. ω = \frac{p}{q}) we consider the corresponding m, n, see (2.4), and set ν = 2p − jq, μ = (2n − jm) for j = 1, …, 7 then convergence occurs for

ε < ( max{a, c} · b · d)^{−1}, with

a =^{def} max{m, q} · max{m|ν|, q|μ|} · max_{j=1,...,7}{|ν||a_{j}|, |μ||a_{j}|},

b = max{q^{2}, 1/((2n − j_{0}m)^{2}|a_{j_{0}}| cos(β_{0}))},

c = C max \{ max \{1, \frac{m|p − q|}{q}\}, max\{1, p\}\},

d = 3 · 2 · 2 · 4 = 48, (A1.15) A1.15 where the first term in a arises from the quantity in the first parenthesis in (A1.4), the second and third arise from the next two factors; the c arises from the term proportional to C in (A1.4); the factor b is due to the propagators attributed to the lines (using ω_{0} = \frac{1}{q} for the first term in the maximum, coming from the lines with label γ = α, β, while the second term comes from the leaves); the factors in d are due to the choices of the values of j_{v} (3), choices of the labels γ (i.e. 2), choices of the nodes representing leaves (2) and choices of the tree structures (4), see [9]. The (A1.15) leads to the table 5.3 (for j_{0} ≠ 2) which is evaluated at C equal to \frac{1}{2}C_{j_{0}} if C_{j_{0}} is the maximal value estimated below which the resonance in question is stable (see Table 5.2 and (5.5)).

p.A1.3 A1.3. About the proof of the other theorems. The technical parts of the proofs of theorems 2.6 and 2.9 (and of theorem in Appendix A5) can be dealt with in the same way, up to slight changes which we leave to the reader.

Appendix A2. Stability of periodic solutions app.A2

p.A2.1 A2.1. Linearization around the periodic solutions. The periodic orbits for our system appear in pairs of stable and unstable orbits: this is a consequence of Poincaré-Birkhoff’s theorem [1]. To detect the stable solutions from the unstable ones a first order computation is enough.

Let us consider the case dealt with through theorem 2.3. The corresponding equations of motion, in the original variables, are

\{ \begin{matrix}θ̇ = ω + Θ, \\ ṫ = 1, \end{matrix}  \{ \begin{matrix}Θ̇ = −ε∂_{θ}g − εC (ω + Θ) , \\ Ṫ = −ε∂_{t}g. \end{matrix}  (A2.1) A2.1

It is clear that for a study of the stability it is enough to consider the first and third equation, i.e. the system

θ̇ = ω + Θ,  Θ̇ = −ε∂_{θ}g − εC (ω + Θ) .  (A2.2) A2.2

The variable T indeed does not enter in the other equations. On the other hand we can study the stability of a periodic orbit by looking at the corresponding periodic point for the Poincaré map defined by t ∈ 2πℤ. This allows us to disregard variations in the inial value of t.

The linearization of (A2.2) around the periodic solution X(t) = (θ(t), Θ(t)) given by the theorem 2.3, leads to

Ξ̇ = L(t) Ξ,  Ξ = (δθ, δΘ) ,  (A2.3) A2.3

where

L(t) = ( \begin{matrix}0 & 1 \\ −ε∂_{θθ}g(θ(t), t) & −εC \end{matrix} ) ,  L(t) = ∑_{k=0}^{∞} ε^{k}L^{(k)}(t),  (A2.4) A2.4

with

L^{(0)}(t) ≡ L^{(0)} = ( \begin{matrix}0 & 1 \\ 0 & 0 \end{matrix} ) ,  L^{(1)}(t) = ( \begin{matrix}0 & 0 \\ −∂_{θθ}g(θ_{0} + ωt, t) & −C \end{matrix} ) ,  (A2.5) A2.5

and so on.

p.A2.2 A2.2. Wronskian matrix. Let us denote by W (t) = ∑_{k=0}^{∞} ε^{k}W^{(k)}(t) the Wronskian matrix, that is the matrix whose columns are two independent solutions of the linearized system (A2.3) (so that Ẇ (t) = L(t)W (t)). Then one has

W^{(0)}(t) = exp [tL^{(0)}] = ( \begin{matrix}1 & t \\ 0 & 1 \end{matrix} ) , (A2.6) A2.6

while W^{(1)}(t) is obtained by solving the system Ẇ^{(1)} = L^{(0)}(t)W^{(1)}(t) + L^{(1)}(t)W^{(0)}(t), i.e.

W^{(1)}(t) = W^{(0)}(t) [W^{(1)}(0) + ∫_{0}^{t} dτ (W^{(0)}(τ))^{−1} L^{(1)}(τ)W^{(0)}(τ)] , (A2.7) A2.7

where one has to take W^{(1)}(0) = 0 in order to have W (0) = 11.

A trivial computation shows that in (A2.7) one has

W^{(0)}(τ)^{−1}L^{(1)}(τ)W^{(0)}(τ) = ( \begin{matrix}1 & −τ \\ 0 & 1 \end{matrix} )( \begin{matrix}0 & 0 \\ x(τ) & −C \end{matrix} )( \begin{matrix}1 & τ \\ 0 & 1 \end{matrix} ) = ( \begin{matrix}−τx(τ) & −τ^{2}x(τ) + Cτ \\ x(τ) & τx(τ) − C \end{matrix} ),

(A2.8) A2.8 where we have set x(t) = −∂_{θθ}g(θ_{0}+ωt, t). Denote by x_{k}(t) the kth primitive of x(t) with x_{k}(0) = 0 (so that ẋ_{k}(t) = x_{k−1}(t), with x_{0}(t) = x(t)). Then, by using that

\begin{matrix}∫_{0}^{T} dτ x(τ) = x_{1}(T),  ∫_{0}^{T} dτ τx(τ) = T x_{1}(T) − x_{2}(T), \\ ∫_{0}^{T} dτ τ^{2}x(τ) = T^{2}x_{1}(T) − 2T x_{2}(T) + 2x_{3}(T), \end{matrix}  (A2.9) A2.9

we obtain

\begin{matrix}∫_{0}^{T} dτ (W^{(0)}(τ))^{−1} L^{(1)}(τ)W^{(0)}(τ) = \\ = ( \begin{matrix}−T x_{1}(T) + x_{2}(T) & −T^{2}x_{1}(T) + 2T x_{2}(T) − 2x_{3}(T) + CT^{2}/2 \\ x_{1}(T) & T x_{1}(T) − x_{2}(T) − CT \end{matrix} ) , \end{matrix} (A2.10) A2.10

so that

W^{(1)}(T) = ( \begin{matrix}x_{2}(T) & T x_{2}(T) − 2x_{3}(T) − CT^{2}/2 \\ x_{1}(T) & T x_{1}(T) − x_{2}(T) − CT \end{matrix} ) . (A2.11) A2.11

p.A2.3 A2.3. Floquet (or Lyapunov) multipliers. The Floquet multipliers around the periodic solution are the eigenvalues of the Wronskian matrix, computed at time T. Hence, for ε = 0, W^{(0)}(T) =

( \begin{matrix}1 & T \\ 0 & 1 \end{matrix} ), and the corresponding Floquet multipliers are equal to 1. To first order one has

\begin{aligned}W (T) & = W^{(0)}(T) + εW^{(1)}(T) + O(ε^{2}) \\ & = ( \begin{matrix}1 + εx_{2}(T) & T + εT x_{2}(T) − 2εx_{3}(T) − εCT^{2}/2 \\ εx_{1}(T) & 1 + εT x_{1}(T) − εx_{2}(T) − εCT \end{matrix} ) + O(ε^{2}). \end{aligned} (A2.12) A2.13

so that, by neglecting the terms O(ε^{2}) in the Wronskian matrix, the corresponding multipliers are λ_{1} = λ_{+}, λ_{2} = λ_{−}, where λ_{±} are the roots of the equation λ^{2} − 2bλ + c = 0 with

\begin{matrix}2b =^{def} [2 + εT (x_{2}(T) − C)] \\ c =^{def} [1 − εCT + εx_{2}(T)ε^{2} (2x_{1}(T)x_{3}(T) − x_{2}^{2}(T) − CT x_{2}(T) + CT^{2}x_{1}(T)/2)] . \end{matrix} (A2.13) A2.15

Taking the divergence of (A2.2) it is easy to see that the phase space contraction in one period is given by εC. This implies that we must have λ_{−}λ_{+} = 1 − Cε + O(ε^{2}), as it can be checked. Therefore the two Floquet multipliers λ_{±} are given by

λ_{±} = 1 + \frac{εT}{2} (x_{1}(T) − C) ± (A2.14) A2.17

± \sqrt{εT x_{1}(T) + \frac{ε^{2}}{4} (T^{2}C^{2} + T^{2}x_{1}^{2}(T) + 4x_{2}^{2}(T) + 4CT x_{2}(T) − 8x_{1}(T)x_{3}(T) − 4T^{2}x_{1}C)}

= 1 ± \sqrt{εT x_{1}(T)} + \frac{εT}{2} (x_{1}(T) − C) + O(ε^{3/2}).

Appendix A3. Andoyer-Deprit angles trigonometry and rigid bodies app.A3

p.A4.3 A3.1. The spin-orbit model. The potential energy of a body B of ellipsoidal shape with polar radius R_{p} and equatorial radius R_{0} and density ν, attracted by a center of mass M_{0} at distance a from its center is, up to irrelevant (time dependent) additive quantities,

−ν ∫_{B} ρ^{2}dρ sin θ dθ dφ \frac{κM_{0}}{a} (\frac{ρ}{a})^{2}P_{2}(sin α sin θ cos φ + cos α cos θ), (A3.1) A3.1

where κ is the universal gravitation constant, P_{2} is the second Legendre polynomial, P_{2}(z) = (3z^{2} − 1)/2, ρ, θ, φ are the polar coordinates of the point x in the body B of density ν and α is the angle between the polar axis of the ellipsoid and the radius from the center of the ellipsoid to the center of attraction; one has a = a_{0}(1 − e^{2})/(1 − e cos λ_{T}), if e is the eccentricity of the body orbit, λ_{T} is its true anomaly on the orbital plane measured from the apocenter and a_{0} is the orbital semiaxis. Therefore we get that aside from an additive periodic term (depending only on a) the potential energy of the body is, using cilindrical coordinates r, z, φ instead of the polar ρ, θ, φ,

− (2πν) \frac{3}{2} \frac{κM_{0}}{a^{3}} ∫_{−R_{p}}^{R_{p}} dz ∫_{0}^{R_{0}\sqrt{1−(z/R_{p})^{2}}} r dr ρ^{2}(\frac{1}{2} sin^{2} α sin^{2} θ + cos^{2} α cos^{2} θ)

= −\frac{3}{2} \frac{κM_{0}}{a^{3}} [ (2πν) R_{0}^{2}R_{p} 2 cos^{2} α ∫_{0}^{1} dz ∫_{0}^{√\overline{1−z^{2}}} r dr \frac{1}{2}(2z^{2} − r^{2})].

(A3.2) A3.2

Since the term in square brackets equals the difference between the moments of inertia I and J with respect to the x–axis (or y–axis) and with respect to the z axis, respectively, we get

\frac{3}{2}J \frac{κM_{0}}{a_{0}^{3}} δ (\frac{1 − e cos λ_{T}}{1 − e^{2}} )^{3} cos^{2} α,  δ = \frac{J − I}{J} .  (A3.3) A3.3

More generally in the case of an ellipsoid with axes R_{z} < R_{x} < R_{x} (note that in such a case the moments of inertia satisfy I_{z} > I_{y} > I_{x}), if β is the angle between the comoving axis i_{1} of the ellipsoid and the radius from the center of the ellipsoid to the center of attraction, we get, instead of (A3.3),

\frac{3}{2}Jω^{2}(λ_{T})(δ_{1} cos^{2} α + δ_{2} sin^{2} α cos^{2} β). (A3.4) A3.4

where

ω(λ_{T}) = (\frac{1 − e cos λ_{T}}{1 − e^{2}} )^{\frac{3}{2}} ω_{T} ,  ω_{T}^{2} =^{def} \frac{κM_{0}}{a_{0}^{3}} ,  (A3.5) A3.5

with J ≡ I_{z}, δ_{1} =^{def}(I_{z} − I_{x})/J, and δ_{2} =^{def}(I_{y} − I_{x})/J.

Define: m the line of intersection between the plane orthogonal to the angular momentum A and the orbit plane, \overline{n} the line of intersection between the plane orthogonal to the symmetry axis and the orbit plane, n the line of intersection between the planes orthogonal to the symmetry axis and to the angular momentum axis; and call \overline{i} the x axis on the orbit plane (arbitrarily prefixed and from which the true anomaly λ_{T} is measured). Let

φ= angle between m and n, ψ= angle between n and i_{1}, γ= angle between m and \overline{i},

then the pairs (K, γ), (A, φ) and (L, ψ) are canonically conjugated variables, cf. p. 318 in [7], and the Hamiltonian for the system becomes

\frac{A^{2}}{2I_{x}} − δ′ \frac{L^{2}}{2J} − \frac{1}{2}δ_{2} \frac{I_{z}}{I_{x}I_{y}} (A^{2} − L^{2}) cos^{2} ψ + \frac{3}{2}Jω^{2}(λ_{T}) (δ_{1} cos^{2} α + δ_{2} sin^{2} α cos^{2} β), (A3.6) A3.6

where we recall the definitions J = I_{z}, δ_{1} = (I_{z} − I_{x})/I_{z}, and δ_{2} = (I_{y} − I_{x})/I_{z}, and define δ′ = (J − I_{x})/I_{x}.

To proceed one needs to express λ_{T} in terms of the mean anomaly λ = ω_{T} t = t in our units and cos α, sin α cos β in terms of the canonical variables. This can be based upon the identities

cos α = − sin \overline{θ} sin(λ_{T} − \overline{φ}),
sin α cos β = cos \overline{ψ} cos(λ_{T} − \overline{φ}) + cos \overline{θ} sin \overline{ψ} sin(λ_{T} − \overline{φ}),
sin α sin β = − sin \overline{ψ} cos(λ_{T} − \overline{φ}) + cos \overline{θ} cos \overline{ψ} sin(λ_{T} − \overline{φ}),

(A3.7) A3.7

expressing the polar angles of the line joining the center of attraction S to the center of the planet, i.e. its the declination α and the longitude β, in the comoving frame, if the latter is determined by its Euler angles (\overline{φ}, \overline{ψ}, \overline{θ}) with respect to the fixed frame. The algebraic work is discussed below.

p.A4.1 A3.2. Spherical trigonometry and Andoyer-Deprit angles. We refer here to figures A3.1, A3.2 and A3.3 below (see also [7], p. 321÷323) and to the well known spherical triangles trigonometrical identities (see figure A3.1):

\frac{sin A}{sin α} =\frac{sin B}{sin β} = \frac{sin C}{sin γ} ,  “sine rule”
cos A = cos B cos C + sin B sin C cos α,  “cosine rule”
sin C cos β = cos B sin A − sin B cos A cos γ,  “analogue rule”
cos A cos γ = sin A cot B − sin γ cot β,  “four parts rule”

(A3.8) A3.8

α γ  β A  B  C

Fig. A3.1. Spherical triangle with sides A, B, C and angles α, β, γ.

The above rules become well known trigonometric properties of a triangle with sides a, b, c if the spherical triangle is imagined drawn on a sphere of radius R so that the A, B, C are respectively a/R, b/R, c/R. The third relation is analogous to the statement that a side is the sum of the orthogonal projections on it of the other two sides and the last rule follows by combining the sine rule and the analogue rule. Other other identities can be obtained by “duality” by remarking the if (α, β, γ; A, B, C) is a spherical triangle also (A, B, π − C; α, β, π − γ) is a spherical triangle.

The Andoyer-Deprit angles are represented in figure A3.2, where (\overline{i}, \overline{j}, \overline{k}) is the “fixed reference system”, with \overline{k} axis orthogonal to the orbital plane, (i, j, k) is the reference system with k axis parallel to the angular momentum A and i axis on the line between the orbit plane and the plane orthogonal to the angular momentum, and (i_{1}, i_{2}, i_{3}) is the reference frame with i_{3} axis parallel to the symmetry axis of the body and i_{1} axis fixed on the equatorial plane of the body.

\overline{i}  \overline{φ} γ i ≡ m  n φ ψ  \overline{n} i_{1} \overline{ψ} i_{2} \overline{k} i_{3} A  \overline{θ} ζ θ  O  \overline{j}

Fig. A3.2. Andoyer-Deprit angles: here \overline{n} is the node line (\overline{i}, \overline{j}) ∩ (i_{1}, i_{2}), n is the node line (i_{1}, i_{2}) ∩ (i, j) and m ≡ i is the node (i, j) ∩ (\overline{i}, \overline{j}). j axis not drawn. The axis i_{1} has to be thought below the (\overline{i}, \overline{j}) plane (for a drawing consistent with the notations used here.

The above trigonometrical identities imply, see figure A3.3, the following relations

cos ζ =\frac{K}{A} ,  cos θ = \frac{L}{A},
cot(\overline{φ} − γ) =(cos φ cos ζ − sin ζ cot θ)/ sin φ,
cot(\overline{ψ} − ψ) =(cos φ cos θ − sin θ cot ζ)/ sin φ,
sin \overline{θ} = sin θ \frac{sin φ}{sin(\overline{φ} − γ)} = − sin ζ \frac{sin φ}{sin(\overline{ψ} − ψ)},

(A3.9) A3.9

immediately from the definitions, see [7], p. 323.

\overline{θ} ζ  θ φ  \overline{φ} − γ ψ − \overline{ψ}

Fig. A3.3. The spherical triangle associated with the Andoyer-Deprit angles.

Therefore writing λ − \overline{φ} ≡ (λ − γ) − (\overline{φ} − γ) and \overline{ψ} = ψ + (\overline{ψ} − ψ) the relations (A3.9) lead to expressing cos α and sin α cos β, sin α sin β in terms of the Andoyer-Deprit angles.

This is interesting because the pairs (A, φ), (L, ψ), (K, γ) are canonically conjugated pairs of action–angle coordinates. However it is clear that several square roots will appear involving quantities like 1 − K^{2}/A^{2} or 1 − L^{2}/A^{2}, which in the interesting cases will be close to 0, and their reciprocals. Therefore the resulting equations of motion will be difficult to discuss unless we can show that such a priori singular expressions do not really appear. We call this the “rationalization” of the Hamiltonian and to show it we study a few identities relating the angles introduced above.

Appendix A4. Rationalization app.A4

p.A4.x A4.1. Rationalization: L ≃ K ≃ A symmetric case. Let c_{x} ≡ cos x, s_{x} ≡ sin x, J = I_{z}, I = I_{x} = I_{y}. Then cos α can be written as

(s_{λ_{T} −γ}(c_{φ}c_{ζ}s_{θ} − s_{ζ}c_{θ}) − c_{λ_{T} −γ}s_{θ}s_{φ}) = (s_{λ_{T} −γ}(c_{φ}(c_{ζ} − 1)s_{θ} − s_{ζ}c_{θ}) + s_{λ_{T} −φ−γ}s_{θ}). (A4.1) A4.1

Define canonically

\begin{matrix}(μ, γ) = ( \begin{matrix}1 & 1 \\ 0 & 1 \end{matrix} ) ( \begin{matrix}φ \\ γ \end{matrix} ) ,  (A, T) = ( \begin{matrix}1 & 0 \\ −1 & 1 \end{matrix} ) ( \begin{matrix}A \\ K \end{matrix} ) , \\ p = \sqrt{−2T} cos γ,  q = \sqrt{−2T} sin γ, \end{matrix}  (A4.2) A4.2

and

\begin{matrix}(ξ, μ) = ( \begin{matrix}1 & 1 \\ 0 & 1 \end{matrix} ) ( \begin{matrix}ψ \\ μ \end{matrix} ) ,  (L, G) = ( \begin{matrix}1 & 0 \\ −1 & 1 \end{matrix} ) ( \begin{matrix}L \\ A \end{matrix} ) , \\ π = \sqrt{2G} cos μ,  κ = \sqrt{2G} sin μ. \end{matrix}  (A4.3) A4.3

If we set s_{θ} = (2G/L)^{\frac{1}{2}} σ_{θ}, s_{ζ} = (−2T/L)^{\frac{1}{2}} σ_{z}, thus defining σ_{z}, σ_{θ}, and use c_{ζ} −1 ≡ −\frac{1}{2} \frac{(−2T)}{L}c_{θ}, c_{θ}− 1 ≡ −\frac{1}{2} \frac{(2G)}{L} then cos α becomes

\frac{σ_{θ}}{L^{\frac{1}{2}}} (s_{λ_{T}} π − c_{λ_{T}} κ) − \frac{σ_{ζ}c_{θ}}{L^{\frac{1}{2}}} (s_{λ_{T}} p − c_{λ_{T}} q) − \frac{σ_{θ}c_{θ}}{2L^{\frac{1}{2}}} \frac{(pπ + qκ)}{L}(s_{λ_{T}} p − c_{λ_{T}} q), (A4.4) A4.4

and, since 2G = π^{2} + κ^{2} and −2T = p^{2} + q^{2} and therefore σ_{θ}, σ_{ζ}, c_{θ}, c_{ζ} are analytic at the origin and have value 1 there we have achieved complete rationalization of the Hamiltonian for T, G close to 0. Since there is no singularity at 0 a simplified Hamiltonian could be

H = \frac{L}{2I} (π^{2} + κ^{2}) + \frac{3Jω^{2}(λ_{T})}{2L} δ_{1}((π − p) s_{λ_{T}} − (κ − q) c_{λ_{T}} )^{2}. (A4.5) A4.5

Note that L is rigorously a constant of motion.

The equation is not useful to study the tidal capture into resonance phenomena in heavenly bodies because the evolution takes place on a time scale proportional to δ_{1}^{−1}i^{−1} if i is the tilt angle: the latter is very small in most cases (except perhaps the Moon and a few other large satellites) and of the order of δ_{1} so that the effects that we study become relevant on time scales much longer than the ones due to even slight equatorial asymmetry. Therefore the interesting model for our cases is the model with A ≃ L ≃ K and equatorial asymmetry (i.e. I_{y} < I_{x}): this is much harder and is discussed below.

p.A4.4 A4.2. Rationalization: L ≃ A ≠ K symmetric case. In this case define canonically

\begin{matrix}(ξ, φ) = ( \begin{matrix}1 & 1 \\ 0 & 1 \end{matrix} ) ( \begin{matrix}ψ \\ φ \end{matrix} ) ,  (L, G) = ( \begin{matrix}1 & 0 \\ −1 & 1 \end{matrix} ) ( \begin{matrix}L \\ A \end{matrix} ) , \\ π = \sqrt{2G} cos φ,  κ = \sqrt{2G} sin φ, \end{matrix}  (A4.6) A4.6

and cos α becomes

s_{λ_{T} −γ}(π \frac{c_{ζ}σ_{θ}}{L^{\frac{1}{2}}} − s_{ζ}c_{θ}) − c_{λ_{T} −γ} \frac{σ_{θ}}{L^{\frac{1}{2}}} κ. (A4.7) A4.7

And calling f_{0}(π, κ, L, K, γ) the function in (A4.7) and H_{0} = \frac{1}{2I} (L+(\frac{1}{2}(π^{2}+κ^{2}))^{2}) the Hamiltonian H = H_{0} + \frac{3}{2}δ_{1} Jω^{2}(λ_{T})f_{0}^{2} can be simplified for G ≃ 0 into

H = \frac{L}{2I} (π^{2} + κ^{2}) + \frac{3}{2}δ_{1} J ω^{2}(λ_{T})(s_{λ_{T} −γ}( − s_{ζ} + \frac{c_{ζ}}{L^{\frac{1}{2}}} π) − c_{λ_{T} −γ} \frac{1}{L^{\frac{1}{2}}} κ)^{2}. (A4.8) A4.8

The equations above can be used to study the corrections to the cruder approximation consisting in considering the Hamiltonian \frac{1}{2J} A^{2} + \frac{3}{2}δ_{1}Jω(λ_{T})^{2}(s_{λ_{T} −γ}s_{ζ})^{2} (i.e. neglecting s_{θ} = (1 − L^{2}/A^{2})^{\frac{1}{2}} ) and subsequently assuming that the adiabatic invariant K is constant. Since s_{ζ} = (1 − K^{2}/A^{2})^{\frac{1}{2}} this leads to the nodal precession model γ̇ = 3δ_{1} \frac{ω(λ_{T})^{2}}{ω_{D}} c_{θ}s_{λ_{T} −γ}^{2} which is a model used in simple theories of the lunar node precession (with ζ identified with the tilt angle, i.e. c_{ζ} = \frac{K}{A} of the symmetry axis over the normal to the orbit and A/J identified with the rotation velocity) which originated the theory of the rotation numbers by Poincaré, see [14], [13].

p.A4.4a A4.3. Rationalization: L ≃ A ≠ K asymmetric case. The nonsymmetric case, which is the only one that interests us in our applications, requires understanding the rationalization of sin α cos β. Algebraic analysis shows that H will take the form, if G ≡ \frac{1}{2}(π^{2} + κ^{2}) and ψ ≡ μ − φ, so that G cos^{2} ψ ≡ \frac{1}{2}(π cos μ + κ sin μ)^{2},

\begin{aligned}H = & \frac{1}{2I_{x}} (L + \frac{1}{2}(π^{2} + κ^{2}))^{2} − δ′ \frac{L^{2}}{2I_{z}} − δ_{2} \frac{I_{z}}{I_{x}I_{y}} (LG + \frac{1}{2}G^{2}) cos^{2} ψ+ \\ & + \frac{3}{2}Jω^{2}(λ_{T}) · (δ_{1}f_{0} + δ_{2}f_{1}), \end{aligned} (A4.9) A4.9

where f_{0}, f_{1} are analytic (see below) in their arguments (p, q), (K, γ), (L, ξ) and nonsingular near π, κ = 0, 0 ≤ K < L. To check the statement we express sin α cos β from the second of (A3.7). Setting \overline{cot}_{\overline{φ}−γ} = (c_{φ}c_{ζ}s_{θ} − s_{ζ}c_{θ}), so that \overline{cot}_{\overline{φ}−γ} = cot(\overline{φ} − γ) s_{θ}s_{φ}, and \overline{cot}_{\overline{ψ}−ψ} = (c_{φ}c_{θ}s_{ζ} − s_{θ}c_{ζ}), so that \overline{cot}_{\overline{ψ}−ψ} = cot(\overline{ψ} − ψ) s_{ζ}s_{φ}, we find

sin α cos β = \frac{1}{s2 ; \overline{θ}} ((\overline{cot}_{\overline{φ}−γ}c_{λ_{T} −γ} + s_{λ_{T} −γ}s_{θ}s_{φ})(\overline{cot}_{\overline{ψ}−ψ}c_{ψ} − s_{ψ}s_{ζ}s_{φ})+

+ c\overline{θ}(\overline{cot}_{\overline{φ}−γ}s_{λ_{T} −γ} − c_{λ_{T} −γ}s_{θ}s_{φ})(\overline{cot}_{\overline{ψ}−ψ}s_{ψ} + c_{ψ}s_{ζ}s_{φ})), (A4.10) A4.10

which is rewritten as

\begin{matrix}\frac{1}{s2 ; \overline{θ}} (((c_{φ}c_{ζ}s_{θ} − s_{ζ}c_{θ})c_{λ_{T} −γ} + s_{λ_{T} −γ}s_{θ}s_{φ})((c_{φ}c_{θ}s_{ζ} − s_{θ}c_{ζ})c_{ψ} − s_{ψ}s_{ζ}s_{φ}) \\ + c\overline{θ}(((c_{φ}c_{ζ}s_{θ} − s_{ζ}c_{θ})s_{λ_{T} −γ} − c_{λ_{T} −γ}s_{θ}s_{φ})((c_{φ}c_{θ}s_{ζ} − s_{θ}c_{ζ})s_{ψ} + c_{ψ}s_{ζ}s_{φ}))), \end{matrix} (A4.11) A4.11

which, by (A3.9), is (sin \overline{θ})^{−2}(a + b cos \overline{θ}), with

\begin{matrix}a = & ((c_{φ}(c_{ζ} − 1)s_{θ} − s_{ζ}c_{θ})c_{λ_{T} −γ} + c_{λ_{T} −γ−φ}s_{θ})((c_{φ}(c_{θ} − 1)s_{ζ} − s_{θ}c_{ζ})c_{ξ−φ} + c_{ξ}s_{ζ}), \\ b = & ((c_{φ}(c_{ζ} − 1)s_{θ} − s_{ζ}c_{θ})s_{λ_{T} −γ} + s_{λ_{T} −γ−φ}s_{θ})((c_{φ}(c_{θ} − 1)s_{ζ} − s_{θ}c_{ζ})s_{ξ−φ} + s_{ξ}s_{ζ}), \end{matrix} (A4.12) A4.12

where ξ = ψ + φ according to (A4.6). Furthermore cos α is given by (A4.7) and cos \overline{θ}, sin \overline{θ} are

cos \overline{θ} = (1 − s_{θ}^{2}s_{φ}^{2}/s_{\overline{φ}−γ}^{2})^{\frac{1}{2}}
= (1 − s_{θ}^{2}s_{φ}^{2}(1 + cot^{2}(\overline{φ} − γ))) = c_{θ}c_{ζ} + s_{θ}s_{ζ}c_{φ},
(sin \overline{θ})^{2} = s_{ζ}^{2} + s_{θ}^{2} − (1 + c_{φ}^{2}) s_{θ}^{2}s_{ζ}^{2} − 2s_{θ}c_{θ}s_{ζ}c_{ζ}c_{φ},

(A4.13) A4.13

and (A4.13) show analyticity because only c_{φ}, s_{φ} appear multiplied by the appropriate power of (2G)^{\frac{1}{2}} . Collecting the above algebraic relations we get

sin α cos β =c_{λ_{T} −γ−ξ} + L^{− \frac{1}{2}} (C_{1}π + C_{2}κ) + O(π^{2} + κ^{2}),
cos α =s_{ζ}s_{λ_{T} −γ} − c_{ζ}L^{− \frac{1}{2}} (πs_{λ_{T} −γ} − κ c_{λ_{T} −γ}) + O(π^{2} + κ^{2}),

(A4.14) A4.14

where C_{1}, C_{2} are (simple) analytic functions of ξ, λ_{T} − γ and analytic terms of second order in p, q are neglected, while no nonanalytic terms appear (to any order in π, κ).

The squares of the r.h.s.of the second in (A4.14) is f_{0} and the r.h.s. of the first equation is f_{1}. f_{0} in (A4.9). The above analysis requires |K|/L < 1 and it is not uniform as K → ±L.

p.A4.5 A4.4. Rationalization: L ≃ A ≃ K asymmetric case. In this case we canonically set

( \begin{matrix}μ \\ γ \end{matrix} ) = ( \begin{matrix}1 & 1 \\ 0 & 1 \end{matrix} ) ( \begin{matrix}φ \\ γ \end{matrix} ) ,  ( \begin{matrix}A \\ T \end{matrix} ) = ( \begin{matrix}1 & 0 \\ −1 & 1 \end{matrix} ) ( \begin{matrix}A \\ K \end{matrix} ) ,  (A4.15) A4.15

and then

( \begin{matrix}ξ \\ μ \end{matrix} ) = ( \begin{matrix}1 & 1 \\ 0 & 1 \end{matrix} ) ( \begin{matrix}ψ \\ μ \end{matrix} ) ,  ( \begin{matrix}L \\ G \end{matrix} ) = ( \begin{matrix}1 & 0 \\ −1 & 1 \end{matrix} ) ( \begin{matrix}L \\ A \end{matrix} ) ,  (A4.16) A4.16

so that (T, γ), (L, ξ), (G, μ) are canonically conjugated variables:

\begin{matrix}T = K − A, & G = A − L,  L = L, \\ γ = γ, & μ = φ + γ,  ξ = φ + γ + ψ. \end{matrix}  (A4.17) A4.17

We define canonically also the pairs (p, q) and (π, κ)

\begin{matrix}p = & (−2T)^{\frac{1}{2}} cos γ, \\ q = & (−2T)^{\frac{1}{2}} sin γ, \end{matrix}  \begin{matrix}π = & (2G)^{\frac{1}{2}} cos μ, \\ κ = & (2G)^{\frac{1}{2}} sin μ. \end{matrix} (A4.18) A4.18

The rationalization, i.e. the absence of the square roots of p^{2} +q^{2} and of π^{2} +κ^{2} in the remainders, is seen by going back to the original expression for sin α cos β in (A4.12) and change the name of the variables to adapt it to the new definitions in (A4.19).

In this case, however, we cannot achieve an analytic expression: the elimination of the square roots is not sufficient because the division by (sin \overline{θ})^{2} introduces denominators which vanish at the origin. Even though whenever the denominators vanish also the numerators do (because sin α cos β is bounded) a singularity of the generic type xy/(x^{2} + y^{2}) cannot be excluded: and therefore the leading term in the Hamiltonian might not necessarily be the leading term in the equations of motion and this makes perturbation analysis very difficult.

Let R(ρ) be a rotation by an angle ρ in a plane, then sin α cos β, by (A3.7), becomes the scalar product

R(ψ) ( \begin{matrix}c_{\overline{ψ}−ψ} \\ s_{\overline{ψ}−ψ} \end{matrix} ) · ( \begin{matrix}1 & 0 \\ 0 & cos \overline{θ} \end{matrix} ) R(λ_{T} − γ) ( \begin{matrix}c_{\overline{φ}−γ} \\ −s_{\overline{φ}−γ} \end{matrix} ) . (A4.19) A4.19

In fact the expression s−2 ; \overline{θ}(c\overline{θ} − 1) is analytic near the origin so that we only have to study the scalar product in (A4.19) with c\overline{θ} replaced by 1. This is

−s−2 ; \overline{θ} ( \begin{matrix}\overline{cot} & \overline{ψ}−ψ \\ s_{ζ}s_{φ} \end{matrix} ) · R(λ_{T} − γ − ψ) ( \begin{matrix}\overline{cot} & \overline{φ}−γ \\ −s_{θ}s_{φ}, \end{matrix} ) (A4.20) A4.20

where it is should be noted that s\overline{θ} is the common value of the length of the vectors. The latter scalar product is

−s−2 ; \overline{θ}R(−φ) ( \begin{matrix}c_{φ}c_{θ}s_{ζ} − s_{θ}c_{ζ} \\ s_{ζ}s_{φ} \end{matrix} ) · R(λ_{T} − ξ) ( \begin{matrix}c_{φ}c_{ζ}s_{θ} − s_{ζ}c_{θ} \\ −s_{θ}s_{φ} \end{matrix} ) (A4.21) A4.21

having applied a rotation R(−φ) to both sides of the scalar product.

Let us denote A, B and C, D the components of the two vectors:

\begin{matrix}A = & c_{φ} (c_{φ}c_{θ}s_{ζ} − s_{θ}c_{ζ}) + s_{φ}s_{ζ}s_{φ}, \\ B = & − s_{φ} (c_{φ}c_{θ}s_{ζ} − s_{θ}c_{ζ}) + c_{φ}s_{ζ}s_{φ}, \\ C = & (c_{φ}c_{ζ}s_{θ} − s_{ζ}c_{θ}), \\ D = & − s_{θ}s_{φ}, \end{matrix} (A4.22) A4.22

and one has, by our construction and by the unitarity of the rotations, s2 ; \overline{θ} = A^{2} + B^{2} = C^{2} + D^{2}. Therefore sin α cos β is such that

− sin α cos β = c_{λ_{T} −ξ} \frac{2(AC + BD)}{A^{2} + B^{2} + C^{2} + D^{2}} + s_{λ_{T} −ξ} \frac{2(BC − AD)}{A^{2} + B^{2} + C^{2} + D^{2}} , (A4.23) A4.23

where the ratios depend only on π, κ, p, q, L, ξ but they are λ_{T} -independent; furthermore by the above analysis the combinations of the functions A, B, C, D appearing in the numerators and denominators of the above expression are analytic at the origin. One has also

\begin{matrix}\frac{2(AC + BD)}{A^{2} + B^{2} + C^{2} + D^{2}} = −1 + Ω_{2}(p, q, π, κ), \\ \frac{2(BC − AD)}{A^{2} + B^{2} + C^{2} + D^{2}} = Ω′_{2}(p, q, π, κ), \end{matrix} (A4.24) A4.24

where Ω_{2}, Ω′_{2} denote functions analytic outside the region \overline{θ} = 0 (which coincides with A^{2} + B^{2} + C^{2} + D^{2} = 0) and with bounded second derivatives, possibly discontinuous on the surface \overline{θ} = 0. That the leading term near the origin is given by the r.h.s. of (A4.24) can realized by staring at figure 4 above and remarking that \overline{θ} → 0 implies φ →, θ − ζ → 0, (\overline{φ} − γ) − (\overline{ψ} − ψ) → 0. That the corrections are of “second order” in the sense of (A4.24) follows from the fact that the analyticity of the expressions like (A4.12) gives as a result a low degree polynomial in p, π, q, κ with coefficients which are analytic in T/L and G/L, i.e. in p^{2} + q^{2}, π^{2} + κ^{2}.

A straightforward calculation yields, in fact, the values of the quantities in (A4.24). For instance introducing the further abbreviation s_{x}^{2}S_{x}^{2} =^{def}(c_{x} − 1 + \frac{1}{2}s_{x}^{2}) the coefficients of c_{λ_{T} −ξ} and s_{λ_{T} −ξ} in (A4.23) are, respectively, the (manifestly rational) expressions

\frac{AC + BD}{C^{2} + D^{2}} = − (1 + \frac{s_{φ}^{2}s_{θ}^{2}s_{ζ}^{2}}{s2 ; \overline{θ}} \frac{(1 − c_{ζ})(1 − c_{θ})}{s_{θ}^{2}s_{ζ}^{2}}(s_{θ}s_{ζ}c_{φ} + \frac{s_{ζ}^{2}s_{θ}^{2}}{(1 − c_{ζ})(1 − c_{θ})}(S_{θ}^{2} + S_{ζ}^{2}))), (A4.25) A4.25
\frac{BC − AD}{C^{2} + D^{2}} =\frac{1}{s2 ; \overline{θ}} (s_{φ}s_{θ}s_{ζ}(1 − c_{ζ}c_{θ}) − s_{φ}c_{φ}(s_{θ}^{2}s_{ζ}^{2}(1 + S_{ζ}^{2} + S_{θ}^{2}) + c_{φ}s_{θ}s_{ζ}(c_{θ} − 1)(c_{ζ} − 1))),

which has the form in (A4.24). This can be seen to be not analytic at the origin: for instance the second derivative of the first in (A4.22) in the direction p_{0}, q_{0}, π_{0}, κ_{0} evaluated at the origin is 0 if (p_{0} − π_{0})^{2} + (q_{0} − κ_{0})^{2} ≠ 0 but on the latter 2–dimensional plane it is bounded and in general not 0. Therefore the Laplacian is bounded and it has a bounded normal limit to the singularity surface and the function in (A4.20) is of class C^{c} for all c < 2 but it is not more regular (in particular it is not even C^{2}): roughly in the variables x = (p − π)a, y = (q − κ)b, z = (p + π)c, w = (q + κ)d it is a sum of terms similar to the following

\frac{(xy)^{2}(x^{2} + y^{2} + z^{2} + w^{2})}{x^{2} + y^{2} + (x^{2} + y^{2})(z^{2} + w^{2})} (A4.26) A4.26

if the constants a, b, c, d > 0 are suitably defined.

Note that the (A4.24) implies that the corrections to the equations of motion due to the terms of higher order in π, κ, p, q are also small in the latter quantities and of almost the first order.

Retaining only the lowest orders a simplified Hamiltonian is

H = \frac{1}{2J} L^{2} + \frac{3}{2}Jω_{T}^{2} (λ_{T})δ_{2} cos^{2}(λ_{T} − ξ), (A4.27) A4.27

with δ_{2} = (I_{y} − I_{x})/I_{z}. Therefore the equations of motion are rather simple in the considered approximation and become ξ̇ = L/J and L̇ = −\frac{3}{2}Jω_{T}^{2} (λ_{T})δ_{2}∂_{ξ} cos^{2}(λ_{T} − ξ) because L, ξ are pairs of conjugated coordinates. Or

ξ̈ = −\frac{3}{2} \frac{ω^{2}(λ_{T})}{ω_{T}^{2}}δ_{2} \frac{∂}{∂ξ} cos^{2}(λ_{T} − ξ), (A4.28) A4.28

in dimensionless units (i.e. time is measured in units of ω_{T}^{−1}). Expanding λ_{T} in powers of the eccentricity e via Kepler’s law λ = (1 − e^{2})^{\frac{3}{2}} ∫_{0}^{λ_{T}} (1 − e cos ζ)^{−2}dζ and setting λ = t (as in our units the average anomaly of the planet and the time coincide) then, once also the friction contribution is taken into account, (A4.26) gives (5.3) if the expansion is truncated to the fifth degree in e. The coefficients a_{j} in (4.4) are computed in table 4.1.

The singularity that remains at the origin of the new variables is much weaker that the one in the original variables (i.e. it is algebraic rather than of square root type and it vanishes to order almost 2). This gives some grounds to argue in favor of the statement that the model might be a good approximation as its widespread use in Celestial Mechanics shows. However the rigorous perturbation analysis necessary to substantiate more quantitatively the statement would pose serious problems. For instance C^{2}–regularity in a three and half degrees of freedom, like the asymmetric model considered here, is not enough (in general) to prove any type of KAM stability result.

Finally it should be stressed that, of course, the (A4.21) are the same equations that one would obtain by studying an asymmetric rigid body constrained to rotate around an axis which moves orthogonally to the orbit plane and intersecting it on a Keplerian elliptic orbit. The point of our analysis has been to check in which cases no nonanalytic terms arise in the corrections of higher orded in s_{ζ}, s_{θ} in a proper system of coordinates.

Appendix A5. Some Generalizations app.A5

p.A5.1 A5.1. More general forms for the friction. In the applications to Celestial Mechanics the friction takes often rather elaborated forms: therefore it is of some interest to study how strongly the analysis in this paper depends on the simple forms assumed for the friction.

So (2.1) will be replaced by

θ̈ + εG(θ, t) + γθ̇r(θ, θ̇, t) + γ′ = 0, (A5.1) A5.1

where r(θ, θ̇, t) > 0 is 2π-periodic in (θ, t), analytic in each variable, and ε, γ, γ′ are parameters. For instance the tidal torque should be proportional to the sine of the phase shift τ and τ is proportional to ξ̇ − 1 only for ξ̇ small. Hence (A5.1) provides flexibility for a rather general class of interesting friction models and it is desirable to extend the theory t such cases.

For γ = γ′ = 0 we can still derive (A5.1) as the Hamilton equations of the system described by the Hamiltonian (2.2). For γ, γ′ ≠ 0, (A5.1) can be seen as a Hamiltonian system with friction given by the equations of motion ^{5}

\{ \begin{matrix}θ̇ = ω + Θ, \\ ṫ = 1, \end{matrix}  \{ \begin{matrix}Θ̇ = −ε∂_{θ}g(θ, t) − γ(Θ + ω) r(θ, Θ + ω, t) − γ′, \\ Ṫ = −ε∂_{t}g(θ, t). \end{matrix}  (A5.2) A5.2

Likewise (2.5) is replaced by

\{ \begin{matrix}α̇ = 1/q + mP, \\ β̇ = −qP, \end{matrix} \{ \begin{matrix}Ȧ = −ε∂_{α}f(α, β) − pγP u(α, β, P) − p^{2}γ/q u(α, β, P) − pγ′, \\ Ḃ = −ε∂_{β}f(α, β) − nγP u(α, β, P) − npγ/q u(α, β, P) − nγ′ & , \end{matrix} (A5.3) A5.3

where f(α, β) = g(θ(α, β), t(α, β)) as in (2.5), and

\begin{matrix}P ≡ Θ = mA − qB = (α̇ − 1/q)/m = −β̇/q, \\ u(α, β, P) = r(θ(α, β), θ̇(A, B), t(α, β)) = r(pα + nβ, p/q + P, qα + mβ); \end{matrix} (A5.4) A5.4

for γ = γ′ = 0 the corresponding Hamiltonian is given by (2.6), and, for ε = 0, the periodic solution X(t) is still transformed into X_{0}(t) = CX(t) = (α_{0} + t/q, β_{0}, 0, 0), with qα_{0} + mβ_{0} = 0.

In terms of the only (α, β) variables, by using that

Ṗ = −ε (m∂_{α}f − q∂_{β}f) − γP u − \frac{p}{q}γu − γ′, (A5.5) A5.5

we can rewrite (A5.3) as

\{ \begin{matrix}α̈ = −mε (m∂_{α}f − q∂_{β}f) − γα̇u(α, β, (α̇ − 1/q)/m) − nγu(α, β, (α̇ − 1/q)/m) − mγ′, \\ β̈ = qε (m∂_{α}f − q∂_{β}f) − γβ̇u(α, β, −β̇/q) + pγu(α, β, −β̇/q) + qγ′. \end{matrix} (A5.6) A5.6

As in Section 2 it will be convenient to write

f(α, β) = ∑_{ν∈ℤ} e^{iνα}f_{ν}(β), u(α, β, P) = ∑_{N=0}^{∞} P^{N}u_{N}(α, β) = ∑_{N=0}^{∞} P^{N} ∑_{ν∈ℤ} e^{iνα}u_{N,ν}(β), (A5.7) A5.7

where, for all ν ∈ ℤ and N ∈ ℤ_{+}, the coefficients f_{ν}(β) and u_{N,ν}(β) are 2π-periodic in β.

We want to study the possibility of continuing some of the periodic solutions X_{0}(t) when ε ≠ 0 and γ, γ′ ≠ 0.

We shall fix γ = Cε and γ′ = Cε, with C and C′ parameters to be varied, and we shall look for a solution which is analytic in ε for ε small enough. This means that we shall write, formally, α(t) and β(t) as in (2.9), with the functions a(ψ, β; ε) and b(ψ, β; ε) given by (2.10).

5 It would be also possible to generalize such equations by considering a potential g depending also on T, Θ and by adding friction terms to all components of the vector field appearing in (A5.2).

Observe that a solution (2.9) of (A5.6) is not necessarily a periodic orbit of (A5.3). In fact α̇ and β̇ do not determine A and B but only P = mA − qB. In term of the original variables this is due to the fact that θ, Θ and t are independent of T, so that (2.9) describes a periodic orbit only if the friction terms do not work in the average; but this would impose conditions on C and C′.

For the model (A5.1) a result analogous to theorem 2.3 holds. More precisely one has the following result (reducing to theorem 2.3 for C′ = 0 and u(α, β, P) = 1). The proof is entirely analogous to the one in Section 3.

p.A5.2 A5.2. Theorem. Fix ω = p/q. If C and C′ are so fixed that there exists β_{0} such that

∂_{β}f_{0}(β_{0}) = \frac{p}{q^{2}} Cu_{0,0}(β_{0}) + \frac{1}{q}C′,  ∂_{β}^{2}f_{0}(β_{0}) − \frac{p}{q^{2}} C∂_{β}u_{0,0}(β_{0}) ≠ 0,  (A5.8) A5.8

then for ε small enough there is a periodic solution of the equation of motions of the form (2.9).

p.A5.3 A5.3. Remarks. (1) Even if ∂_{β}f_{0}(β) is identically vanishing the condition (A5.8) can be still verified, provided that one chooses C and C′ such that

\frac{p}{q^{2}} Cu_{0,0}(β_{0}) + \frac{1}{q}C′ = 0,  C∂_{β}u_{0,0}(β_{0}) ≠ 0;  (A5.9) A5.9

in particular the function u_{0,0}(β) cannot be identically constant. Note however that such a choice of the constants C and C′ has no particular physical meaning, and it has to be considered just as a mathematical curiosity.

(2) If we are interested in positive values for C and C′, and u_{0,0}(β) is a function weakly varying around its positive average value, then a (reasonable) sufficient condition for the first equation in (A5.8) to be satisfied is to require that one has

\frac{p}{q^{2}} C max_{β∈[0,2π]} u_{0,0}(β) + \frac{1}{q}C′ ≤ max_{β∈[0,2π]} ∂_{β}f_{0}(β), (A5.10) A5.10

note that such a condition reduces to (2.11) for r(θ, θ̇, t) ≡ 1 (so that u_{0}(β) ≡ 1) and γ′ = 1 (so that C′ ≡ 0) in (A5.1).

(3) Of course we could also consider if and how theorems 2.6 and (2.9) extend to the more general case: we leave their formulation to the reader.

References

B[1] G.D. Birkhoff, Dynamical systems, American Mathematical Society Colloquium Publications Vol. IX, Providence, R.I. 1966.

BKM[2] P.I. Bakulin, E.V. Kononovich, V.I. Moroz, Handbook of general astronomy, Nauka, Moscow, 1983.

BFG[3] G.C. Benettin, F. Fassò, M. Guzzo, , preprint, 2003.

BGG[4] M.V. Bartuccelli, G. Gentile, K.V. Georgiou, On the dynamics of a vertically driven damped planar pendulum, R. Soc. Lond. Proc. Ser. A Math. Phys. Eng. Sci. 457 (2001), no. 2016, 3007–3022.

C[5] A. Celletti, Analysis of resonances in the spin-orbit problem in celestial mechanics: the synchronous resonance (part I), Z. Angew. Math. Phys. 41 (1990), no. 2, 174–204. Analysis of resonances in the spin-orbit problem in celestial mechanics: higher order resonances and some numerical experiments (part II), Z. Angew. Math. Phys. 41 (1990), no. 4, 453–479.

CC[6] A. Celletti, L. Chierchia, Construction of stable periodic orbits for the spin-orbit problem of celestial mechanics, Regul. Chaotic Dyn. 3 (1998), no. 3, 107–121.

G1[7] G. Gallavotti, The elements of mechanics, Springer, Berlin, 1983.
G2[8] G. Gallavotti, Foundations of fluid dynamics, Springer, Berlin, 2002.

GG[9] G. Gallavotti, G. Gentile, Hyperbolic low-dimensional invariant tori and summation of divergent series, Comm. Math. Phys. 227 (2002), no. 3, 421–460.

GP[10] P. Goldreich, S. Peale, Spin-orbit coupling in the solar system, Astronom. J. 71 (1966), no. 6, 425–438.

BJ[11] T.J. Burns, C.K.R.T. Jones, A mechanism for capture into resonances, Phys. D 69 (1993), no. 1–2, 85–106.

LN[12] J. Laskar, O. Néron de Surgy, On the long term evolution of the spin of the Earth, Astronom. and Astrophys. 318 (1997), 975–989.

Le[13] E.E. Levi, Sur les équations différentielles périodiques, C. R. Acad. Sci. Paris Sér. A-B, 153 (1911). Reprinted in Opere, Vol. 1, p. 224, Cremonese, Roma, 1959.

LC[14] T. Levi-Civita, Sur les équations linéaires à coefficients periodiques et sur le mouvement moyen du noeud lunaire, Ann. Sci. École Norm. Sup. 28 (1911), 325–376. Reprinted in Opere Matematiche. Memorie e note, Vol. 3, p. 205, Zanichelli, Bologna, 1956.

N[15] A.I. Neishtadt, On probabilistic phenomena in perturbed systems, Selecta Math. Sov. 12 (1993), no. 3, 195–210.

PMC[16] C.C. Plummer, D. McGeary, D.H. Carlson, Physical geology, McGraw-Hill, 2003.

Po[17] H. Poincaré, Sur les courbes définies par les équations differentielles, J. Math. Pures Appl. 1 (1885), 167–244. Reprinted in Oeuvres Complètes, Vol. 1, paper no. 57, p. 156, Gauthiers–Villars, Paris, 1934–1954, and in Mémoire sur les courbes définies par une équation differentielle, Ch. XV, p. 243, Gabay, Paris, 1993.

R[18] G. Ranalli, Rheology of the Earth, Chapman & Hall, London, 1987.

W[19] M.M. Woolfson, The origin and evolution of solar system, Graduate series in Astronomy, Institute of Physics Publishing, Bristol, 2000.

Momenti d’inerzia di un ellissoide con assi c < b < a e densità ρ:

C = \frac{ρV}{5} (a^{2} + b^{2}), B = \frac{ρV}{5} (c^{2} + a^{2}), A = \frac{ρV}{5} (c^{2} + b^{2})

where V = \frac{4π}{3} a b c and ρ is the density. It is C > B > A.

Assi solidali si ottengono applicando le rotazioni R_{z}(\overline{φ})R_{x}(\overline{θ})R_{z}(\overline{ψ}) agli assi fissi. Si ha

( \begin{matrix}cos \overline{φ} & − sin \overline{φ} & 0 \\ sin \overline{φ} & cos \overline{φ} & 0 \\ 0 & 0 & 1 \end{matrix} ) ( \begin{matrix}1 & 0 & 0 \\ 0 & cos \overline{θ} & − sin \overline{θ} \\ 0 & sin \overline{θ} & cos \overline{θ} \end{matrix} ) ( \begin{matrix}cos \overline{ψ} & − sin \overline{ψ} & 0 \\ sin \overline{ψ} & cos \overline{ψ} & 0 \\ 0 & 0 & 1 \end{matrix} )

e quindi

(cos \overline{φ} cos \overline{ψ} − sin \overline{φ} sin \overline{ψ} cos \overline{θ}) (− cos \overline{φ} sin \overline{ψ} − sin \overline{φ} cos \overline{ψ} cos \overline{θ}) sin \overline{θ} sin \overline{φ}
(sin \overline{φ} cos \overline{ψ} + cos \overline{φ} sin \overline{ψ} cos \overline{θ}) (− sin \overline{φ} sin \overline{ψ} + cos \overline{φ} cos \overline{ψ} cos \overline{θ}) − sin \overline{θ} cos \overline{φ}
sin \overline{ψ} sin \overline{θ} cos \overline{ψ} sin \overline{θ} cos \overline{θ}

e le colonne sono i vettori i_{1}, i_{2}, i_{3}. Trigonometria sferica: Sia (A, B, C, α, β, γ) un triangolo sferico; allora

cos α = − cos β cos γ + sin β sin γ cos A

Posto: P =^{def} \frac{1}{2}(A + B + C) si trova

sin \frac{1}{2}α =( sin(P − B) sin(P − C)/(sin B sin C))^{\frac{1}{2}}
cos \frac{1}{2}α =( sin P sin(P − A)/(sin B sin C))^{\frac{1}{2}}
tan \frac{1}{2}α =( sin(P − B) sin(P − C)/(sin P sin(P − A)))^{\frac{1}{2}}

Se σ =^{def} \frac{1}{2}(α + β + γ)

sin \frac{1}{2}A =( − cos σ cos(σ − α)/(sin β sin γ))^{\frac{1}{2}}
cos \frac{1}{2}A =( cos(σ − β) cos(σ − γ)/(sin β sin γ))^{\frac{1}{2}}
tan \frac{1}{2}α =( − cos σ cos(σ − α)/(sin(σ − β) sin(σ − γ)))^{\frac{1}{2}}

Inoltre (regole di Nepero)

\frac{sin \frac{1}{2}(α + β)}{sin \frac{1}{2}(α − β)} = \frac{tan \frac{1}{2}C}{tan \frac{1}{2}(A − B)}, \frac{sin \frac{1}{2}(A + B)}{sin \frac{1}{2}(A − B)} = \frac{cot \frac{1}{2}γ}{tan \frac{1}{2}(α − β)}

\frac{cos \frac{1}{2}(α + β)}{cos \frac{1}{2}(α − β)} = \frac{tan \frac{1}{2}C}{tan \frac{1}{2}(A + B)}, \frac{cos \frac{1}{2}(A + B)}{cos \frac{1}{2}(A − B)} = \frac{cot \frac{1}{2}γ}{tan \frac{1}{2}(α + β)}

Infine se ε =^{def}α + β + γ − π l’area del triangolo è πε e

tan \frac{1}{4}ε = ( tan \frac{1}{2}P tan \frac{1}{2}(P − A) tan \frac{1}{2}(P − B)\frac{1}{2} tan(P − C))^{\frac{1}{2}}

19/novembre/2010; 10:26 33 Altre identità sin \frac{1}{2}(α + β) cos \frac{1}{2}C = cos \frac{1}{2}(A − B) cos \frac{1}{2}γ cos \frac{1}{2}(α + β) cos \frac{1}{2}C = cos \frac{1}{2}(A + B) sin \frac{1}{2}γ sin \frac{1}{2}(α − β) sin \frac{1}{2}C = sin \frac{1}{2}(A − B) cos \frac{1}{2}γ cos \frac{1}{2}(α − β) sin \frac{1}{2}C = sin \frac{1}{2}(A + B) sin \frac{1}{2}γ

tan \frac{1}{4}ε = \frac{sin \frac{1}{2}(α + β) − sin \frac{1}{2}(π − γ)}{cos \frac{1}{2}(α + β) + cos \frac{1}{2}(π − γ)}

e inoltre se

r =^{def}( sin(P − A) sin(P − B) sin(P − C)/ sin P )^{\frac{1}{2}}
ρ =^{def}(\frac{− cos σ}{cos(σ−α) cos(σ−β) cos(σ−γ)})^{\frac{1}{2}}

si trova

tan \frac{1}{2}α = \frac{r}{sin(P − A)}, tan \frac{1}{2}A = ρ cos(σ − α)

Da: [CS] Carmichael, R.D., Smith, E.R.: Plane and spherical trigonometry, Ginn, Boston, 1930. Teorema (dualità sferica) Se (α, β, γ; A, B, C) è un triangolo sferico anche (A, B, π−C; α, β, π−γ) e’ un triangolo sferico.