Periodic and quasi-periodic attractors for the spin-orbit evolution of Mercury with a realistic tidal torque

TL;DR

This paper investigates Mercury's spin-orbit resonance using a realistic tidal torque model, revealing that Mercury's librations are quasi-periodic and providing detailed numerical and analytical insights into its dynamics.

Contribution

It introduces a detailed analysis of Mercury's spin-orbit evolution with a realistic, non-smooth tidal torque model, highlighting the quasi-periodic nature of librations.

Findings

Mercury's librations in the 3:2 resonance are quasi-periodic.

The realistic tidal torque model affects the stability and nature of Mercury's spin states.

Numerical methods must carefully handle non-smooth functions for accurate simulations.

Abstract

Mercury is entrapped in a 3:2 resonance: it rotates on its axis three times for every two revolutions it makes around the Sun. It is generally accepted that this is due to the large value of the eccentricity of its orbit. However, the mathematical model originally introduced to study its spin-orbit evolution proved not to be entirely convincing, because of the expression commonly used for the tidal torque. Only recently, in a series of papers mainly by Efroimsky and Makarov, a different model for the tidal torque has been proposed, which has the advantages of being more realistic, and of providing a higher probability of capture in the 3:2 resonance with respect to the previous models. On the other hand, a drawback of the model is that the function describing the tidal torque is not smooth and consists of a superposition of kinks, so that both analytical and numerical computations turn out to be rather delicate: indeed, standard perturbation theory based on power series expansion cannot be applied and the implementation of a fast algorithm to integrate the equations of motion numerically requires a high degree of care. In this paper, we make a detailed study of the spin-orbit dynamics of Mercury, as predicted by the realistic model: In particular, we present numerical and analytical results about the nature of the librations of Mercury's spin in the 3:2 resonance. The results provide evidence that the librations are quasi-periodic in time.

Tables

Table 1. Table 4: The initial conditions for the existing periodic solutions and the corresponding eigenvalues. The eigenvalues λ 1 , 2 subscript 𝜆 1 2 \lambda_{1,2} for each periodic solution either form a conjugate pair, in which case | λ 1 , 2 | − 1 subscript 𝜆 1 2 1 |\lambda_{1,2}|-1 is given, or are both real, in which case both are given.

Resonance	$\theta \left(0\right)$	$\dot{\theta }\left(0\right)/n$	Eigenvalues
$-1:1$	0.36289190131044645472	$-1.00004242365858443089$	$2.024\times {10}^{-9}$
$-1:2$	0.14580421354300878946	$-0.50006849412399051400$	$2.783\times {10}^{-9}$
1:2	3.14129563170348761883	0.49980635331803679181	$0.9669,1.0342$
1:1	3.14151499384565687042	0.99986201340697665762	$-4.461\times {10}^{-4}$
3:2	3.14150380436395113505	1.50005973350740330252	$1.055\times {10}^{-4}$
2:1	$3.26027930307144126711\times {10}^{-5}$	2.00012557558534916792	$1.786\times {10}^{-3}$
5:2	3.14140519201664595044	2.50012075040501328073	$-3.628\times {10}^{-4}$
3:1	3.14109199137670843320	3.00009814397107114853	$-2.636\times {10}^{-5}$
7:2	3.14027640889440704126	3.50007711111008245662	$-3.835\times {10}^{-6}$
4:1	3.13797190712320535390	4.00006157245270746253	$-6.337\times {10}^{-7}$
$-1:1$	1.20792006104664609582	$-0.99995757575039987029$	$0.9992,1.0008$
$-1:2$	1.42500286020411552411	$-0.49993150584317966306$	$0.9986,1.0014$
1:2	1.57112385469851460569	0.50019364882055637631	$2.342\times {10}^{-6}$
1:1	1.57068938450889863242	1.00013792675908729505	$0.9048,1.1042$
3:2	1.57075984135159670901	1.49994030293249049891	$0.9185,1.0889$
2:1	1.57099968204819540739	1.99987444617026058657	$0.9433,1.0638$
5:2	1.57101812013673537458	2.49987925336853351100	$0.9613,1.0395$
3:1	1.57130265033260668261	2.99990185551468461907	$0.9760,1.0246$
7:2	1.57211353266178141100	3.49992288814339436128	$0.9854,1.0147$
4:1	1.57441706775605802984	3.99993842708145177608	$0.9914,1.0087$

Equations

$$C\ddot{\theta }={\text{mathpzc}T}_z^{(\mathrm{TRI})}+{\text{mathpzc}T}_z^{(\mathrm{TIDE})},$$

$$\ddot{\theta }=-\epsilon G(\theta ,t)-\epsilon \gamma F(\dot{\theta }),$$

$$\ddot{\theta }=-\zeta G(\theta ,t)-\eta F(\dot{\theta }),$$

$$G(\theta ,t)=\sum _{k\in \mathcal{K}}{A}_k\sin (2\theta -knt),\mathcal{K}=\{-2,-1,0,1,2,3,4,5,6,7,8\},$$

$$F(\dot{\theta }):=\sum _{k\in \mathcal{Q}}{A}_k^2\Xi ({\Omega }_k),\mathcal{Q}=\{1,2,3,4,5,6,7,8,9\},$$

$${\Omega }_k:=nk-2\dot{\theta },\Xi (\omega )=\mathrm{sgn}(\omega )\frac{\mathcal{I}(|\omega |)|\omega |}{(\mathcal{R}(|\omega |)+\mathcal{A}|\omega |{)}^2+{\mathcal{I}}^2(|\omega |)},$$

$$\mathcal{I}(\omega )$$

$$\mathcal{R}(\omega )$$

$$\zeta$$

$$\eta$$

$$\ddot{\theta }=-\epsilon G(\theta ,t)-\epsilon \gamma F(\dot{\theta }),$$

$$\theta (t)=\bar{\theta }+n{\omega }_{0}t+\bar{y}t-\epsilon {\int }_0^t\mathrm{d}\tau {\int }_0^{\tau }\mathrm{d}{\tau }^{\prime }\left[G(\theta ({\tau }^{\prime }),{\tau }^{\prime })+\gamma F(\dot{\theta }({\tau }^{\prime }))\right],$$

$${\theta }_0(t)={\bar{\theta }}_0+n{\omega }_{0}t,$$

$${\theta }_1(t)$$

$$⟨f⟩:=\frac{n}{2\pi q}{\int }_0^{2\pi q/n}\mathrm{d}tf(t).$$

$$u_1(t):=G({\bar{\theta }}_0+n{\omega }_{0}t,t)+\gamma F(n{\omega }_0),U_1^{\prime }(t):=u_1(t),{\mathcal{U}}_1^{\prime }(t):=U_1(t),$$

$${\int }_0^{\tau }\mathrm{d}{\tau }^{\prime }\left[G({\bar{\theta }}_0+n{\omega }_0{\tau }^{\prime },{\tau }^{\prime })+\gamma F(n{\omega }_0)\right]=U_1(\tau )-U_1(0)$$

$$\gamma F(n{\omega }_0)+\frac{n}{2\pi q}{\int }_0^{2\pi q/n}\mathrm{d}tG({\bar{\theta }}_0+n{\omega }_{0}t,t)=0,$$

$$A_{k_0}\sin 2{\bar{\theta }}_0=-\gamma F(n{\omega }_0),k_0=\frac{2p}{q}.$$

$$\bar{\theta}_{0}^{(1)}=\frac{1}{2}{\rm arcsin}\Biggl{(}-\frac{\gamma\,F(n\omega_{0})}{A_{k_{0}}}\Biggr{)}\in\left(-\frac{\pi}{4},\frac{\pi}{4}\right),\quad\bar{\theta}_{0}^{(2)}=\frac{\pi}{2}-\bar{\theta}_{0}^{(1)},\quad\bar{\theta}_{0}^{(3)}=\bar{\theta}^{(1)}-\pi,\quad\bar{\theta}_{0}^{(4)}=\bar{\theta}_{0}^{(2)}-\pi,$$

$$U_1(t)=-\sum _{k\in {\mathcal{K}}_0}\frac{A_k}{(2{\omega }_0-k)n}\cos (2{\bar{\theta }}_0+(2{\omega }_0-k)nt),$$

$${\mathcal{U}}_1(t)=-\sum _{k\in {\mathcal{K}}_0}\frac{A_k}{(2{\omega }_0-k{)}^2{n}^2}\sin (2{\bar{\theta }}_0+(2{\omega }_0-k)nt).$$

$${\theta }_1(t)={\bar{\Theta }}_1+n{\omega }_{0}t+\epsilon \sum _{k\in {\mathcal{K}}_0}\frac{A_k}{(2{\omega }_0-k{)}^2{n}^2}\sin (2{\bar{\theta }}_0+(2{\omega }_0-k)nt),$$

$${\theta }_2(t)={\bar{\theta }}_2+n{\omega }_{0}t+{\bar{y}}_{2}t-\epsilon {\int }_0^t\mathrm{d}\tau {\int }_0^{\tau }\mathrm{d}{\tau }^{\prime }\left[G({\theta }_1({\tau }^{\prime }),{\tau }^{\prime })+\gamma F({\dot{\theta }}_1({\tau }^{\prime }))\right],$$

$$\frac{n}{2\pi q}{\int }_0^{2\pi q/n}\mathrm{d}{\tau }^{\prime }\left[G({\bar{\Theta }}_1+n{\omega }_0{\tau }^{\prime }+{\xi }_1({\tau }^{\prime }),{\tau }^{\prime })+\gamma F(n{\omega }_0+{\dot{\xi }}_1({\tau }^{\prime }))\right]=0,$$

$$\mathcal{G}({\bar{\Theta }}_1):=\frac{n}{2\pi q}{\int }_0^{2\pi q/n}\mathrm{d}{\tau }^{\prime }G({\bar{\Theta }}_1+n{\omega }_0{\tau }^{\prime }+{\xi }_1({\tau }^{\prime }),{\tau }^{\prime })=-\frac{n}{2\pi q}{\int }_0^{2\pi q/n}\mathrm{d}{\tau }^{\prime }\gamma F(n{\omega }_0+{\dot{\xi }}_1({\tau }^{\prime })),$$

$$G({\bar{\Theta }}_1+n{\omega }_0{\tau }^{\prime }+{\xi }_1({\tau }^{\prime }),{\tau }^{\prime })=G({\bar{\Theta }}_1+n{\omega }_0{\tau }^{\prime },{\tau }^{\prime })+\Gamma ({\tau }^{\prime },\epsilon ),\Gamma ({\tau }^{\prime },\epsilon )={\partial }_{\theta }G({\bar{\Theta }}_1+n{\omega }_0{\tau }^{\prime },{\tau }^{\prime }){\xi }_1({\tau }^{\prime })+O({\epsilon }^2),$$

$$\mathcal{G}({\bar{\Theta }}_1)=\frac{n}{2\pi q}{\int }_0^{2\pi q/n}\mathrm{d}{\tau }^{\prime }G({\bar{\Theta }}_1+n{\omega }_0{\tau }^{\prime },{\tau }^{\prime })+O(\epsilon )=A_{k_0}\sin 2{\bar{\Theta }}_1+O(\epsilon ).$$

$$J(n{\omega }_0):=\frac{n}{2\pi q}{\int }_0^{2\pi q/n}\mathrm{d}{\tau }^{\prime }\gamma F(n{\omega }_0+{\dot{\xi }}_1({\tau }^{\prime })),$$

$$A_{k_0}\sin 2{\bar{\Theta }}_1=-J(n{\omega }_0),k_0=\frac{2p}{q},$$

$$\bar{\Theta}_{1}^{(1)}=\frac{1}{2}{\rm arcsin}\Biggl{(}-\frac{J(n\omega_{0})}{A_{k_{0}}}\Biggr{)}\in\left(-\frac{\pi}{4},\frac{\pi}{4}\right),\quad\bar{\Theta}_{1}^{(2)}=\frac{\pi}{2}-\bar{\Theta}_{1}^{(1)},\quad\bar{\Theta}_{1}^{(3)}=\tilde{\theta}^{(1)}-\pi,\quad\bar{\Theta}_{1}^{(4)}=\bar{\Theta}_{1}^{(2)}-\pi.$$

Full text

1 Introduction

The reason that Mercury is entrapped in the 3:2 resonance has been extensively investigated in the literature. It is commonly accepted that this is due to the high value of Mercury’s eccentricity ($0.2056$); however, there is no universal consensus about the mechanism by which the entrapment has occurred.

The Mercury-Sun system is usually studied as a satellite-planet system, with the satellite described as an ellipsoidal body orbiting around its primary in a Keplerian orbit. If $\theta$ denotes the sidereal angle, that is, the angle that the longest axis of the satellite forms with respect to the line of apsides of the orbit, the time evolution of $\theta$ is described by the second order ordinary differential equation

[TABLE]

with $C$ being the maximal moment of inertia of the satellite, and where ${\mathpzc T}_{z}^{\rm(TRI)}$ and ${\mathpzc T}_{z}^{\rm(TIDE)}$ are traditionally called the triaxiality-caused (or simply triaxial) torque and the tidal torque.

There is general agreement in the literature as to the expression for ${\mathpzc T}_{z}^{\rm(TRI)}$ (see for instance Danby [17]): the triaxial torque is written as an infinite Fourier series, which is usually truncated, since only a few resonances are really relevant. By contrast, the expression for ${\mathpzc T}_{z}^{\rm(TIDE)}$ is a much more delicate issue. In the first paper devoted to the problem, by Goldreich and Peale [28], the MacDonald model was mainly used for the tidal torque. This entailed a simple form for ${\mathpzc T}_{z}^{\rm(TIDE)}$, well suited to both analytical and numerical computations. However, the results obtained with such a model were rather disappointing: the probability of capture in the 3:2 resonance was found to be very small (7%), and only assuming a chaotic evolution of the eccentricity of Mercury — as demonstrated much later by Correia and Laskar [15] — can it become much higher (55%). Here, the probability of capture is defined as the probability of the satellite being trapped in a resonance when crossing it; if one takes into account the strong variations which the eccentricity underwent in the past, the probability is highly enhanced because multiple crossings become possible.

Moreover, tidal models such as MacDonald’s, based on a constant time lag (CTL), lead to the existence of a stable quasi-synchronous solution, which in the case of Mercury is characterised by a spin rate $\dot{\theta}\approx 1.26n$, where $n$ is the mean motion of Mercury; such a solution turns out to attract most of trajectories in the case of constant eccentricity. More generally, the CTL model produces rather nonphysical results when applied to satellite-primary systems, as the most common resonance for satellites is the synchronous one (1:1).

Recently, the physical validity of tidal models based on constant time lag was strongly questioned by Efroimsky and Makarov [18, 19, 21, 33]. A more realistic model has been introduced by Efroimsky [19], based on the Darwin-Kaula expansion of the tidal torque [30], which takes into account both the rheology and the self-gravitation of Mercury. By relying on such a model, Makarov showed that the probability of capture in the 3:2 resonance is 100%, that is to say when Mercury crosses the 3:2 resonance it is inevitably entrapped in it [31]. Later on, Noyelles et al. studied the case of non-constant eccentricity and found not only that trapping in the 3:2 resonance is the most probable outcome of the time evolution, but also that the trapping time is much smaller than that predicted by previous theories [38]. We refer to [38] for more details and for a very clear discussion of the existing results in the literature; in the following we shall refer to the system with the tidal model used in [38] as the NFME model. We note that, in the NFME model, the tidal torque is not a smooth function (it is only $C^{1}$) and it appears as a superposition of kinks: this makes both the numerical and the analytical investigations rather subtle.

In general, because of the small value of the dissipation, integrating the equations of motion requires very long times. Thus, it may be convenient for practical purposes to make some assumptions on the initial data of the system: usually one fixes the initial velocity and considers a large sample of initial phases (say, 1000). Therefore, a very high probability of capture (even 100%) in a given resonance does not necessarily imply that every trajectory ends up in that resonance, because, for that to happen, one needs the trajectory not to have been trapped earlier by other resonances that it has crossed during its time evolution. Since the initial condition is not known, it may be important to investigate a larger sample of initial conditions by varying $\dot{\theta}$ as well as $\theta$, randomly distributed in the phase space. A fast numerical integration method was proposed by Bartuccelli et al. [5, 6], which allows consideration of a larger number of initial data (say, 50 000). Then one can evaluate the probability of capture in a given resonance as the fraction of trajectories which are eventually attracted into that resonance: for the NFME model it was found that the 3:2 resonance is still the most probable final state, since it attracts about 42% of the trajectories with initial conditions $(\theta,\dot{\theta})$ inside the set $[0,2\pi]\times[0,5n]$, and even more if one takes into account only initial conditions above the 3:2 resonance (see Section 7).

In this paper, we study more closely, both numerically and analytically, the nature of the attractors for the NFME model — something which is still missing in the literature. Indeed, as noted in [34], the tidal torque in the NFME model “leaves little room for analytical applications” and usually very strong simplifying assumptions on the equation are made in order to obtain analytical expressions for approximate solutions. In fact, when one speaks of a resonance $p$:$q$, usually one simply means that the solution $\theta(t)$ to (1.1) is such that $\dot{\theta}(t)\approx pn/q$, with $n$ being the mean motion, but it is not obvious at all whether the solution is periodic, i.e. has frequency commensurate with $n$.

Numerically, one finds a finite number of periodic attractors. The exact number depends on the truncation of the triaxial torque ${\mathpzc T}_{z}^{\rm(TRI)}$. It is important to note that, if we enlarged the number of harmonics included in ${\mathpzc T}_{z}^{\rm(TRI)}$ with respect to the truncation used in [38], it is true that new periodic attractors would appear, but in fact they would attract only a few trajectories: as a consequence, the general scenario would remain essentially unchanged. Using the same truncation as in [38], the periodic solutions in Figure 1 are found numerically. Not all such solutions are stable (see Section 5 for details); the unstable ones do not correspond to attractors. In particular, the periodic solutions which are unstable include the 3:2 resonance. This may seem a bit surprising, since the 3:2 resonance is expected to be the dominant one.

In fact, an attracting solution with $\dot{\theta}(t)\approx 3n/2$ is found numerically. However, if we look more carefully at such a solution, we realise that it does not appear to be a periodic solution. More precisely, if we write $\theta(t)=3nt/2+z(t)$ and plot $\dot{z}(t)/n$ versus $z(t)$, we obtain the curve in Figure 2 (the function $z(t)$ describes the librations of the spin rate). So, the attracting solution has a much more complicated structure with respect to the periodic solution with $\dot{\theta}(t)\approx 3n/2$ depicted in Figure 1. Apparently, the solution is characterised by two frequencies: there is a fast oscillating motion superimposed on a slow oscillation. This is confirmed by a Fourier Transform analysis: the dynamics involves two frequencies $n$ and $\omega$, with $n\approx 73.9\omega$ (see Section 4). We term such a solution a quasi-periodic attractor (see Remark 6 in Section 4, though).

So, from a numerical point of view, we find that the main attractor of the system (1.1) does not really correspond to what one usually means by a resonance, that is, a periodic solution with frequency commensurate with the forcing frequency. As a matter of fact, while the mechanism by which periodic attractors appear in a periodically perturbed system is rather clear, as follows from Melnikov’s theory [29], the appearance of quasi-periodic solutions is much less standard and deserves further investigation.

We aim to provide an analytical description of the attractors represented in both Figures 1 and 2. Essentially, we shall use perturbation theory, but with some caveats, since the tidal torque is not a smooth function and has very rapid variations. We shall see that a few steps of perturbation theory are sufficient to provide an analytical expression for the solutions which is in very good agreement with the numerical results. However, before entering into the mathematical details, let us discuss briefly — and informally — what kind of solutions may be expected.

We can rewrite (1.1) as

[TABLE]

where $\varepsilon$ and $\gamma$ are positive numbers; the explicit expression of the functions $F$ and $G$ will be given in Section 2. The number $\varepsilon$ is small, so it plays the role of a perturbation parameter. Let us consider first the case in which $\varepsilon\neq 0$ and $\gamma=0$ (conservative limit). In that case, for $\varepsilon$ small enough, one has a quasi-integrable system and hence KAM theorem applies [2]: most invariant tori persist, while the resonant ones are destroyed. Therefore, most of the solutions are quasi-periodic. This does not mean that periodic solutions do not exist: what happens is that of each entire resonant torus only a finite number of trajectories survive in the presence of the perturbation. When the dissipative term is also present ($\gamma\neq 0$), the scenario changes drastically: while the periodic solutions persist, the quasi-periodic solutions disappear almost completely. Moreover, the periodic solutions assume a pivotal role, since they become attractors. In practice, both $\varepsilon$ and $\gamma$ are different from zero, so that the scenario may be somewhat different. In particular, according to the exact form of the dissipative term, quasi-periodic attractors are still possible for the full system. For instance, this is what happens when the MacDonald torque is considered [12, 4]: the quasi-periodic solution corresponds to the quasi-synchronous solution which is also found numerically. While no other quasi-periodic solutions are observed, analytically the existence of the quasi-synchronous solution is rather tricky to prove. Indeed, in order to implement a KAM-like scheme, one needs to assume that the frequencies of the solution are strongly non-resonant — in practice Diophantine (see for instance [12, 5, 36]). However, the frequencies depend continuously on the parameters and the non-resonance condition is not necessarily satisfied when the parameters are varied.

In the case of the realistic tidal model considered by Efroimsky and Makarov, stable quasi-synchronous solutions are not possible [33]. The attracting solutions are not necessarily periodic though. Indeed, for $\varepsilon$ small enough, periodic solutions corresponding to resonances are expected to exist and be attractive, as follows from Melnikov’s theory, but, when increasing the value of $\varepsilon$, bifurcation phenomena may occur: a periodic solution may become unstable and a new kind of solution may appear nearby. For the corresponding Poincaré map — or a suitable iteration of it — the periodic solution gives a fixed point; then the new solution looks like a curve, which emerges from the fixed point by Hopf bifurcation when $\varepsilon$ crosses a threshold value [35].

A natural question is whether one can account analytically for the attracting solutions that we have just described. As far as the solution is periodic, one can apply Melnikov’s theory for subharmonic solutions; see for instance [25] and references therein. Care is in order in that case, because the function $F$ in (1.2) is only $C^{1}$ and has rapid variations, so one cannot apply perturbation theory as usually implemented in the case of smooth functions.

A mathematical description of the solution appearing by Hopf bifurcation is more demanding. To understand how to proceed, let us make a step back and consider the conservative system with $\gamma=0$. In that case, as we have said, the resonant tori are destroyed. In addition, the closest tori also undergo the same fate, so that narrow gaps appear in phase space where the resonant tori have disappeared. It is in these gaps that the periodic orbits are located. Nevertheless, not all the motions inside the gaps are chaotic. In fact, there are a lot of Lagrangian tori [2, 36, 11]: if one considers the Poincaré map, the tori look like closed curves encircling the fixed point which corresponds to the periodic solution. Of course, such tori are not KAM tori (that is, they are not deformations of the tori of the unperturbed system), but they still correspond to quasi-periodic solutions.

Now, let us consider what happens when the dissipation is taken into account. If the dissipation is large enough, then all the tori disappear and the fixed point turns out to be asymptotically stable: such a scenario corresponds to the periodic solutions found by Melnikov’s theory. On the other hand, if the dissipation is small, the force may prevail and push away from the fixed point (which still exists, but becomes unstable, this being an effect of the force dominating the dissipation) and one of the curves nearby survives and is stable. Such a curve is the closed curve which arises by Hopf bifurcation from the fixed point.

To summarise, we find an attracting solution close to each resonance: whether such a solution is periodic or has a more complicated structure depends on the values of the parameters. For fixed values of the parameters, it may happen that some periodic solutions are stable (and hence are attractors) and others are unstable (and so a quasi-periodic solution appears nearby); in the case of Mercury-Sun, for the physical values of the parameters — see Section 2 — we find that the attracting solution corresponding to the 3:2 resonance is quasi-periodic. Both kinds of solutions can be found by perturbation theory. The periodic solutions to the full system (1.2) are obtained by starting from the periodic solutions of the unperturbed system ($\varepsilon=0$). By contrast, the quasi-periodic solutions are not perturbations of solutions of the unperturbed system: one has to perform first a change of variables (a ‘normal form’, in the language of KAM theory), which allows us to write the system as a perturbation of a new unperturbed system which is essentially a pendulum. Then the unperturbed solution to be continued is an oscillatory solution of the pendulum.

In the forthcoming sections, we use perturbation theory to compute approximations to the periodic and quasi-periodic solutions described above. To overcome the lack of smoothness of the vector field, we use an iteration method based on the same idea as the Picard approximants. In principle, one should prove that the iteration scheme converges, but we do not address this issue here. We confine ourselves to computing a few steps of the iteration and show that the results are in very good agreement with the numerical simulations. The problem of convergence is certainly non-trivial, particularly in the case of quasi-periodic solutions. In fact, as pointed out previously, in KAM-like problems one usually assumes strong resonance conditions on the frequencies, whereas, in the case of dissipative systems, the frequencies are expected to depend continuously on the parameters. This means that, for some values of the parameters the motions can be periodic, and still invariant tori exist. We mention that existence of quasi-periodic solutions in both conservative and dissipative systems, without assuming a non-resonance condition on the frequencies, was proved in different contexts [8, 9, 10, 26], but the results do not apply to the class of systems we are considering.

Of course, without discussing the issue of convergence of the iteration scheme, we cannot conclude that the solution we look for really exists. Indeed, the solution is defined as the limit — if it exists — of the successive approximations found along the iteration. When computing the quasi-periodic solution corresponding to the 3:2 resonance, we can claim only a posteriori, after comparing with the results of the numerical investigation, that such a solution really exists. The main advantages of the analytical approach, with respect the numerical results, are that:

we obtain an explicit formula that approximates reasonably well the quasi-periodic attracting solution; 2. 2.

we provide some insight into the nature of such a solution and on the mechanism of its creation.

The rest of the paper is organised as follows. In Section 2 we define the model, and give the explicit expressions of the functions $F$ and $G$ appearing in (1.1). In Section 3 we study analytically the existence of periodic solutions and find explicit expressions for them. These match very well the solutions that are found numerically, represented in Figure 1. In Section 4 we study the quasi-periodic attractor close to the resonance 3:2 (the most interesting one, since it corresponds to the dominant attractor), by applying perturbation theory after a preliminary step which redefines the unperturbed system: once more we find good agreement with the numerical results as given in Figure 2. In Section 5 we discuss the stability of the periodic solutions and the appearance of the quasi-periodic solutions by Hopf bifurcation, by presenting numerical results which provide further support to the analytical ones. In Section 6 we study analytically the dynamics far from the resonances; we show that several simplifications can be made, which we justify a posteriori by comparison with the numerical results: the approximate analytical solutions fit the numerical ones closely and reveal a slow, almost linear decreasing of the velocity. In Section 7, inspired by the results of the previous sections, we revisit the problem of determining the probability of capture of the existing attractors; we find that, assuming an originally faster rotating Mercury, the probability of capture in the $3$:$2$ resonance is higher than $50\%$. Finally, Section 8 is devoted to the conclusions and a discussion of open problems.

2 The spin-orbit model with a realistic tidal torque

The spin-orbit model describes an asymmetric ellipsoidal celestial body which moves in a Keplerian elliptic orbit around a central body and rotates around an axis orthogonal to the orbit plane [17, 28, 37].

The ordinary differential equation governing the dynamics of the system is

[TABLE]

where $\theta\in\mathds{T}=\mathds{R}/2\pi\mathds{Z}$ denotes the angle between the longest axis of the body and the line of apsides, the parameters $\zeta,\eta\in\mathds{R}_{+}$ are small and the dots denote derivatives with respect to the time $t$. On the right hand side of (2.1), the term with $G(\theta,t)$, that we call the force in what follows, represents the triaxial torque acting on the system, while the term with $F(\dot{\theta})$ models the dissipation due to the tidal torque. In the literature, $\ddot{\theta}^{\rm(TRI)}:=-\zeta\,G(\theta,t)$ and $\ddot{\theta}^{\rm(TIDE)}:=-\eta F(\dot{\theta})$ sometimes are referred to, slightly improperly, as the triaxial acceleration and the tidal acceleration, respectively. In this paper we focus on the Mercury-Sun system, even though the ideas could be applied to any system formed by a satellite orbiting its primary — of course the tidal model to use strongly depends on the system one is interested in, as forcefully pointed out in the recent literature [20, 39, 18, 19, 32, 21, 33, 22, 23, 24, 34].

The function $G(\theta,t)$ has the form [17, 28, 37]

[TABLE]

where $n$ is the forcing frequency and the coefficients $A_{k}$, which depend on the eccentricity, are related to the Hansen coefficients [37] $G_{20q}$ by $A_{k}=G_{20q}$, with $k=q+2$. For the Mercury-Sun system, for which $n=26.0879\,\hbox{ yr}^{-1}$, the coefficients $A_{k}$ for $k\neq 0$ are given in Table 1, with $A_{0}=0$.

The function $F(\dot{\theta})$ is given by [20, 19, 31, 32, 38]

[TABLE]

where

[TABLE]

with the functions $\mathcal{I}(\omega)$ and $\mathcal{R}(\omega)$ given by

[TABLE]

where $\Gamma$ is the Gamma function. The values of the constants appearing in (2.4) and (2.5) are $\alpha=0.2$, $\tau_{M}=\tau_{A}=500$ and $\mathcal{A}=38\pi\mu R^{4}/3GM^{2}=15.51726$, where $\mu=7.967\times 10^{28}$ km*-1yr-2*, $R=2.44\times 10^{3}$ km, $M=3.301\times 10^{23}$ kg and $G=6.646\times 10^{-5}$ kg*-1km3yr-2* are the unrelaxed rigidity, the radius of Mercury, the mass of Mercury and the gravitational constant, respectively. The parameters $\tau_{M}$ and $\tau_{A}$ are the Maxwell and Andrade times. The tidal acceleration $\ddot{\theta}^{\rm(TIDE)}=-\eta\,F(\dot{\theta})$ is plotted, on two different scales, in [6, Figures 3 and 4]: there are five pronounced kinks where the function changes sign, the three most significant being at $\dot{\theta}/n\approx 1,3/2,2$.

Finally, the parameters $\zeta$ and $\eta$ in (2.1) are

[TABLE]

where $A<B<C$ are the moments of inertia with respect to the $x$, $y$ and $z$ axes of Mercury. We refer to [38, 6] and the references therein for further details of the spin-orbit model. All the values of the parameters are taken from [38].

3 Method of successive approximations for the periodic attractors

We set $\zeta=\varepsilon$ and $\eta=\gamma\varepsilon$ in (2.1), with $\gamma=0.3243$, and write (2.1) as

[TABLE]

where there is only one parameter. For $\varepsilon=0$ the equation is trivially solved.

In this section we study the existence of periodic solutions to (3.1), by treating $\varepsilon$ as a perturbation parameter. Since the tidal torque is only $C^{1}$, instead of the usual perturbation theory based on power series expansions in $\varepsilon$, we shall rely on a Picard-like iteration method [13], more suited for differential equations with $C^{1}$ vector fields, to find successive approximations to the periodic solutions.

3.1 Zeroth approximation

We look for a periodic solution which continues the unperturbed one with frequency $\omega_{0}=p/q\in\mathds{Q}$, that is, a solution which reduces to $\theta_{0}(t):=\bar{\theta}_{0}+n\omega_{0}t$ as $\varepsilon\to 0$. Let us write (3.1) as an integral equation:

[TABLE]

where the constants $\bar{\theta}$ and $\bar{y}$ have to be fixed by requiring the solution $\theta(t)$ to be periodic with period $2\pi q/n$.

Let us consider as the zeroth approximation the solution $\theta_{0}(t)$ of the equation obtained by setting $\varepsilon=0$, $\theta_{0}(t)=\bar{\theta}_{0}+n\omega_{0}t+\bar{y}_{0}t$, which yields $\bar{y}_{0}=0$, so that

[TABLE]

with $\bar{\theta}_{0}$ being arbitrary. This is the unperturbed solution. To fix the value of $\bar{\theta}_{0}$ we have to consider the equation for the first approximation $\theta_{1}(t)$, which is the solution to the integral equation

[TABLE]

For any periodic function $f$ with period $2\pi q/n$ we denote its average by

[TABLE]

By defining the functions $u_{1}(t)$, $U_{1}(t)$ and ${\mathcal{U}}_{1}(t)$ as

[TABLE]

one needs $\langle u_{1}\rangle=0$ for the function

[TABLE]

to be periodic (and not to grow linearly with $\tau$). This leads to the requirement that

[TABLE]

which, for $G$ as in (2.2), yields $2p/q\in{\mathcal{K}}$ and

[TABLE]

Note that (3.5) can be satisfied if and only if (i) $q=2$ and $p\in{\mathcal{K}}$ and (ii) $\gamma|F(n\omega_{0})/A_{k_{0}}|<1$. If these conditions are fulfilled, the constant $\bar{\theta}_{0}$ is fixed to one of the $4$ values

[TABLE]

which, at least for small values of $\varepsilon$, correspond to two stable and two unstable solutions (attractors and repellers, respectively); this is a consequence of the Poincaré-Birkhoff theorem [2, 14]. We note at this point that bifurcation phenomena may occur when increasing the value of $\varepsilon$ (see Section 5 below).

Therefore the zeroth approximation is given by (3.2), with $\bar{\theta}_{0}$ given by one of the four values in (3.6).

3.2 First approximation

Now we want to compute the first approximation $\theta_{1}(t)$. Once $\bar{\theta}_{0}$ has been fixed in such a way that (3.4) is satisfied, we can compute $U_{1}$ in (3.3). If we require $\langle U_{1}\rangle$ to vanish (for the function ${\mathcal{U}}_{1}(t)$ to be bounded and hence periodic), we obtain, for $G$ as in (2.2),

[TABLE]

where ${\mathcal{K}}_{0}:={\mathcal{K}}\setminus\{k_{0}\}$. Then, by fixing $\bar{y}_{1}=-U_{1}(0)$, and choosing ${\mathcal{U}}_{1}(t)$ so that its average $\langle{\mathcal{U}}_{1}\rangle$ vanishes, one finds $\theta_{1}(t)=\bar{\Theta}_{1}+n\omega_{0}t-\varepsilon{\mathcal{U}}_{1}(t)$, where $\bar{\Theta}_{1}:=\bar{\theta}_{1}-\varepsilon{\mathcal{U}}_{1}(0)$ and, for $G$ as in (2.2),

[TABLE]

Thus, we obtain

[TABLE]

with $\bar{\Theta}_{1}$ to be determined by requiring that the second order $\theta_{2}(t)$ has period $2\pi q/n$ as well.

The second approximation $\theta_{2}(t)$ is the solution to the equation

[TABLE]

where $\theta_{1}(t)$ is the solution (3.7) found at the first iterative step. Therefore, for the solution to be periodic, we need

[TABLE]

where $\xi_{1}(t):=-\varepsilon{\mathcal{U}}_{1}(t)$. If we rewrite (3.9) as

[TABLE]

we see we have to invert the function $\mathcal{G}$ to find $\bar{\Theta}_{1}$. By expanding and using the fact that $\xi_{1}=O(\varepsilon)$, we find

[TABLE]

so that we can write

[TABLE]

By neglecting the corrections of order $\varepsilon$ in (3.11) and defining

[TABLE]

we find

[TABLE]

where again one must have $q=2$ and $p\in{\mathcal{K}}$, so that

[TABLE]

The values of the integrals $J(n\omega_{0})$ in (3.12) are given in Table 2, where they are also compared to the values $\gamma F(n\omega_{0})$. It is evident that there are appreciable discrepancies between the two values for $\omega=1$, 3/2 and 2, which correspond to the major kinks of the function $F$.

The values of the constants $\bar{\Theta}_{1}$ in the interval $(-\pi/4,\pi/4)$ are given in Table 3, where the corresponding values $\bar{\theta}_{0}$ are also given. Once more the difference between the two values $\bar{\theta}_{0}$ and $\bar{\Theta}_{1}$ is larger for resonances such as $\omega_{0}=1$, $\omega_{0}=3/2$ and $\omega_{0}=2$, where the kinks of the function $F$ are more pronounced.

Remark 1**.**

One may wonder why in (3.10) we Taylor-expand the function $G$ but not the function $F$. This is due to the fact that $G(\theta,t)$ is a smooth function (in fact it is analytic), while $F$ is only $C^{1}$. Moreover, the first derivative of $F$ is very large at some resonances such as $\omega_{0}=1$ and $\omega_{0}=3/2$; an explicit computation gives $\eta\partial_{\dot{\theta}}F(n)=24.8421$ and $\eta\partial_{\dot{\theta}}F(3n/2)=13.2493$. Thus, since the size of $\xi_{1}(t)$ is of order $\varepsilon$, for $\varepsilon=\zeta$ fixed as in (2.6), the function $G(\bar{\Theta}_{1}+n\omega_{0}t+\xi_{1}(t),t)$ is well approximated by $G(\bar{\Theta}_{1}+n\omega_{0}t,t)$, while $F(n\omega_{0}+\dot{\xi}_{1}(\tau^{\prime}))$ can be appreciably different from $F(n\omega_{0})$.

---

Remark 2**.**

The reason why we neglect the corrections of order $\varepsilon$ in (3.11) is that, then, the implicit function equation (3.10) is easily solved. Of course, in so doing, an error is introduced. However, we argue that such an error is of the same order of magnitude of the terms we are disregarding by stopping Picard’s iteration at the first step. Indeed, as the difference between $\theta_{1}(t)$ and $\theta_{0}(t)$ is of size $\varepsilon$, so the difference between $\theta_{2}(t)$ and $\theta_{1}(t)$ is expected to be of size $\varepsilon^{2}$. Taking into account the terms $O(\varepsilon)$ in (3.11) would lead to values of $\bar{\Theta}_{1}$ which differ from those given by (3.13) by terms of size $\varepsilon^{2}$ — the same size as the corrections to be found in the next iterative step.

---

In conclusion, the first approximation $\theta_{1}(t)$ is given by (3.7), with the constants $\bar{\theta}_{0}$ and $\bar{\Theta}_{1}$ taken from Table 3. By defining the libration $z(t):=\theta_{1}(t)-n\omega_{0}t$, one finds that, in the $(z,\dot{\theta})$-plane, the periodic attractors are described by the curves in Figure 3. There is an attractor for any value of $\omega_{0}=p/q$, with $q=2$ and $p\in{\mathcal{K}}\setminus\{0\}$. Each attractor is characterised by the property that it describes an oscillation $\theta(t)$ with $\dot{\theta}(t)$ close to $n\omega_{0}$, while $\theta(t)-n\omega_{0}t$ moves around $\bar{\Theta}_{1}$. Note that the amplitude of the libration is essentially determined by the triaxial torque, as expected for relatively cold celestial bodies with high viscosity and large Maxwell time $\tau_{M}$ [34]; however its centre is fixed by the tidal torque, through (3.12).

Remark 3**.**

The first approximation is not a first order perturbation theory solution. Indeed, if one considered the Taylor expansion in $\varepsilon$ of the solution and kept only the terms up to the first order, then one should replace $J(n\omega)$ with $\gamma\,F(n\omega)$ in (3.10), so that $\bar{\Theta}_{1}$ would become a correction of order $\varepsilon$ of $\bar{\theta}_{0}$. Thus, the solution would be still of the form (3.7), but with a different value for $\bar{\Theta}_{1}$. An explicit computation shows that the value inside the interval $(-\pi/4,\pi/4)$ would be appreciably smaller than $\bar{\Theta}_{1}$, as given in Table 3, so that the curves would appear slightly shifted in the $z$-direction with respect to those represented in Figure 3. This effect would be noticeable in all periodic attractors.

---

3.3 Comparison with the numerical results

For all the numerical computations — both in this section and in the following ones — we rely on the fast numerical integrator introduced in [5, 6], to which we refer for details. In order to compute the attractors, we follow the time evolution of a large number of initial conditions in phase space (more than 50 000) so as to determine the asymptotic behaviour of the corresponding trajectories — see also [3]. By using a large number of initial conditions we argue that even if, besides those that have been detected, other attractors existed, they would be irrelevant, since they would attract only a negligible fraction of trajectories. Of course, if other periodic attractors exist, they are not obtained by continuation from the unperturbed solutions of the form (3.2). The periodic solutions which are found numerically are represented in Figure 1. Not all of them are attractors (see Section 5 below): the periodic solutions which are unstable may be detected by a fixed point method by considering the corresponding Poincaré section.

A comparison with the analytical results shows that there is a very good agreement so that we conclude that the first order approximation $\theta_{1}(t)$ provides an accurate description of the solution.

Actually, a shift is observed in the case of the resonances $\omega_{0}=-1/2$ and, particularly, $\omega_{0}=-1$. The reason behind that is very likely due to the smallness of the coefficients $A_{k}$ corresponding to the retrograde resonances ($\omega_{0}=-1/2$ and $\omega_{0}=-1$), which makes the dependence of $\bar{\Theta}_{1}$ very sensitive to the exact value of the integral $J(n\omega_{0})$. For instance an error in the fourth decimal digit, which has no effect for the prograde resonances, is able to produce an appreciable shift for the retrograde resonances.

4 The quasi-periodic attractor close to the 3:2 resonance

We study now the existence of quasi-periodic solutions to (2.1) describing invariant tori. Consider the ordinary differential equation (2.1), with the functions $G(\theta,t)$ and $F(\dot{\theta})$ as in (2.2) and (2.3), respectively. We look for a solution $\theta(t)$ with $\dot{\theta}(t)\approx 3n/2$, which suggests setting

[TABLE]

In terms of $\xi$ the equation of motion becomes

[TABLE]

which we rewrite as

[TABLE]

where $\omega:=\sqrt{2\zeta A_{3}}=0.3534$ and ${\mathcal{K}}_{0}={\mathcal{K}}\setminus\{k_{0}\}$, with $k_{0}=3$.

If in (4.2) we set $\zeta=\varepsilon$ and $\eta=\gamma\varepsilon$, without affecting $\omega$, (4.2) becomes

[TABLE]

where we have defined $B_{i}:=A_{i+3}=G_{20i+1}$ and ${\mathcal{I}}_{0}=\{\pm 1,\pm 2,\pm 3,\pm 4,\pm 5\}$. We would like to study (4.3) by considering $\omega$ as a parameter independent of $\varepsilon$. Thus, (4.3) could be seen as a perturbation of the pendulum equation, $\ddot{\xi}=-\omega^{2}\sin\xi$, to which it reduces as $\varepsilon=0$. However, since we look for a solution to (4.3) with $\xi$ close to [math], to simplify the analysis we proceed in a slightly different way. Starting from (4.2), we make a linear approximation of the unperturbed system by splitting $\sin\xi$ into two terms and re-writing (4.2) as

[TABLE]

where again $\omega$ is to be considered fixed, while $\varepsilon$ can be varied. Of course, eventually, we have to fix $\varepsilon$ to the value such that $2\varepsilon A_{3}=\omega^{2}$. Note that a similar approach has been proposed in [40, 27] to study the secondary resonances in non-dissipative spin-orbit models with $\varepsilon$ close to a rational number and small values of eccentricity — a situation very far from that of Mercury.

In general, quasi-periodic solutions to (4.4) will have two frequencies $\omega_{L}$ and $n$, with $n$ fixed and $\omega_{L}$ close to $\omega$ — see [3, 25, 41]. We call $\omega_{L}$ the low frequency and $n$ the high frequency, because $\omega_{L}\approx\omega$ and $n=73.82\omega$. Hence we write the solution as $\xi(t)=X(\boldsymbol{\omega}_{L}t)$, where $\boldsymbol{\omega}_{L}:=(\omega_{L},n)$. To take into account the possible dependence of frequency on the perturbation we write (4.4) as

[TABLE]

where $\omega^{2}=\omega_{L}^{2}+\mu\,\varepsilon$, with $\omega_{L}$ and $\mu$ constants to be determined. As a further simplification, we approximate $\sin\xi\approx\xi-\xi^{3}/6$, so that (4.5) becomes

[TABLE]

We look for approximate solutions to (4.6) using a Picard-like iteration scheme as in Section 3. The strategy is the following. The zeroth approximation is the solution to (4.6) with $\varepsilon=0$. Suppose that at step $k\geq 0$ we have found an approximation $\xi_{k}(t)=X_{k}(\boldsymbol{\omega}_{L}t)$, with frequency vector $\boldsymbol{\omega}_{L}=(\omega_{L},n)$, for some $\omega_{L}$ close to $\omega$. As we shall see such a solution depends on an arbitrary parameter $C_{k}$. Then we compute the $(k+1)-$th approximation as the solution of

[TABLE]

where we consider $\omega_{L}$ as a parameter related to $\mu$ through the relation $\omega^{2}=\omega_{L}^{2}+\mu\varepsilon$, with $\omega$ given. For such a solution to be bounded and hence quasi-periodic with frequency vector $\boldsymbol{\omega}_{L}$, we have to fix the parameters $C_{k}$, $\mu$ and $\omega_{L}$, the latter two so as to also satisfy the constraint that $\omega^{2}=\omega_{L}^{2}+\mu\varepsilon$; the value of the low frequency $\omega_{L}$ will be determined by an implicit function problem. The corresponding solution will be a quasi-periodic function $\xi_{k+1}(t)=X_{k+1}(\boldsymbol{\omega}_{L}t)$, with $\boldsymbol{\omega}_{L}$ slightly different from the value found at the previous step, depending on a new arbitrary constant $C_{k+1}$ to be fixed at the next iteration step to a value close to $C_{k}$.

4.1 Zeroth approximation

The zeroth approximation $\xi_{0}(t)$ is the solution to the equation

[TABLE]

obtained from (4.6) by setting $\varepsilon=0$. The solution to (4.7) is $\xi_{0}(t)=\bar{\xi}_{0}\cos\omega_{L}t+\bar{y}_{0}\sin\omega_{L}t$, where $\bar{\xi}_{0}$ and $\bar{y}_{0}$ are related to the initial data through the relation $\bar{\xi}_{0}=\xi(0)$ and $\bar{y}_{0}=\dot{\xi}(0)/\omega$.

For any initial datum, we can write

[TABLE]

where $C_{0}>0$ and $\varphi_{0}$ is the initial phase. By varying the initial phase $\varphi_{0}$ the corresponding trajectories in the space $(\xi,y,t)$ describe an invariant torus which appears as a right circular cylinder with axis along the $t$-axis and radius $C_{0}$.

For $\varepsilon\neq 0$ we expect the cylinder to persist, albeit deformed with respect to the unperturbed case. This allows us to fix arbitrarily the phase $\varphi_{0}$, which plays no role: we can put it equal to [math] for convenience. Then the torus is parameterised in terms of the value $C_{0}\in\mathds{R}_{+}$ at which it crosses the positive $\xi$-axis (corresponding to the phase $\varphi_{0}=0$). We expect that the values of $C_{0}$ and $\omega_{L}$, which are arbitrary for $\varepsilon=0$, will be fixed when we set $\varepsilon\neq 0$ and take into account the dissipation by the requirement that the perturbed solution is still bounded.

Remark 4**.**

Since we are making a linear approximation, all unperturbed solutions have the same frequency $\omega_{L}$. Therefore all the approximate quasi-periodic solutions (4.8), by construction, will have the same frequency vector $(\omega_{L},n)$. As we shall see, if we wish to compute the higher order approximations, we need to take into account the change of frequency $\omega_{L}$ with respect to the linearised unperturbed system — see [7, 5]. In this regard, we note that, even though usually perturbation theory computations are easier in terms of action-angle variables, in our case it is more convenient to work with Cartesian coordinates, because we are looking for a solution around the origin, where the action-angle variables are singular — see [16] and references therein for similar comments.

---

4.2 First approximation

The first order approximation is obtained as the bounded solution to the equation

[TABLE]

with $\xi_{0}(t)$ given by (4.8). The nonhomogeneous linear equation (4.9) can be written as a first order differential equation in $\mathds{R}^{2}$,

[TABLE]

where $\xi_{0}(t)=C_{0}\sin\omega_{L}t$ and $\dot{\xi}_{0}(t)=C_{0}\omega_{L}\cos\omega_{L}t$. More generally, the approximation at step $k$ is defined as the solution to

[TABLE]

where $\xi_{k-1}(t)$ is the approximation found at step $k-1$. At any step $\omega_{L}$ is to be considered a free parameter, close to $\omega$, to be fixed by requiring the solution to be bounded.

The general solution to (4.11) is

[TABLE]

where $(\bar{\xi}_{k},\bar{y}_{k})$ is the initial condition. This leads to

[TABLE]

For $k=1$ we write $\bar{\xi}_{1}\cos\omega_{L}t+\bar{y}_{1}\sin\omega_{L}t=C_{1}\,\sin(\omega_{L}t+\varphi_{1})$ and set $\varphi_{1}=0$ as done for $k=0$. If we define

[TABLE]

we can write the solution $\xi_{1}(t)$ as

[TABLE]

For the functions in (4.14) to be bounded one needs

[TABLE]

The two limits can be computed as

[TABLE]

since the time average of a quasi-periodic function equals the average over the torus. After expanding

[TABLE]

and using

[TABLE]

by parity, and

[TABLE]

and

[TABLE]

we obtain

[TABLE]

The two equations (4.16) are coupled, since, even though $\mu$ does not appear explicitly in (4.16a), the frequency $\omega_{L}$ must satisfy the relation $\omega^{2}=\omega_{L}^{2}+\mu\,\epsilon$. If $C_{0}\neq 0$, (4.16b) gives $\mu=A_{3}C_{0}^{2}/4$, while $C_{0}=0$ would not fix $\mu$ to any value. We fix $C_{0}$ in such a way to make the integral

[TABLE]

vanish. One has

[TABLE]

which shows that $C_{0}=0$ is a zero of $I_{1}(C_{0})$. To study the existence of other zeroes, we compute numerically the integral $I_{1}(C_{0})$, by approximating $\omega_{L}$ with $\omega$. The result, given in Figure 4, shows that there are no other zeroes.

This fixes the value of $C_{0}=0$ in the zeroth order approximation (4.8). With $C_{0}=0$, the value of $\mu$ in (4.16b) is left undetermined and can be set equal to [math]. With the two values $C_{0}$ and $\mu$ being fixed, the right hand side of (4.10) becomes determined and $\Xi_{1}(t)$ in (4.14) can be computed explicitly. One has

[TABLE]

so that, by computing the integrals explicitly, one finds

[TABLE]

where

[TABLE]

and $\omega_{L}$ has to be fixed at the value $\omega$. Since one has $\rho=5.242\times 10^{-5}$ and $2\varepsilon B/\omega=-8.819\times 10^{-3}$, for $\omega_{L}=\omega$, the last contribution in (4.18) is negligible, so that we can approximate

[TABLE]

In (4.20) the constant $C_{1}$ is still arbitrary. In order to obtain the full expression for the first approximation, we have to study the equation for the second approximation and impose the requirement that the solution be bounded.

Of course, if we are interested in the first order contribution, we have to set $\omega_{L}=\omega$. However, as explained before, if we want to compute the second order approximation we have to leave $\omega_{L}$ as a free parameter, to be fixed together with $C_{1}$ by requiring the second order approximation to remain bounded.

4.3 Second approximation

The second approximation is obtained as the solution to the equation

[TABLE]

where the function $\xi_{1}(t)$ will be approximated by (4.20) in the following. Therefore $\xi_{2}(t)$ is given by (4.2) for $k=2$, where we write once more

[TABLE]

with $\varphi_{2}$ fixed to be zero and $C_{2}>0$ to be determined by imposing that no secular terms appear when computing the third approximation. Therefore, we can write the second approximation as

[TABLE]

where $M_{2,c}(t)$ and $M_{2,s}(t)$ are defined in (4.13) with $k=2$. For the function (4.21) to be bounded one needs

[TABLE]

We can write $I_{2,c}=I_{2,c,1}+I_{2,c,2}+I_{2,c,3}+I_{2,c,4}$ and $I_{2,s}=I_{2,s,1}+I_{2,s,2}+I_{2,s,3}+I_{2,s,4}$, where

[TABLE]

and, analogously,

[TABLE]

with

[TABLE]

Note that $b(nt)=\dot{\beta}(nt)$, where the derivative is with respect to $t$. In particular $\beta(\psi)$ is odd and hence $b(\psi)$ is even, so that $I_{2,s,2}=0$. Moreover $I_{2,c,3}=I_{2,c,4}=0$ and, since for $\alpha\neq 0$ one has

[TABLE]

also $I_{2,c,1}=0$. As a consequence one has $I_{2,c}=I_{2,c,2}$ and $I_{2,s}=I_{2,s,1}+I_{s,2,3}+I_{2,s,4}$. The integrals $I_{2,s,3}$ and $I_{2,s,4}$ can be easily computed, up to corrections, and give

[TABLE]

where terms of order $\alpha\|\beta\|_{\infty}$ have been neglected in the latter, since $\|\beta\|_{\infty}=1.644\times 10^{-4}$.

In conclusion (4.22) leads to the equations

[TABLE]

The integral $I_{2}(a)$ can be computed numerically, by approximating $\omega_{L}=$ with $\omega$. The outcome is given in Figure 4. The results suggest the existence of three simple zeroes; besides $a=0$, there are two other zeroes $a=\pm 3.174\times 10^{-3}$. Since we are taking $C_{1}>0$ and hence $a>2\varepsilon B=-0.3116\times 10^{-3}$, only the positive zero has to be considered. The corresponding value of $C_{1}$, computed once more at $\omega_{L}=\omega$, is

[TABLE]

which gives $\alpha=8.981\times 10^{-3}$ in (4.20).

Remark 5**.**

In principle one should consider also the zero $a=0$ (which would give $C_{1}=8.817\times 10^{-2}$). However such a value makes $\mu$ disappear from the equations (4.25) and hence does not fix $\omega_{L}$. To study the fate of such a solution, in particular to see whether it does correspond to a solution of the full equation, one should go to higher orders, where we expect such a solution to disappear — see also Section 8. Here we focus on the solution corresponding to $a\neq 0$, since the latter correctly describes the attractor relevant for the dynamics, as shown by the comparison below with the numerics.

---

In (4.25b) we can approximate the integral by expanding

[TABLE]

where the other terms of order equal to or higher than three have been neglected. The only non-zero contribution to the integral is

[TABLE]

where the explicit form of $\beta(\psi_{2})$ has been used. Only the contributions with $j=\pm i$ are non-zero, so that eventually we obtain

[TABLE]

An explicit computation, with $\omega_{L}=\omega$, gives

[TABLE]

Inserting (4.28) into (4.27) yields

[TABLE]

and hence $\mu\varepsilon=1.034\times 10^{-5}$. This fixes $\omega_{L}$ to a value such that $\omega^{2}-\omega_{L}^{2}=\mu\varepsilon$, so that $\omega_{L}-\omega\approx 1.46\times 10^{-5}$; hence $\omega_{L}$ is very close to $\omega$. For the value of $C_{1}$ in (4.26) and $\omega_{L}$ approximated by $\omega$, we consider the solution

[TABLE]

with $\alpha/2=4.490\times 10^{-3}$, and define the libration as $z(t):=\theta(t)-(3/2)nt$. Then $(z(t),\dot{z}(t)/n)$ is as plotted in Figure 5.

4.4 Comparison with the numerical results

The approximate solution found in the previous sections has to be compared with that found by numerical analysis, see Figure 2. Fast Fourier Transform analysis produces the spectrum of the derivative of the numerical solution in Figure 6: the figure to the left, on a larger scale, shows the presence of peaks at multiples of the frequency $n$, while the figure to the right, on a smaller scale, shows that there are peaks as well at multiples of the frequency $\omega_{L}=n/73.9034=0.3530$. Therefore, with respect to the approximate analytical solution, one has $\omega-\omega_{L}=4\times 10^{-4}$ and hence $\omega^{2}-\omega_{L}^{2}=0.3534^{2}-0.3530^{2}=2.8256\times 10^{-4}$, which is larger than the analytical value. However, apart from that, the agreement between analytical approximations and numerical results is very good. Indeed, we expect such a difference to be negligible as far as we are interested in the amplitude of the solution. Nevertheless the phase shift will become appreciable on longer timescales.

It seems likely to us that, in order to obtain the correct value of the frequency $\omega_{L}$, one has to go to higher orders of approximation: because of the rapid variations of the function $\Phi$, the computations depend very sensitively on the approximation $\dot{\xi}_{k-1}$ appearing in its argument in (4.2). More precisely, when studying the third approximation $\xi_{3}=C_{3}\sin\omega_{L}t+\Xi_{3}(t)$, we expect the following to happen: requiring the solution to remain bounded should fix $C_{2}$ to a value very close to $C_{1}$, while the corrections $\mu$ to the slow frequency should be appreciably larger than the value found at the previous step. Then the second approximation $\xi_{2}(t)$ in (4.21) should provide not only the right amplitude, but the right frequency as well — see Section 8 for further comments.

Remark 6**.**

It is important to point out that, even though we are using the term ‘quasi-periodic’ for the attractor that we have studied in this section, in principle the two frequencies $\omega_{L}$ and $n$ could be commensurate. In fact, since we are determining $\omega_{L}$ numerically, it is not possible to exclude the possibility that the ratio $\omega_{L}/n$ might be a rational number. A better definition could be ‘multi-periodic’ solution, since the solution appears as a superposition of two oscillations, each with its own frequency. However, if we let the parameters change slightly, we find the ratio to be irrational most of the time (see Section 5 below). Thus, we can conclude that the attractor is very likely to be genuinely quasi-periodic. Note also that, even if the ratio were rational and hence the solution were periodic, its period would be very large, since in any case the ratio would not be close to any rational $p/q$, with $q$ small.

---

5 Numerical study of the stability of the attractors

The stability of the periodic attractors for the system described by (2.1) can be studied by considering the corresponding Poincaré map and computing the eigenvalues of the linearised system around its fixed points [35, 29].

As shown in Section 3, there are four periodic orbits for each resonance and hence four fixed points for a suitable iteration of the Poincaré map. Among the fixed points, for $\varepsilon$ small enough, two are stable and two are unstable. We can confine ourselves to the fixed points corresponding to the values of the initial phase of the first approximation $\bar{\Theta}_{1}^{(1)}$ and $\bar{\Theta}_{1}^{(2)}$; indeed, by symmetry, the points corresponding to the phases $\bar{\Theta}_{1}^{(3)}$ and $\bar{\Theta}_{1}^{(4)}$ have the same stability of $\bar{\Theta}_{1}^{(1)}$ and $\bar{\Theta}_{1}^{(2)}$, respectively.

For the values of the parameters given in Section 2, we find the eigenvalues given in Table 4. Numerically, one observes that, when the periodic attractors are both unstable, quasi-periodic attractors such as that studied in Section 4, appear. This happens for the resonances 1:2, 3:2 and 2:1, besides the retrograde resonances which, however, do not play a relevant role in the dynamics of the system. The stable quasi-periodic attractors appear by a Hopf bifurcation, occurring when changing the parameters of the system [35].

As in Section 4, in the following we concentrate on the resonance 3:2, which is the most interesting one in the case of Mercury, but a similar analysis could be easily extended to the other cases. However, as a comparison with the results in Table 4 shows, only quasi-periodic attractors with $\dot{\theta}/n\approx 1/2$, $3/2$ and $2$ are found to exist among the prograde ones, so that we could confine ourselves to these three cases.

Writing (2.1) as

[TABLE]

where $S=\lambda=1$; in terms of $\xi$, defined according to (4.1), the equation becomes

[TABLE]

where the set ${\mathcal{I}}_{0}$ and the coefficients $B_{i}$ are defined as after (4.3). We study the transitions in the dynamics of the system described by (5.1) when we vary either the parameter $\lambda$ (at fixed $S=1$) or the parameter $S$ (at fixed $\lambda=1$).

The bifurcation diagram with the parameter $S$ in Figure 7 shows that the periodic solution with frequency $3n/2$ is stable up to the value $S\approx 0.134$, where the solution loses stability and a stable quasi-periodic solution appears. By increasing $S$, the amplitude of the oscillations increases as well. At $S=1$ we have the quasi-periodic solution studied in Section 4, while the periodic solution with frequency $3n/2$ is unstable. The numerical investigation in [6] demonstrates that the solution with velocity $\dot{\theta}$ close to $3n/2$ is the main attractor for the values of the parameters as in Section 2: more than $42\%$ of the initial conditions $(\theta,\dot{\theta})\in[0,2\pi]\times[0,5n]$ are captured by such an attractor.

In principle, the attractor studied in Section 4, which we call quasi-periodic, could be periodic: indeed, if $\omega_{L}/n$ is a rational number $P/Q$, then the trajectory closes after a suitable time $T$ (large, since $Q$ would be large). Hence, we investigate numerically how the slow frequency $\omega_{L}$ changes when varying a parameter of the system. In Figure 8 we plot the non-fixed frequency $\omega_{L}$ of the quasi-periodic solution as a function of the parameter $S$. The apparent continuity of the curve suggests that, up to a zero-measure set of values of $S$, the two frequencies $\omega$ and $n$ are incommensurate, so that the motion is genuinely quasi-periodic.

On the Poincaré section obtained by sampling the system evolution at integer multiples of $T_{0}=2\pi/n$, the periodic solution appears as a finite set of points (3 for the 3:2 resonance). By contrast, the quasi-periodic solution has support on a torus which intersects the section along a closed curve. In terms of the dynamics on the section we have a Hopf bifurcation, with the transition from a fixed point to a closed curve [35]. A plot of the amplitude $A$ of the quasi-periodic solution versus the parameter $S$ is also given in Figure 8: a fit obtained by looking at the leftmost points in the figure gives $A=A_{0}(S-S_{0})^{\kappa}$, with $A_{0}=7.976\cdot 10^{-5}$, $S_{0}=0.134$ and $\kappa=0.5016$; the numerical value of the exponent $\kappa$ is in agreement with the critical exponent $1/2$ typical of the Hopf bifurcation [35].

From a physical point of view, the bifurcation diagram with parameter $S$ in Figure 7 can be interpreted as follows. For $S=0$ there is a fixed point at the origin, which is asymptotically stable because of the presence of dissipation. When $S$ becomes large enough (that is, if the forcing terms not included in the unperturbed system can not be neglected any longer), the fixed point loses stability and an asymptotically stable closed curve appears nearby by Hopf bifurcation.

A similar phenomenology is observed at $S=1$, when varying the parameter $\lambda$; see Figure 7. Actually we need a dissipation much larger than the physical value $\lambda=1$. Again the transition from the periodic to the quasi-periodic solutions is described by a Hopf bifurcation for the dynamics of the corresponding Poincaré section. The bifurcation diagram with the parameter $\lambda$ in Figure 7 may also be expected on physical grounds. If there is no dissipation, the Poincaré map corresponding to the equation (5.2) has both a fixed point near the origin and a large measure of invariant tori encircling such a point [2, 36, 11]. As soon as $\lambda>0$, all but one of the tori are destroyed, this one being attractive; the fixed point exists as well but it is unstable. By taking larger values of $\lambda$, at some point the torus is destroyed, while the fixed point becomes stable and attracts all trajectories starting from initial data nearby.

6 Dynamics far from the attractors

We now look at the pre-capture dynamics, that is, the period during which the satellite is decelerating but before it has been captured. This will enable us to estimate times to capture, and we accomplish this by making approximations that greatly simplify the dissipation term.

Starting from (2.1), and assuming that (i) $\dot{\theta}>n$ and (ii) $\dot{\theta}$ is not close to any kink, we make the approximation $\eta\,F(\dot{\theta})=a-b\dot{\theta}$. This yields $a,b>0$ — see [6, Figure 3]. For $\dot{\theta}<n$, $a$ would be negative, but the argument below, suitably adapted, would still work. In practice, both $a$ and $b$ are small. For instance, expanding $F(\dot{\theta})$ around $\dot{\theta}=1.75n$ gives $a\approx 1.1\times 10^{-5}$ and $b\approx 1.3\times 10^{-7}$.

It is convenient to rescale time by $\tau=bt$, so that $\dot{\theta}=b\theta^{\prime}$, where the prime denotes the derivative with respect to $\tau$. Using the above, we can approximate (2.1) as

[TABLE]

6.1 The splitting argument

We now make an estimate of the rate at which $\dot{\theta}$ decreases with time, over long time scales, typically of the order of $10^{6}T_{0}$, where $T_{0}=2\pi/n$. We start by defining $s(\tau)$ to be the solution of $s^{\prime\prime}=s^{\prime}-a/b^{2}$ with initial conditions $s(0)=s_{0}$ and $s^{\prime}(0)=s_{1}$, which gives

[TABLE]

We write $\theta(\tau)=s(\tau)+f(\tau)$, in which $s(\tau)$ accounts for the slow decay of $\dot{\theta}$ apparent on a large time scale and $f(\tau)$ represents a small amplitude, high frequency correction term, which is visible only on a smaller time scale — see Figure 9. In other words, we split $\theta(\tau)$ into fast and slow components. Substituting this expression for $\theta(\tau)$ in (6.1), and bearing in mind the ODE obeyed by $s(\tau)$, we find that

[TABLE]

where $\phi=2s_{0}$, $\omega_{k}:=2s_{1}-kn/b$ and where we have used the Taylor series to $O\left(\tau^{2}\right)$ for $s(\tau)$.

We now make the following assumption: since $\omega_{k}$ is $O(1/b)$ and is therefore large, we neglect the term $2f(\tau)$ and so obtain

[TABLE]

for an approximation to the ODE that defines $f(\tau)$. Conveniently, this ODE is linear with constant coefficients, and so is straightforward to solve. With initial conditions $f(0)=f_{0}$ and $f^{\prime}(0)=f_{1}$, we find

[TABLE]

We now consider the initial conditions, noting that we are free to choose initial values for $\theta(0)=\theta_{0}=s_{0}+f_{0}$ and $\theta^{\prime}(0)=\theta_{1}=s_{1}+f_{1}$. Once $\theta_{0}$, $\theta_{1}$ are specified, any values of $s_{0}$, $s_{1}$, $f_{0}$ and $f_{1}$ that satisfy these constraints can be chosen.

Recall that in order to derive (6.2), we assumed that $|f(\tau)|\ll 1$. To be consistent with this, we therefore choose $f_{0}=0$, from which we immediately deduce that $s_{0}=\theta_{0}$, so that $\phi=2\theta_{0}$.

We need to take a little more care over the choice of $f_{1}$. Since $\omega_{k}$ is $O(1/b)$, one has $|\omega_{k}|\gg 1$. Hence, in order to guarantee the smallness of $f(\tau)$, we need if possible to choose $f_{1}$ so as to cancel out the largest terms in the sum in (6.3), these being $\omega_{k}e^{\tau}\cos 2\theta_{0}$ and $-\omega_{k}\cos 2\theta_{0}$, both of which are $O(1/b)$. This cancellation can be accomplished by setting

[TABLE]

whereupon the expression for $f(\tau)$ becomes

[TABLE]

Since $\tau=bt$ with $b\sim 10^{-7}$ and the timescale that we consider is $t\sim 10^{6}$, $e^{\tau}$ is $O(1)$.

Finally we consider $s_{1}$. Since $\theta_{1}$ is given and we have already specified $f_{1}$, we immediately find that

[TABLE]

Hence, since $\theta(\tau)=s(\tau)+f(\tau)$, we have

[TABLE]

We now approximate this expression, starting from the fact that $|\omega_{k}|\gg 1$. Additionally, we express everything in terms of $t$ rather than $\tau$, where $\tau=bt$, so that $b\theta_{1}=\dot{\theta}(0)$. Then (6.5) becomes

[TABLE]

where $\Omega_{k}=b\omega_{k}$, so

[TABLE]

Note that this equation defines $\Omega_{k}$ implicitly — this is a consequence of the way the initial conditions have been assigned. For $n$ sufficiently large, the approximation $\Omega_{k}\approx\dot{\theta}(0)-kn$ will be good. We investigate this approximation in practice in the next subsection.

Differentiating (6.6) we find

[TABLE]

where we have neglected a term of order $b$ in the sum.

6.2 The approximations in practice

Several assumptions have been used to derive (6.6) and (6.8), so we now investigate numerically how good these approximations are in practice.

First of all, we compare ‘exact’ numerical solution $\Omega_{k}$ to (6.7) with the approximation $\Omega_{k}\approx\Omega_{k}^{\rm app}:=\dot{\theta}(0)-kn$, for the two values $\dot{\theta}(0)=1.95n$ and $3.25n$. We find

• •

For $\dot{\theta}(0)=1.95n$, $\max_{\theta_{0}\in\mathbb{T},\,k\in{\mathcal{Q}}}|\Omega_{k}-\Omega_{k}^{\rm app}|\approx 0.073$;

• •

For $\dot{\theta}(0)=3.25n$, $\max_{\theta_{0}\in\mathbb{T},\,k\in{\mathcal{Q}}}|\Omega_{k}-\Omega_{k}^{\rm app}|\approx 0.010$.

Note that $\min_{\theta_{0}\in\mathbb{T},\,k\in{\mathcal{Q}}}\left|\Omega_{k}\right|\approx 6.5$, so we are justified in using the approximation (6.7), since it leads to a relative error of no more than 1.1%.

We concentrate first on (6.8) as an estimate of the dynamics leading up to capture; we can also use this to estimate the time to capture in an orbit of a given spin rate. Despite the fact that the triaxial acceleration has been neglected and the tidal acceleration has been replaced by a simple linear approximation, this simplification gives surprisingly good results — see Figure 10, which shows a large timescale comparison, for times of order $10^{7}$, for two different values of the initial spin rate, $\dot{\theta}(0)=1.95n$ and $3.25n$.

Time-to-capture estimates can be made as follows:

If $\dot{\theta}(0)=1.95n$, then we compute that $a=1.423\times 10^{-5}$ and $b=1.894\times 10^{-7}$. Neglecting the oscillatory term in the last square brackets in (6.8) and defining

[TABLE]

the time taken for $\dot{\theta}(t)$ to decay from $\dot{\theta}(0)$ to $\dot{\theta}=1.5n$ is estimated to be $b^{-1}\ln[(R-0.45n)/R]$, which leads to the value $2.08\times 10^{6}$, to be compared with the numerical value (see Figure 10) of $2.14\times 10^{6}$. 2. 2.

If $\dot{\theta}(0)=3.25n$, then $a=6.733\times 10^{-6}$, $b=2.056\times 10^{-8}$ and $R=-242.70$. In this case, the estimated time to reach $\dot{\theta}=1.5n$ is $8.38\times 10^{6}$, to be compared with the numerical result, which is also $8.38\times 10^{6}$.

For behaviour on a small timescale, of order 1, again see Figure 10, in which, on the right, we plot $\dot{\theta}(t)$ over a time interval of width 2. Only the order of magnitude of the approximation is correct, but this is not surprising given that the approximation completely neglects the triaxial acceleration.

It is interesting to note that the triaxial torque appears to make very little difference during deceleration: Figure 10 gives evidence that it only has an important role to play very close to capture; see also [22, 23].

7 Probability of capture revisited

Since the seminal paper by Goldreich and Peale [28], the probability of capture of a satellite in a resonance $p$:$q$ has usually been studied as the probability $P(p/q)$ for the satellite to be trapped in that resonance when its rotation velocity $\dot{\theta}$ approaches the value $pn/q$. However, it may happen that the satellite never comes close to a given resonance, because it has previously been captured in another one. Therefore, it may be useful to redefine the probability of capture of a given resonance as the fraction of initial conditions whose trajectories are attracted by that resonance. Of course the original rotational state of the satellite is not even approximately known, so that one has to fix a suitable region in phase space from which the initial conditions are taken. For the spin-orbit model defined in Section 2, all the attractors turn out to be contained inside the region $[0,2\pi]\times[-1.5n,4.5n]$, so that, by assuming an initially prograde satellite (and exploiting the $\pi$-periodicity of the equations of motion by considering $\theta$ only in the interval $[0,\pi]$), it is reasonable to confine the choice of initial conditions to the region ${\mathcal{Q}}:=[0,\pi]\times[0,4.5n]$. This is essentially what has been done in [6], where the probability of capture in the $3$:$2$ resonance for Mercury has been estimated to be about $42\%$. Indeed, even if the initial velocity $\dot{\theta}$ is much higher, eventually it decreases almost linearly, up to small oscillations — as suggested by the analysis performed in Section 6 — until it enters the region ${\mathcal{Q}}$.

Nevertheless, there is no reason why the basins of attraction should be uniformly distributed in phase space. In fact, Makarov’s results [31], providing an estimated value $P(3/2)=1$, hint that all solutions starting above the $3$:$2$ resonance never reach the attractors contained in the region $\{(\theta,\dot{\theta}):\dot{\theta}<1.5n\}$: apparently there is a barrier. Hence, it may be worthwhile to study in more detail the distribution of the basins of attraction. Since they look like sets of points with no apparent structure, rather than a picture of the basins it is more illuminating to divide the region ${\mathcal{Q}}$ into nine strips ${\mathcal{Q}}_{0},{\mathcal{Q}}_{1},\ldots,{\mathcal{Q}}_{8}$, each of width $0.5n$ in the velocity direction, and compute the fraction of the basins of attraction which falls inside each strip. The results are given in Table 5, where only the prograde resonances have been considered; we never observed a trajectory with positive initial velocity being attracted by a retrograde resonance.

Table 5 confirms the existence of a barrier associated with the $3$:$2$ resonance: the trajectories starting with $\dot{\theta}>1.5n$ cannot reach the attractors below the $3$:$2$ resonance and, vice versa, the trajectories starting with $\dot{\theta}<1.5n$ cannot reach the attractors above. Surprisingly, an analogous barrier exists associated with the resonance $1$:$1$. More precisely, the basins of attraction of the resonances above $3$:$2$ are contained in the region ${\mathcal{Q}}_{a}:=\{(\theta,\dot{\theta}):\dot{\theta}>1.5n\}$ and, similarly, the basins of attraction of the resonances below $1$:$1$ are contained in the region ${\mathcal{Q}}_{b}:=\{(\theta,\dot{\theta}):\dot{\theta}<n\}$. Moreover, all trajectories with initial velocity $\dot{\theta}\in(n,1.5n)$ are attracted either by the $3$:$2$ resonance or by the $1$:$1$ resonance. Obviously, we cannot claim that the barriers completely obstruct the passage of trajectories. In principle, it is possible that holes appear in the barriers; in any case, crossings of the barriers — if they occur at all — appear to be very rare events.

Another interesting point seen in the results of Table 5 is that the basin of attraction of the $3$:$2$ resonance dominates the strip ${\mathcal{Q}}_{3}:=\{(\theta,\dot{\theta}):1.5n<\dot{\theta}\leq 2n\}$, the area of its intersection with the strip being about $99\%$ of the whole strip, and becomes more or less uniformly distributed above the resonance $2$:$1$. Indeed about $51\%$ of the initial conditions in each strip with $\dot{\theta}>2.5n$ end up being captured in the $3$:$2$ resonance. Therefore, if we assume the initial velocity $\dot{\theta}$ to be high enough (as explained above, in practice this means it is sufficient to fix the initial condition above the highest ($4$:$1$) resonance; usually in the literature one takes a value $\dot{\theta}\approx 4.4n$; see for instance [13, 38]), then the probability of capture in the $3$:$2$ resonance is more than $50\%$ — a value higher than that given in [38], and comparable with that found by Correia and Laskar [13] for the CTL model with varying eccentricity.

8 Conclusions

We have studied the attractors of the spin-orbit model with the realistic tidal torque used in [38]. First, we have investigated numerically the dynamics of the system: besides the periodic attractors for which the frequency is locked in a resonance with Mercury’s mean motion, quasi-periodic attractors are also detected. Which attractors arise actually depends on the values of the parameters; in particular quasi-periodic attractors bifurcate from periodic solutions when the latter become unstable. For the physical values of the parameters, the main attractor of the Mercury-Sun system is a quasi-periodic attractor close to the resonance 3:2. This means that, according to the NMFE model, the librations of the spin rate are quasi-periodic in time in the case of Mercury.

Thereafter, by using a suitable iteration scheme based on perturbation theory for non-smooth systems, we have provided an explicit analytical expression for the attracting solutions: such expressions, despite being obtained after a few steps of the iteration and hence being only approximate, match closely the numerical solutions. So we deduce a posteriori that the perturbative approach provides a reliable description of the dynamics.

However, there are a few aspects which the analysis we have performed does not account for:

The quasi-periodic attractor close to the 3:2 resonance has two frequencies: the fast one is the mean motion $n$, while the slow one depends on the parameters. We expect the latter to be slightly different from the frequency $\omega$ obtained by taking into account just a single harmonic from the Fourier expansion of the triaxial torque — see (4.2). The second approximation computed in Section 4 is correct in describing such a phenomenon, but, from a quantitative point of view, does not provide the right value — that is to say, the value found numerically (see Section 4.4). 2. 2.

When constructing the second approximation, we have to impose the condition that a certain integral vanishes — see (4.25a) — in order to fix the amplitude $C_{1}$ of the leading term. This leaves two values: one of them corresponds to the solution which correctly describes the quasi-periodic attractor found numerically, while the other one has been discarded (see Remark 5 in Section 4.3).

As far as the slow frequency is concerned, the first attempt would be to study the third approximation and check how the value of the slow frequency changes. However, this is non-trivial because of the form of the tidal torque, which makes the analytical computations rather tricky: going beyond the second approximation requires handling equations which involve functions expressed as integrals, over very long times, of non-smooth functions depending quasi-periodically on time.

As to the discarded second approximation, we conjecture that it does not correspond to any real solution to the equation because either it disappears at some higher step of the iteration or the iteration scheme does not converge in such a case — both possibilities are likely to be difficult to check in practice.

Another issue that deserves further investigation is the presence of barriers in phase space discussed in Section 7. Indeed, it is not obvious why the trajectories cannot cross the resonance $3$:$2$ from above and the resonance $1$:$1$ from below. Even if the $3$:$2$ resonance corresponds to a quasi-periodic solution (see Sections 4 and 5), such a solution is not a KAM torus, which may create an obstruction in phase space (in contrast to what happens in the case of the quasi-synchronous resonance arising in the CTL model); quite the reverse, it is very localised. Moreover nothing similar happens for the other quasi-periodic attractors, such as the $2$:$1$ resonance (the second dominant attractor); in addition the $1$:$1$ resonance corresponds to a periodic solution. All this suggests that the occurrence of the barriers is not due to quasi-periodicity. Rather, it is likely that the phenomenon is related to the amplitude of the peaks appearing in the tidal torque: indeed the largest peaks correspond to the $3$:$2$ and $1$:$1$ resonances.