Secular dynamics of a planar model of the Sun-Jupiter-Saturn-Uranus system; effective stability into the light of Kolmogorov and Nekhoroshev theories

TL;DR

This paper analyzes the long-term stability of the Sun-Jupiter-Saturn-Uranus system using advanced mathematical models, demonstrating that the system remains stable over timescales exceeding the Solar System's lifetime within a refined secular framework.

Contribution

It introduces a refined secular model with second-order mass effects and constructs a Kolmogorov normal form to identify a stable KAM torus, extending stability estimates via Nekhoroshev theory.

Findings

Identifies an invariant KAM torus closely approximating the secular orbits.

Shows the stability time exceeds the Solar System's lifetime within a specific neighborhood.

Determines the neighborhood size is about ten times smaller than observational uncertainties.

Abstract

We investigate the long-time stability of the Sun-Jupiter-Saturn-Uranus system by considering a planar secular model, that can be regarded as a major refinement of the approach first introduced by Lagrange. Indeed, concerning the planetary orbital revolutions, we improve the classical circular approximation by replacing it with a solution that is invariant up to order two in the masses; therefore, we investigate the stability of the secular system for rather small values of the eccentricities. First, we explicitly construct a Kolmogorov normal form, so as to find an invariant KAM torus which approximates very well the secular orbits. Finally, we adapt the approach that is at basis of the analytic part of the Nekhoroshev's theorem, so as to show that there is a neighborhood of that torus for which the estimated stability time is larger than the lifetime of the Solar System. The size of such a neighborhood, compared with the uncertainties of the astronomical observations, is about ten times smaller.

Tables

Table 1. Table 1: Masses m j subscript 𝑚 𝑗 m_{j} and initial conditions for Jupiter, Saturn and Uranus in our planar model. We adopt the AU as unit of length, the year as time unit and set the gravitational constant G = 1 𝐺 1 G=1\, . With these units, the solar mass is equal to ( 2 ​ π ) 2 superscript 2 𝜋 2 (2\pi)^{2} . The initial conditions are expressed by the usual heliocentric planar orbital elements: the semi-major axis a j subscript 𝑎 𝑗 a_{j}\, , the mean anomaly M j subscript 𝑀 𝑗 M_{j}\, , the eccentricity e j subscript 𝑒 𝑗 e_{j} and the perihelion longitude ω j subscript 𝜔 𝑗 \omega_{j}\, . The data are obtained from those reported in Table IV of  [ 50 ] by projecting them on the invariant plane (that is perpendicular to the total angular momentum) in a standard way.

	Jupiter ($j=1$)	Saturn ($j=2$)	Uranus ($j=3$)
$m_j^{}$	${\left(2\pi \right)}^2/1047.355$	${\left(2\pi \right)}^2/3498.5$	${\left(2\pi \right)}^2/22902.98$
$a_j$	$5.20463727204700266$	$9.54108529142232165$	$19.2231635458410572$
$M_j$	$3.04525729444853654$	$5.32199311882584869$	$0.19431922829271914$
$e_j$	$0.04785365972484999$	$0.05460848595674678$	$0.04858667407651962$
${\omega }_j$	$0.24927354029554571$	$1.61225062288036902$	$2.99374344439246487$

Table 2. Table 2: Angular velocities ν ¯ ¯ 𝜈 {\underline{\nu}} and initial conditions ( x ¯ ​ ( 0 ) , y ¯ ​ ( 0 ) ) ¯ 𝑥 0 ¯ 𝑦 0 ({\underline{x}}(0),{\underline{y}}(0)) for our planar secular model about the motions of Jupiter, Saturn and Uranus. The frequency vector ν ¯ ¯ 𝜈 {\underline{\nu}} refers to the harmonic oscillators approximation of the Hamiltonian H ( 𝒟 ) superscript 𝐻 𝒟 H^{(\mathcal{D})} (written in ( 13 )) and its values are given in rad/year.

	$j=1$	$j=2$	$j=3$
${\nu }_j^{}$	$-1.1212724892\times {10}^{-4}$	$-1.9688444678\times {10}^{-5}$	$-1.1134564418\times {10}^{-5}$
$x_j\left(0\right)$	$1.5407573458\times {10}^{-2}$	$-3.0574059274\times {10}^{-2}$	$1.1186486403\times {10}^{-2}$
$y_j\left(0\right)$	$-2.5320810665\times {10}^{-2}$	$-5.2728862107\times {10}^{-3}$	$6.0669645406\times {10}^{-3}$

Equations

$$3n_1-5n_2-7n_3+\left[(3-j)g_1+6g_2+j{g}_3\right],\mathrm{with}j=0,1,2,3,$$

$$F(\tilde{\underline{r}},\underline{r})=T^{(0)}(\tilde{\underline{r}})+U^{(0)}(\underline{r})+T^{(1)}(\tilde{\underline{r}})+U^{(1)}(\underline{r}),$$

$$\begin{array}[]{rclrcl}T^{(0)}(\tilde{{\underline{r}}})&=&\frac{1}{2}\mathrel{\mathop{\kern 0.0pt\sum}\limits_{j=1}^{3}}\frac{m_{0}+m_{j}}{m_{0}m_{j}}\,\|\tilde{{\underline{r}}}_{j}\|^{2}\ ,&T^{(1)}(\tilde{{\underline{r}}})&=&\frac{1}{m_{0}}\Big{(}\tilde{{\underline{r}}}_{1}\cdot\tilde{{\underline{r}}}_{2}+\tilde{{\underline{r}}}_{1}\cdot\tilde{{\underline{r}}}_{3}+\tilde{{\underline{r}}}_{2}\cdot\tilde{{\underline{r}}}_{3}\Big{)}\ ,\cr\cr U^{(0)}({\underline{r}})&=&-G\mathrel{\mathop{\kern 0.0pt\sum}\limits_{j=1}^{3}}\frac{m_{0}\,m_{j}}{\|{\underline{r}}_{j}\|}\ ,&U^{(1)}({\underline{r}})&=&-G\left(\frac{m_{1}\,m_{2}}{\|{\underline{r}}_{1}-{\underline{r}}_{2}\|}+\frac{m_{1}\,m_{3}}{\|{\underline{r}}_{1}-{\underline{r}}_{3}\|}+\frac{m_{2}\,m_{3}}{\|{\underline{r}}_{2}-{\underline{r}}_{3}\|}\right)\ .\cr\end{array}$$

$${\Lambda }_j$$

$$F(\underline{\Lambda },\underline{\lambda },\underline{\xi },\underline{\eta })=F^{(0)}(\underline{\Lambda })+\mu F^{(1)}(\underline{\Lambda },\underline{\lambda },\underline{\xi },\underline{\eta }),$$

$$L_j={\Lambda }_j-{\Lambda }_j^{\ast },\mathrm{for}j=1,2,3.$$

$${\mathcal{H}}^{(\mathcal{T})}(\underline{L},\underline{\lambda },\underline{\xi },\underline{\eta })={\underline{n}}^{\ast }\cdot \underline{L}+\sum _{j_1=2}^{\infty }{h}_{j_1,0}^{(\mathrm{Kep})}(\underline{L})+\mu \sum _{j_1=0}^{\infty }\sum _{j_2=0}^{\infty }{h}_{j_1,j_2}^{(\mathcal{T})}(\underline{L},\underline{\lambda },\underline{\xi },\underline{\eta }),$$

$$\sum _{j=1}^3{n}_j^{\ast }\frac{\partial {\chi }_1^{(\mathcal{O}2)}}{\partial {\lambda }_j}+\sum _{j_2=0}^6{\lceil h_{0,j_2}^{(\mathcal{T})}\rceil }_{\underline{\lambda };8}(\underline{\lambda },\underline{\xi },\underline{\eta })=0.$$

$$\sum _{j=1}^3{n}_j^{\ast }\frac{\partial {\chi }_2^{(\mathcal{O}2)}}{\partial {\lambda }_j}+\sum _{j_2=0}^6{\lceil {\hat{h}}_{1,j_2}^{(\mathcal{O}2)}\rceil }_{\underline{\lambda };8}(\underline{L},\underline{\lambda },\underline{\xi },\underline{\eta })=0.$$

$${\mathcal{H}}^{(\mathcal{O}2)}(\underline{L},\underline{\lambda },\underline{\xi },\underline{\eta })={\underline{n}}^{\ast }\cdot \underline{L}+\sum _{j_1=2}^{\infty }{h}_{j_1,0}^{(\mathrm{Kep})}(\underline{L})+\mu \sum _{j_1=0}^{\infty }\sum _{j_2=0}^{\infty }{h}_{j_1,j_2}^{(\mathcal{O}2)}(\underline{L},\underline{\lambda },\underline{\xi },\underline{\eta };\mu ).$$

$$c_{\underline{j},\underline{k},\underline{r},\underline{s}}{L}_1^{j_1}{L}_2^{j_2}{L}_3^{j_3}{\xi }_1^{r_1}{\xi }_2^{r_2}{\xi }_3^{r_3}{\eta }_1^{s_1}{\eta }_2^{s_2}{\eta }_3^{s_3}\genfrac{}{}{0ex}{}{\sin }{\cos }(k_1{\lambda }_1+k_2{\lambda }_2+k_3{\lambda }_3).$$

$$\left\langle\mathcal{H}^{(\mathcal{O}2)}\big{|}_{{\underline{L}}={\underline{0}}}\,\right\rangle_{{\underline{\lambda}}}=\mu\sum_{j_{2}=0}^{\infty}\big{\langle}h_{0,j_{2}}^{(\mathcal{O}2)}({\underline{\xi}},{\underline{\eta}};\mu)\big{\rangle}_{{\underline{\lambda}}}\ ,$$

$${\mu }^2{h}_{j_1,j_2}^{(q.r.)}(\underline{L},\underline{\lambda },\underline{\xi },\underline{\eta })$$

$$\sum _{j=1}^3{n}_j^{\ast }\frac{\partial {\chi }_1^{(q.r.)}}{\partial {\lambda }_j}+\sum _{j_2=0}^9{h}_{0,j_2}^{(q.r.)}(\underline{\lambda },\underline{\xi },\underline{\eta })=0.$$

$${\mathcal{H}}^{(\sec )}(\underline{\xi },\underline{\eta })$$

$$H^{(\mathcal{D})}(\underline{x},\underline{y})=\sum _{j=0}^3\frac{{\nu }_j}{2}\left(x_j^2+y_j^2\right)+H_2^{(\mathcal{D})}(\underline{x},\underline{y})+H_4^{(\mathcal{D})}(\underline{x},\underline{y})+H_6^{(\mathcal{D})}(\underline{x},\underline{y})+\dots ,$$

$$H^{(\mathrm{I})}(\underline{I},\underline{\phi })=\underline{\nu }\cdot \underline{I}+\sum _{s\ge 1}{f}_{2s}^{(\mathrm{I})}(\underline{I},\underline{\phi }),$$

$$f_{2s}^{(\mathrm{I})}(\underline{I},\underline{\phi })$$

$$H^{(\infty)}({\underline{p}},{\underline{q}})={\underline{\omega}}\cdot{\underline{p}}+\mathcal{O}\big{(}\|{\underline{p}}\|^{2}\big{)}\ ,$$

$$\sum _{j=1}^3{\nu }_j\frac{\partial B^{(\mathrm{II})}}{\partial {\phi }_j}+f_2^{(\mathrm{I})}-⟨f_2^{(\mathrm{I})}{⟩}_{\underline{\phi }}=0.$$

$$H^{(\mathrm{II})}(\underline{I},\underline{\phi })=\underline{\nu }\cdot \underline{I}+f_2^{(\mathrm{II})}(\underline{I})+\sum _{s\ge 2}{f}_{2s}^{(\mathrm{II})}(\underline{I},\underline{\phi }),$$

$$H^{(\mathrm{III})}(\underline{I},\underline{\phi })=\underline{\nu }\cdot \underline{I}+f_2^{(\mathrm{III})}(\underline{I})+f_4^{(\mathrm{III})}(\underline{I})+\sum _{s\ge 3}{f}_{2s}^{(\mathrm{III})}(\underline{I},\underline{\phi }),$$

$$\sum _{j=1}^3{\nu }_j\frac{\partial B^{(\mathrm{III})}}{\partial {\phi }_j}+f_4^{(\mathrm{II})}-⟨f_4^{(\mathrm{II})}{⟩}_{\underline{\phi }}=0;$$

$$I_j=p_j+I_j^{\ast },{\phi }_j=q_j,j=1,\dots ,n.$$

$${\nu }_j+\frac{\partial f_2^{(\mathrm{III})}}{\partial I_j}(\underline{I})+\frac{\partial f_4^{(\mathrm{III})}}{\partial I_j}(\underline{I})={\omega }_j,j=1,\dots ,n.$$

$$H^{(\mathrm{IV})}(\underline{p},\underline{q})=\underline{\omega }\cdot \underline{p}+\sum _{s\ge 0}\sum _{l\ge 0}{f}_l^{(\mathrm{IV},s)}(\underline{p},\underline{q}),$$

$$H^{(r-1)}(\underline{p},\underline{q})=\underline{\omega }\cdot \underline{p}+\sum _{s\ge 0}\sum _{l\ge 0}{f}_l^{(r-1,s)}(\underline{p},\underline{q}),$$

$$\sum _{j=1}^n{\omega }_j\frac{\partial X^{(r)}}{\partial q_j}(\underline{q})+\sum _{s=1}^r{f}_0^{(r-1,s)}(\underline{q})=0,C^{(r)}{\underline{\xi }}^{(r)}\cdot \underline{p}+f_1^{(r-1,0)}(\underline{p})=0,$$

$${\hat{H}}^{(r)}(\underline{p},\underline{q})=\underline{\omega }\cdot \underline{p}+\sum _{s\ge 0}\sum _{l\ge 0}{\hat{f}}_l^{(r,s)}(\underline{p},\underline{q}),$$

$${\hat{f}}_l^{(r,s)}=f_l^{(r-1,s)}(\underline{p},\underline{q}),\mathrm{for}l\ge 0\mathrm{and}s\ge 0.$$

Full text

1 Introduction

The problem that we would like to investigate in this paper is concerned with the long-time stability of the Solar System. However taking into account all planets (possibly including also the satellites) is obviously an overwhelming task. Therefore we limit our efforts to considering a simplified model including three giant planets, namely studying the Sun-Jupiter-Saturn-Uranus (SJSU) system. Moreover we make a further simplification by considering the planar model. In the spirit of a long-time research that we are conducting, we investigate the stability in Nekhoroshev sense in the neighbourhood of an approximated KAM torus close to the real orbits of the planets.

To this end, we first make a short summary of the relevant information concerning KAM and Nekhoroshev theories. Then we give an account of our results.

1.1 General framework

In his celebrated article appeared in 1954 (see [16]), Kolmogorov proved the persistence of quasi-periodic motions on invariant tori for nearly-integrable Hamiltonian systems. In the same paper he pointed out the possible applications of his theorem to the problem of the stability of the Solar System. More generally, Celestial Mechanics appeared as a main field of applications of Kolmogorov theory. Nowadays, there is no doubt that the suggestion of Kolmogorov was prophetic.

Perhaps the first application is concerned with the stability of the equilateral Lagrangian points in the planar circular restricted three-body problem (hereafter PCR3BP, see [18]). A little later the papers of Moser and Arnold (see [41] and [1]) gave the first detailed proof of Kolmogorov’s theorem, that was only sketched in the original 1954 paper. The three above mentioned works marked the beginning of the so-called KAM theory.

The relevance of KAM theory for Celestial Mechanics was emphasized also by Moser (see [42]). In particular he recalled the old idea, strongly supported by Weierstrass, that the trigonometric perturbation series of Celestial Mechanics should be actually convergent, thus proving that the motion of the planets is quasi-periodic. That is, essentially described by the old-fashioned epicycles of the Greek astronomy.

KAM theory has been used to prove the stability of some particular systems: the spin-orbit problem, the motion of the asteroid Victoria in the framework of the PCR3BP, the secular dynamics of the Sun-Jupiter-Saturn system, etc. (see [6], [7], [36], respectively). All these models have two degrees of freedom, therefore their stability is proved by a topological argument. The key point is that an invariant torus is a two-dimensional manifold that separates the three-dimensional energy surface. Therefore, orbits starting in the region between two invariant KAM tori on the same energy surface, remain perpetually confined there. For more degrees of freedom, $n>2$, one is confronted with the problem that the set of invariant KAM tori has positive Lebesgue measure, but empty interior. Therefore, the energy manifold is ($2n-1$)-dimensional, while the invariant tori are $n$-dimensional, so that they cannot act as barriers for the motion. This opens the problem of the so-called Arnold diffusion, the generic existence of which has been recently proven near a resonance of codimension one (see [3] and references therein).

A different approach to the problem of stability may be based on the theorem of Nekhoroshev (see [44] and [45]). The difference may be summarized as follows. In KAM theory one proves perpetual stability only for a large (Cantor-like) set of initial conditions, thus renouncing to consider all possible configurations. In Nekhoroshev theory one looks for a result valid for all initial data (e.g., in an open set), but accepts stability for a finite time, asking that time to be very long. Actually, this kind of approach was already proposed by Moser and Littlewood (see [40], [34] and [35]) in a local approach around an equilibrium point. Their theory essentially constitutes the so-called analytic part of the proof of Nekhoroshev’s theorem. The remarkable contribution of Nekhoroshev is represented by the geometric part, which allows him to describe the dynamics in the whole phase space.

Let us rephrase the statement of Nekhoroshev’s theorem in an informal way. The theorem is concerned with the general problem of dynamics, so-named by Poincaré. One considers an analytic Hamiltonian $H({\underline{p}},{\underline{q}})=h({\underline{p}})+\varepsilon f({\underline{p}},{\underline{q}})$ (in action-angle variables $({\underline{p}},{\underline{q}})\in\mathcal{U}\times\mathbb{T}^{n}$, with open $\mathcal{U}\subset\mathbb{R}^{n}$).

If $\varepsilon$ is small enough and the unperturbed Hamiltonian $h({\underline{p}})$ satisfies an appropriate steepness hypothesis, then for every orbit with initial condition \big{(}{\underline{p}}(0),{\underline{q}}(0)\big{)}\in\mathcal{U}\times\mathbb{T}^{n} the bound $\|{\underline{p}}(t)-{\underline{p}}(0)\|<\varepsilon^{b}$ holds true for $|t|\leq T(\varepsilon)$ where $T(\varepsilon)\sim\exp(1/\varepsilon^{a})$, for some constants $a,b\in(0,1)$.

Let us add some considerations concerning the application of the concept of stability over finite, large time to physical systems. In this respect, the theory of Moser, Littlewood and Nekhoroshev can be seen as an evolution of the old-fashioned adiabatic theory. The underlying idea is that the perpetual stability in Lyapunov sense is a too strong property, which can hardly be proved. In physical applications, including a planetary system, it is enough to prove stability for a time interval comparable with the lifetime of the system itself. E.g., for the Solar System the estimated age of the Universe is enough. This means that we should investigate stability up to a time of about $10^{10}$ years.

Giving the adjective “long” a definite mathematical sense is of course a more difficult task. We can only rely on the dependence of $T(\varepsilon)$ on the perturbation parameter $\varepsilon$. With reference to history we can collect a short list of attempts.

1. i.

Adiabatic theory means $T(\varepsilon)\sim 1/\varepsilon$. This is a concept that has been widely investigated in physics and played a relevant role in the development of Quantum Theory. 2. ii.

Birkhoff complete stability means $T(\varepsilon)\sim 1/\varepsilon^{r}$ for some $r>1$. The concept has been proposed by Birkhoff in [5]. 3. iii.

Exponential stability means $T(\varepsilon)\sim\exp(1/\varepsilon^{a})$ as in the statement above. 4. iv.

Super-exponential stability means T(\varepsilon)\sim\exp\bigl{(}\exp(1/\varepsilon^{a})\bigr{)}. This has been proposed by Morbidelli and one of the authors in [39]. A further improvement by the same authors shows that in appropriate boxed subsets of the phase space one finds T(\varepsilon)\sim\exp\bigl{(}\exp\bigl{(}\ldots\exp(1/\varepsilon^{a})\bigr{)}\bigr{)} with an increasing number of exponentials. The limit of the boxed subsets appears to be connected with the set of invariant KAM tori, see [11].

The possible applications to the real world deserve a careful detailed discussion. Unlike the mathematical approach, we must face the fact that the size of the perturbation parameter is fixed by Nature. Therefore the question may be formulated as follows: given that $T(\varepsilon)$ has some definite behaviour as $\varepsilon$ goes to [math], what can we say for a specific system with a given value of $\varepsilon$? We refer to this concept as effective stability, in the sense that we aim to prove that $T(\varepsilon)$, for that given value of $\varepsilon$, is large enough to cover some characteristic time of the physical system, e.g., its lifetime, as we have already said.

The relevant fact in this connection is that the analytic form of $T(\varepsilon)$ provided by analytical theories appears to be just a smoothing of a more complex behavior of the estimated stability time, as it can be found with the help of computer algebra or similar methods. Let us explain this fact in the framework, e.g., of the stability of an elliptic equilibrium. The procedure goes through the calculation of the Birkhoff normal form up to a finite order $r$. Birkhoff’s theory of complete stability states that in a neighbourhood of radius $\varrho>0$ of the equilibrium we get $T(\varrho)\sim 1/(C_{r}\varrho^{r})$, with a constant $C_{r}$ that Birkhoff did not evaluate. Here the natural perturbation parameter is the distance $\varrho$ from the equilibrium, which takes the place of $\varepsilon$ in the general statements above. The crucial problem is the relation between $\varrho$ and $r$. What we can actually find is a function $\widetilde{T}(\varrho,r)$ depending on both parameters $\varrho$ and $r$. On the other hand, having fixed $\varrho$, we are allowed to make the best choice of $r$ as a function of $\varrho$ so as to maximize $\widetilde{T}(\varrho,r)$. In the case of the elliptic equilibrium one usually finds $C_{r}\sim(r!)^{c}$ with $c\geq 1$, i.e., $\widetilde{T}(\varrho,r)\sim((r!)^{c}\varrho^{r})^{-1}$, and the optimal choice $r\sim(1/\varrho)^{1/c}$ produces the exponential estimate $T(\varrho)\sim\exp(1/\varrho^{a})$ with $a=1/c$. More precisely one finds that there is an increasing sequence $\varrho_{1},\varrho_{2},\varrho_{3},\ldots$ of values of $\varrho$ such that in every interval $(\varrho_{j},\varrho_{j+1})$ one gets $T(\varrho)\sim 1/\varrho^{\,j}$. In this respect the theory of Moser, Littlewood and Nekhoroshev appears to bound the latter function $T(\varrho)$ from below.

In a practical application, if we are able to explicitly calculate the Birkhoff normal form up to some maximal order $r$, e.g., using computer algebra, then we can actually draw the function $T(\varrho)$ as the sequence of optimal values of $\widetilde{T}(\varrho,r_{j})$ in different intervals, as we do later in Figure 2.

Let us give a short historical account on the applications of the methods above in Celestial Mechanics. Most of them are concerned with the dynamics of Trojan asteroids. In [15], a few asteroids have been shown to be effectively stable (over the age of the Universe) in the framework of the PCR3BP, where Sun and Jupiter played the role of the primary bodies on circular orbits. Such an approach is not limited to models having two degrees of freedom; in fact, it has been extended to the spatial case (see [49]) and to the elliptic one (see [31], where the dynamics is represented by a four-dimensional symplectic map instead of using a continuous Hamiltonian flow). Let us recall that all these results are based on the explicit construction of the Birkhoff normal form using computer algebra. By the way the case of the Lagrangian points is precisely the one studied by Littlewood where an estimate similar to the exponential one by Nekhoroshev has been found. He commented: “while not eternity, this is a considerable slice of it”.

A similar approach allows us also to extend the theory of Lagrange and Laplace for the secular motion of the longitudes of the perihelia and nodes of the planets. Indeed, following the traditional approach, we may introduce the secular approximation by assuming that three semi-major axes remain invariant up to order two in the masses. Then the Hamiltonian for the eccentricities and inclinations, with the conjugate longitudes of the perihelia and nodes, may be written in Poincaré variables as a system around an elliptic equilibrium (see [46]). The equilibrium in this case corresponds to planar circular orbits. Therefore we may investigate the effective stability of the planets by just extending the method used for the Lagrangian equilibria. It should be noted that the complexity of the problem becomes much larger in view of the increasing number of degrees of freedom. We stress however that the difficulty is a mere technical one, due to the limitations of memory and computer power.

This method has been applied by the authors to a planar model of the SJSU system (see [48]). However the results appear neither realistic nor very promising due to the actual values of the eccentricities of the planets. Indeed, in that paper we have shown that the eccentricities of the planets are too large (about twice), in order to ensure the effective stability of the secular dynamics over the age of the Universe. We remark that this is precisely the model that we investigate in the present paper, with an additional improvement that we are going to describe.

A productive combination of KAM and Nekhoroshev theories consists in applying the usual, local theory for an elliptic equilibrium to the neighbourhood of an invariant Kolmogorov torus. In such a neighbourhood, Kolmogorov procedure produces a Hamiltonian that may be given the form $H({\underline{p}},{\underline{q}})={\underline{\omega}}\cdot{\underline{p}}+h_{1}({\underline{p}},{\underline{q}})+h_{2}({\underline{p}},{\underline{q}})+\ldots$, expanded in power series of the actions ${\underline{p}}$ ($h_{i}$ being of degree $i+1$). This is indeed similar to the form of the Hamiltonian in a neighbourhood of an elliptic equilibrium in action-angle variables. The apparently strong difference is that in the case of elliptic equilibrium we are dealing with trigonometric polynomials $h_{i}$ of increasing orders, while in the case of the torus $h_{i}$ is an infinite trigonometric series. However this represents a minor problem, indeed we can suitably arrange the Hamiltonian as an expansion in trigonometric polynomials, exploiting the exponential decay of Fourier coefficients. Therefore the analytic theory for the elliptic equilibrium applies almost verbatim to the neighbourhood of a torus. This remark suggests that the long-time stability of an elliptic equilibrium and of a torus can be investigated using the same method.

We have applied the latter idea to the Sun-Jupiter-Saturn system (see [12]), where we used a previous result on the existence of a torus of the SJS system (see [38]) based on the explicit expansion of the Hamiltonian and on the explicit application of Kolmogorov method up to a finite, not too low order. Then we worked out a Birkhoff normalization and showed that there is a domain of effective stability, which is centered around an invariant KAM torus. The results were close to realistic ones.

In view of the previous experience, we decided to work out the application to the planar SJSU system, exploiting the same idea of making expansions around an approximated KAM torus. This is indeed the goal of the present paper.

1.2 Plan of the work

As it is well known, the major problems in perturbation theory arise from the existence of resonances in the trigonometric expansions. It goes without saying that a reliable theory should positively take into account these resonant terms. In the case of the SJSU system the main resonances have been described by Murray and Holman. With the numerical exploration made in [43], they have pointed out the dynamical mechanism inducing a slightly chaotic component in the motion of the major bodies of our planetary system. Actually, this phenomenon is due to the overlap of some resonances involving three or four bodies. An example is given by the resonances

[TABLE]

where $n_{i}$ stands for the mean motion frequency of the $i$-th planet, $g_{i}$ means the (secular) frequency of its perihelion argument and the labels $1,\,2,\,3$ refer to Jupiter, Saturn and Uranus, respectively. In fact, during the planetary motion each angle corresponding to the resonances above jumps from libration to rotation and vice versa. Many other resonances analogous to the previous ones are located in the vicinity of the real orbit of the system including the Jovian planets, some of them involving also Neptune and the frequencies related to the longitudes of the nodes. In the same article a rather simplistic argument is provided so as to evaluate the time needed by these resonances to eject Uranus from the Solar System, that is estimated to be about $10^{18}$ years. One should also recall that the dynamics of the terrestrial planets is much more chaotic: collisions between Mercury, Mars or, even, Venus with the Earth could take place in about $3\,\times 10^{9}$ years (see [29]).

Our procedure is essentially based on two steps. First, we explicitly perform the construction of the Kolmogorov normal form for the planar secular model of the SJSU system. In this first step, the expansions of the Hamiltonians introduced by the normalization algorithm are computed by using a software specially designed for doing algebraic manipulations (see [14]). In the second step, we avoid the explicit expansions (due to memory and power limits of our computers) by setting up a suitable scheme of estimates for the norms of the functions. Precisely, we replace the explicit construction of the Birkhoff normal form around the invariant KAM torus with a recursive scheme of estimates on the norms. The results so produced are a little worse with respect to an explicit computation of the series, but nevertheless the final results are acceptably close to be realistic. On the other hand associating to every function a norm (i.e., a number) instead of a trigonometric polynomial obviously makes the calculation definitely faster, while allowing to reach much higher orders.

The whole procedure, including the optimization of the estimates with an optimal choice of the expansion order $r$, allows us to evaluate a lower bound for the stability time $T(\varrho_{0})$. The final result is the following. We find a ball of effective stability over the age of the Universe having radius $\varrho_{0}$ and center on the previously found KAM torus. The value of $\varrho_{0}$ is meaningful from a physical point of view. Indeed, considering a ball of initial conditions that takes into account the uncertainties of the astronomical observations, we find a value of $\varrho_{0}$ which is about ten times smaller. Our result is not supported by a rigorous computer-assisted proof (see, e.g., [7]), but we think that this could be done with some additional effort, similarly to what we did in the past (see [36]).

The paper is organized as follows. In order to make the work rather self-consistent, the model is introduced in section 2, where we put a particular care in pointing out some technical difficulties. Section 3 is devoted to the construction of the Kolmogorov normal form, so as to find an invariant KAM torus which approximates very well the secular orbits of our planetary model. In section 4 we perform the search for stability in the neighbourhood of the KAM torus.

2 Settings for the definition of the Hamiltonian model

For the sake of definiteness, in the present section we recall the basic steps that are necessary to introduce the same planar secular model of the Sun-Jupiter-Saturn-Uranus system that was already studied in [48]. We defer to sects. 2 and 3 of that paper for more details.

2.1 Classical expansion of the planar

planetary Hamiltonian

Let us consider four point bodies $P_{0},\,P_{1},\,P_{2},\,P_{3}$, with masses $m_{0},\,m_{1},\,m_{2},\,m_{3}$, mutually interacting according to Newton’s gravitational law. Hereafter the indexes $0,\,1,\,2,\,3$ will correspond to Sun, Jupiter, Saturn and Uranus, respectively. We basically follow the formalism introduced by Poincaré (see, e.g., [26] and [30] for a modern exposition). We remove the motion of the center of mass by using heliocentric coordinates ${\underline{r}}_{j}=\mathrel{\mathop{\kern 0.0ptP_{0}P_{j}}\limits^{\longrightarrow}}\,$, with $j=1,\,2,\,3\,$. Denoting by $\tilde{{\underline{r}}}_{j}$ the momenta conjugate to ${\underline{r}}_{j}$, the Hamiltonian of the system has $6$ degrees of freedom, and reads

[TABLE]

where

[TABLE]

The plane set of Poincaré’s canonical variables is defined as

[TABLE]

for $j=1\,,\,2\,,\,3\,$, where $a_{j}\,,\>e_{j}\,,\>M_{j}$ and $\omega_{j}$ are the semi-major axis, the eccentricity, the mean anomaly and the perihelion longitude, respectively, of the $j$-th planet. Let us remark that both $\xi_{j}$ and $\eta_{j}$ are of the same order of magnitude as the eccentricity $e_{j}\,$. Using Poincaré’s variables (2), the Hamiltonian $F$ can be rearranged so that one has

[TABLE]

where $F^{(0)}=T^{(0)}+U^{(0)}$, $\mu F^{(1)}=T^{(1)}+U^{(1)}$. Here, the small dimensionless parameter $\mu=\max\{m_{1}\,/\,m_{0}\,$, $\,m_{2}\,/\,m_{0}\,$, $\,m_{3}\,/\,m_{0}\,\}$ has been introduced in order to highlight the different size of the terms appearing in the Hamiltonian. According to the common language in Celestial Mechanics, in the following we will refer to ${\underline{\lambda}}$ and to their conjugate actions ${\underline{\Lambda}}$ as the fast variables, while $({\underline{\xi}},{\underline{\eta}})$ will be called secular variables.

We proceed now by expanding the Hamiltonian (3) in order to construct the first basic approximation of Kolmogorov normal form. We pick a value ${\underline{\Lambda}}^{*}$ for the fast actions and perform a translation $\mathcal{T}_{{\underline{\Lambda}}^{*}}$ defined as

[TABLE]

This is a canonical transformation that leaves the coordinates ${\underline{\lambda}}\,$, ${\underline{\xi}}$ and ${\underline{\eta}}$ unchanged. The transformed Hamiltonian $\mathcal{H}^{(\mathcal{T})}=F\circ\mathcal{T}_{{\underline{\Lambda}}^{*}}\,$ can be expanded in power series of ${\underline{L}},\,{\underline{\xi}},\,{\underline{\eta}}$ around the origin. Thus, forgetting an unessential constant we rearrange the Hamiltonian of the system as

[TABLE]

where the functions $h_{j_{1},j_{2}}^{(\mathcal{T})}$ are homogeneous polynomials of degree $j_{1}$ in the actions ${\underline{L}}$ and of degree $j_{2}$ in the secular variables $({\underline{\xi}},{\underline{\eta}})\,$. The coefficients of such homogeneous polynomials do depend analytically and periodically on the angles ${\underline{\lambda}}\,$. The terms $h_{j_{1},0}^{({\rm Kep})}$ of the Keplerian part are homogeneous polynomials of degree $j_{1}$ in the actions ${\underline{L}}\,$, the explicit expression of which can be determined in a straightforward manner. In the latter equation the term which is both linear in the actions and independent of all the other canonical variables (i.e., ${\underline{n}}^{*}\cdot{\underline{L}}$) has been separated in view of its relevance in perturbation theory, as it will be discussed in the next subsection. We also expand the coefficients of the power series $h_{j_{1},j_{2}}^{(T_{F})}$ in Fourier series of the angles ${\underline{\lambda}}\,$, by following a traditional procedure in Celestial Mechanics. We work out these expansions for the case of the planar SJSU system using a specially devised algebraic manipulation (see [14] for an introduction to the main ideas that have been translated in our codes).

We now describe how to determine the fixed values ${\underline{\Lambda}}^{*}$ that allows us to perform the expansion (5) of the Hamiltonian as a function of the canonical coordinates $({\underline{L}},{\underline{\lambda}},{\underline{\xi}},{\underline{\eta}})$. To this end we perform a long-term numerical integration of Newton’s equations starting from the initial conditions related to the data reported in Table 1. After having computed the average values $(a_{1}^{*}\,,\,a_{2}^{*}\,,\,a_{3}^{*})$ of the semi-major axes during the evolution, we determine the values ${\underline{\Lambda}}^{*}$ via the first equation in (2). In our calculations we truncate the expansion as follows. (a) The Keplerian part is expanded up to the quadratic terms. The terms $h_{j_{1},j_{2}}^{(\mathcal{T})}$ include: (b1) the linear terms in the actions ${\underline{L}}\,$, (b2) all terms up to degree $18$ in the secular variables $({\underline{\xi}},{\underline{\eta}})\,$, (b3) all terms up to the trigonometric degree $16$ with respect to the angles ${\underline{\lambda}}\,$. Our choice of the limits will be fully motivated in the next subsection.

2.2 The secular model

We now introduce the rather accurate description of the secular dynamics provided by the average of the Hamiltonian up to order two in the masses (see, e.g., [25], [27], [36], [17], [32] and [33]). To this end we follow the approach described in [38], carrying out two “Kolmogorov-like” normalization steps in order to eliminate the main perturbation terms depending on the fast angles ${\underline{\lambda}}\,$. We concentrate our attention on the quasi-resonant angles $2\lambda_{1}-5\lambda_{2}\,$, $\lambda_{1}-7\lambda_{3}$ and $3\lambda_{1}-5\lambda_{2}-7\lambda_{3}$, which are the most relevant ones for the dynamics. The procedure is a little cumbersome, and requires two main steps that we describe in the following subsections.

2.2.1 Partial reduction of the perturbation

We emphasize that the Fourier expansion of the Hamiltonian (5) is generated just by terms due to two-body interactions, and so harmonics including more than two fast angles cannot appear. Thus, at first order in the masses, only harmonics with the quasi-resonant angles $2\lambda_{1}-5\lambda_{2}$ and $\lambda_{1}-7\lambda_{3}$ do occur. Actually, harmonics with the quasi-resonant angle $3\lambda_{1}-5\lambda_{2}-7\lambda_{3}$ are generated by the first Kolmogorov-like transformation, but are of second order in the masses, and should be removed by the second Kolmogorov-like transformation described in the next section.

Let us go into details. We denote by \big{\lceil}f\big{\rceil}_{{\underline{\lambda}};K_{F}} the Fourier expansion of a function $f$ truncated so as to include only its harmonics ${\underline{k}}\cdot{\underline{\lambda}}$ satisfying the restriction $0<|{\underline{k}}|\leq K_{F}\,$, with some fixed $K_{F}$, being $|{\underline{k}}|=|k_{1}|+|k_{2}|+|k_{3}|\,$. We also denote by $\langle\cdot\rangle_{{\underline{\lambda}}}$ the average with respect to the angles $\lambda_{1}\,$, $\lambda_{2}$ and $\lambda_{3}\,$. The canonical transformations are using the Lie series algorithm (see, e.g., [10]).

We set $K_{F}=8$ and transform the Hamiltonian (5) as $\hat{\mathcal{H}}^{(\mathcal{O}2)}=\exp\mathcal{L}_{\mu\,\chi_{1}^{(\mathcal{O}2)}}\,\mathcal{H}^{(\mathcal{T})}$ with the generating function $\mu\,\chi_{1}^{(\mathcal{O}2)}({\underline{\lambda}},{\underline{\xi}},{\underline{\eta}})$ determined by solving the equation

[TABLE]

Notice that, by definition, \big{\langle}\big{\lceil}f\big{\rceil}_{{\underline{\lambda}};K_{F}}\big{\rangle}_{{\underline{\lambda}}}=0\,, which assures that equation (6) can be solved provided the frequencies $n_{1}^{*}\,,\,n_{2}^{*}$ and $n_{3}^{*}$ are non-resonant up to order $8\,$, as it actually occurs in our planar model of the SJSU system. The Hamiltonian $\hat{\mathcal{H}}^{(\mathcal{O}2)}$ has the same form of $\mathcal{H}^{(\mathcal{T})}$ in (5), with the functions $h_{j_{1},j_{2}}^{(\mathcal{T})}$ replaced by new ones, that we denote by $\hat{h}_{j_{1},j_{2}}^{(\mathcal{O}2)}$, generated by expanding the Lie series $\exp\mathcal{L}_{\mu\,\chi_{1}^{(\mathcal{O}2)}}\,\mathcal{H}^{(\mathcal{T})}$ and by gathering all the terms having the same degree both in the fast actions and in the secular variables.

Now we perform a second canonical transformation $\mathcal{H}^{(\mathcal{O}2)}=\exp\mathcal{L}_{\mu\,\chi_{2}^{(\mathcal{O}2)}}\,\hat{\mathcal{H}}^{(\mathcal{O}2)}\,$, where the generating function $\mu\,\chi_{2}^{(\mathcal{O}2)}({\underline{L}},{\underline{\lambda}},{\underline{\xi}},{\underline{\eta}})$ (which is linear with respect to ${\underline{L}}$) is determined by solving the equation

[TABLE]

Again, the Hamiltonian $\mathcal{H}^{(\mathcal{O}2)}$ can be written in a form similar to (5), namely

[TABLE]

where the new functions $h_{j_{1},j_{2}}^{(\mathcal{O}2)}$ are calculated as previously explained for $\hat{h}_{j_{1},j_{2}}^{(\mathcal{O}2)}\,$. Moreover, they still have the same dependence on their arguments as $h_{j_{1},j_{2}}^{(\mathcal{T})}$ in (5).

If terms of second order in $\mu$ are neglected, then the Hamiltonian $\mathcal{H}^{(\mathcal{O}2)}$ possesses the secular three-dimensional invariant torus ${\underline{L}}={\underline{0}}$ and ${\underline{\xi}}={\underline{\eta}}={\underline{0}}$. Thus, in a small neighborhood of the origin of the translated fast actions and for small eccentricities the solutions of the system with Hamiltonian $\mathcal{H}^{(\mathcal{O}2)}$ differ from those of its average $\langle\mathcal{H}^{(\mathcal{O}2)}\rangle_{{\underline{\lambda}}}$ by a quantity $\mathcal{O}(\mu^{2})$. In this sense the average of the Hamiltonian (8) approximates the real dynamics of the secular variables up to order two in the masses, and due to the choice $K_{F}=8$ takes into account the quasi-resonances $5:2$ between Jupiter and Saturn and $7:1$ between Jupiter and Uranus.

In this part of the calculation we produce a truncated series which is represented as a sum of monomials

[TABLE]

The truncated expansion of $\mathcal{H}^{(\mathcal{O}2)}$ contains $94\,109\,751$ such monomials. We truncate our expansion at degree 16 in the fast angles ${\underline{\lambda}}$ (keeping all harmonics) and at degree 18 in the slow variables ${\underline{\xi}},\,{\underline{\eta}}$ (we shall justify this choice at the end of the next section).

2.2.2 Second approximation and reduction to the secular

Hamiltonian

Since we plan to consider the secular system, we perform a partial average by keeping only the main terms that contain the quasi-resonant angle $3\lambda_{1}-5\lambda_{2}-7\lambda_{3}$. More precisely, we first consider the reduced Hamiltonian

[TABLE]

namely we set ${\underline{L}}={\underline{0}}$ and average $\mathcal{H}^{(\mathcal{O}2)}$ by removing all the Fourier harmonics depending on the angles. Next, we select in $\mathcal{H}^{(\mathcal{O}2)}$ the Fourier harmonics that contain the wanted quasi-resonant angle $3\lambda_{1}-5\lambda_{2}-7\lambda_{3}$ and add them to the Hamiltonian (9). Finally, we perform on the resulting Hamiltonian the second Kolmogorov-like step. With more detail, this is the procedure, which is an adaptation of a scheme already used in [36]. For $(j_{1},j_{2})\in\mathbb{N}^{2}$ we select the quasi-resonant terms

[TABLE]

Actually, this means that in our expression we just remove all monomials but the ones containing the wanted quasi-resonant angle. Using the selected terms we determine a generating function $\mu^{2}\chi_{1}^{({\rm q.r.})}({\underline{\lambda}},{\underline{\xi}},{\underline{\eta}})$ by solving the equation

[TABLE]

Here we make the calculation faster by keeping only terms up to degree $9$ in $({\underline{\xi}},{\underline{\eta}})\,$, this allows us to keep the more relevant quasi-resonant contributions. Then, we calculate only the interesting part of the transformed Hamiltonian $\exp\mathcal{L}_{\mu^{2}\,\chi_{1}^{({\rm q.r.})}}\,\mathcal{H}^{(\mathcal{O}2)}\,$, namely we keep in the transformation only the part which is independent of all the fast variables $({\underline{L}},{\underline{\lambda}})\,$. This produces the secular Hamiltonian $\mathcal{H}^{({\rm sec})}$, independent of ${\underline{\lambda}}$, which satisfies the formal equation \big{\langle}\exp\mathcal{L}_{\mu^{2}\,\chi_{1}^{({\rm q.r.})}}\,\mathcal{H}^{(\mathcal{O}2)}\big{\rangle}_{{\underline{\lambda}}}=\mathcal{H}^{({\rm sec})}+\mathcal{O}(\|{\underline{L}}\|)+o(\mu^{4})\,, where

[TABLE]

Here, we denoted by $\left\{\cdot,\cdot\right\}_{{\underline{L}},{\underline{\lambda}}}$ and $\left\{\cdot,\cdot\right\}_{{\underline{\xi}},{\underline{\eta}}}$ the terms of the Poisson bracket involving only the derivatives with respect to the conjugate variables $({\underline{L}},{\underline{\lambda}})$ and $({\underline{\xi}},{\underline{\eta}})$, respectively.

The Hamiltonian so constructed is the secular one, describing the slow motion of eccentricities and perihelia. In view of the d’Alembert rules $\mathcal{H}^{({\rm sec})}$ contains only terms of even degree and so the lowest order significant term has degree 2. By using our specially designed algebraic manipulator, we have computed the power series expansion of the secular Hamiltonian up to degree 18 in the slow variables. In order to allow comparisons with other expansions, our results up to degree 4 in $({\underline{\xi}},{\underline{\eta}})$ are reported in appendix A of [48].

Let us close this section with a few remarks which justify our choice of the truncation orders. The limits on the expansions in the fast actions ${\underline{L}}$ have been illustrated at points (a) and (b1) at the end of sect. 2.1, and they are the smallest ones that are required in order to make the Kolmogorov-like normalization procedure significant. Since we want to keep the quasi-resonant angles $2\lambda_{1}-5\lambda_{2}\,$, $\lambda_{1}-7\lambda_{3}$ and $3\lambda_{1}-5\lambda_{5}-7\lambda_{3}\,$, we should set the truncation order for Fourier series to $16$, which is enough. The choice to truncate the polynomial expansion at degree $18$ in the secular variables $({\underline{\xi}},{\underline{\eta}})$ is somehow subtler. In view of d’Alembert rules the harmonics $2\lambda_{1}-5\lambda_{2}$ and $\lambda_{1}-7\lambda_{3}$ have coefficients of degree at least 3 and 6, respectively, in the secular variables. Thus, we decided to calculate the generating functions $\chi_{1}^{(\mathcal{O}2)}$ and $\chi_{2}^{(\mathcal{O}2)}$ up to degree $6$ (recall equations (6) and (7)). Furthermore, the quasi-resonant angle $3\lambda_{1}-5\lambda_{5}-7\lambda_{3}$ does not appear initially in the Hamiltonian, but is generated by Poisson bracket between the harmonics $2\lambda_{1}-5\lambda_{2}$ and $\lambda_{1}-7\lambda_{3}$, which produces monomials of degree $9$ in $({\underline{\xi}},{\underline{\eta}})$. Therefore, we decided to calculate the generating function $\chi_{1}^{({\rm q.r.})}$ up to degree $9$ (recall equation (11)). Finally, in the second Kolmogorov-like step we want to keep the secular terms generated by the harmonic $3\lambda_{1}-5\lambda_{5}-7\lambda_{3}$, which are produced by Poisson bracket between monomials containing precisely this harmonic, and then the result has maximum degree $18$ in $({\underline{\xi}},{\underline{\eta}})$. This justifies the final truncation order for the slow variables.

3 Normalization algorithm constructing invariant tori close to

an elliptic equilibrium point

The lowest order approximation of the secular Hamiltonian $\mathcal{H}^{({\rm sec})}$, namely its quadratic term, is essentially the one considered in the theory first developed by Lagrange (see [19]), later extended by Laplace (see [22], [23] and [24]) and further improved by Lagrange himself (see [20], [21]). In modern language, we say that the origin of the reduced phase space (i.e., $({\underline{\xi}},{\underline{\eta}})=({\underline{0}},{\underline{0}})\,$) is an elliptic equilibrium point (for a review using a modern formalism, see sect. 3 of [4], where a planar model of our Solar System is considered).

It is well known that (under mild assumptions on the quadratic part of the Hamiltonian which are satisfied in our case) one can find a linear canonical transformation $({\underline{\xi}},{\underline{\eta}})=\mathcal{D}({\underline{x}},{\underline{y}})$ which diagonalizes the quadratic part of the Hamiltonian, so that we can write $\mathcal{H}^{({\rm sec})}$ in the new coordinates as

[TABLE]

where $\nu_{j}$ are the secular frequencies in the small oscillations limit and $H_{2s}^{(\mathcal{D})}$ is a homogeneous polynomial of degree $2s+2$ in $({\underline{x}},{\underline{y}})\,$. The calculated values of $\nu_{1},\nu_{2}$ and $\nu_{3})$ are reported in Table 2.

At this point, it is convenient to introduce a further preliminary canonical transformation $({\underline{x}},{\underline{y}})=\mathcal{A}({\underline{I}},{\underline{\varphi}})$ so as to introduce action-angle coordinates, namely $x_{j}=\sqrt{2I_{j}}\cos\varphi_{j}$ and $y_{j}=\sqrt{2I_{j}}\sin\varphi_{j}$, for $j=1,\,2,\,3$. Then, the expansion of the new Hamiltonian $H^{({\rm I})}=H^{(\mathcal{D})}\circ\mathcal{A}$ is (sometimes said) of d’Alembert type, i.e., it can be written as

[TABLE]

where the functions $f_{2s}^{({\rm I})}$ are homogeneous polynomials of degree $2s+2$ in the square root of the actions and trigonometric polynomials of degree $2s+2$ in the angles. Moreover, for any fixed index $m\in\{1,\,\ldots,\,n\}$ (obviously being, in our case, the number of degrees of freedom $n=3$) and for all term appearing in the expansion of $f_{2s}^{({\rm I})}$ the $m$-th component of the Fourier harmonics is not greater than the corresponding degree in $\sqrt{I_{m}}$ and they have the same parity, i.e.,

[TABLE]

where $c_{i_{1},\ldots,i_{n},j_{1},\ldots,j_{n}}^{({\rm I})}$ and $d_{i_{1},\ldots,i_{n},j_{1},\ldots,j_{n}}^{({\rm I})}$ are real coefficients.

In principle, the form of the expansion (14) of the Hamiltonian $H^{({\rm I})}$ would be perfectly suitable to perform the procedure constructing invariant tori near an elliptic equilibrium point as it is described in [36]. Unfortunately, in our model of the planetary problem, the initial conditions reported in Table 2 are far enough from the equilibrium point to induce some effects of numerical instability that are mainly due to the determination of the actions translations (which will be properly defined in the following); this can prevent the convergence of the algorithm constructing the KAM torus. In order to circumvent such an obstruction, it is better to adapt to the present context the normalization scheme introduced in [9], where the previous and more usual approach is refined so as to produce an algorithm that is effective also in a region of the phase space not so close to the elliptic equilibrium point. In order to make this work rather self-consistent, it is convenient to first recall in subsect. 3.1 the algorithm described in [36]. Finally, in subsect. 3.2 we show the modifications that will allow us to construct invariant tori in a wider neighbourhood of the equilibrium point.

3.1 Kolmogorov normalization near an elliptic

equilibrium point

The goal is to introduce a suitable sequence of canonical transformations leading the Hamiltonian (14) in Kolmogorov normal form111By a little abuse of notation we denote by ${\underline{\omega}}$ the angular velocity vector, as usual in KAM theory. This should not be confused with the longitudes of perihelia, as usual in Celestial Mechanics., i.e.,

[TABLE]

where $({\underline{p}},{\underline{q}})\in\mathbb{R}^{n}\times\mathbb{T}^{n}$ are action-angle coordinates; thus, the surface ${\underline{p}}={\underline{0}}$ is invariant with respect to the flow induced by $H^{(\infty)}$ and the motion over that torus has angular velocities equal to the entries of a prescribed vector ${\underline{\omega}}\in\mathbb{R}^{n}$. The algorithm consists of a sequence of canonical transformations that we describe in three separated steps.

(i) Birkhoff normalization up to a finite degree

We first determine a generating function $B^{({\rm II})}$ by solving the equation

[TABLE]

The expansion of the transformed Hamiltonian $H^{({\rm II})}=\exp\mathcal{L}_{B^{({\rm II})}}\>H^{({\rm I})}$ can be written as

[TABLE]

where the occurrence of the new normal form term $f_{2}^{({\rm II})}=\langle f_{2}^{({\rm I})}\rangle_{{\underline{\varphi}}}$ has been highlighted, by separating it from the series of the perturbing terms. The recursive expression of $f_{2s}^{({\rm II})}$ as a function of $B^{({\rm II})}$ and $f_{2l}^{({\rm I})}$ (with $l\leq s$) can be computed by just collecting the homogeneous polynomials having the same degree in the square root of the actions. Thus, it is easy to check that the functions $f_{2s}^{({\rm II})}$ are of the same type as in (15).

We stress that the Birkhoff normal form up to degree $2$ in the actions is enough to start the following construction of the Kolmogorov normal form. On the other hand, by performing the Birkhoff normalization up to a degree higher than $3$, we can improve the numerical stability of the calculation of the coefficients appearing in the expansions generated by the algorithm. Here, we have computed the Birkhoff normalization up to the third degree in ${\underline{I}}$ which is good enough for our purposes, in the framework of the model we are studying. Therefore, the final Hamiltonian is given by $H^{({\rm III})}=\exp\mathcal{L}_{B^{({\rm III})}}\>H^{({\rm II})}$, that is

[TABLE]

where (a) the functions $f_{2s}^{({\rm III})}$ are homogeneous polynomials of degree $2s+2$ in the square root of the actions $I$ and are of type (15); (b) the generating function $B^{({\rm III})}$ is defined by the equation

[TABLE]

(c) $f_{4}^{({\rm III})}=\langle f_{4}^{({\rm II})}\rangle_{{\underline{\varphi}}}\,$. As a whole, the canonical transformation $\mathcal{B}$ inducing the Birkhoff normalization up to degree three in actions is explicitly given by $\mathcal{B}({\underline{I}},{\underline{\varphi}})=\exp\mathcal{L}_{B^{({\rm III})}}\,\circ\,\exp\mathcal{L}_{B^{({\rm II})}}\>({\underline{I}},{\underline{\varphi}})\,$, indeed it is easy to check that H^{({\rm III})}({\underline{I}},{\underline{\varphi}})=H^{({\rm I})}\big{(}\mathcal{B}({\underline{I}},{\underline{\varphi}})\big{)} using the exchange theorem for Lie series.

(ii) Initial translation of the actions

The canonical transformation $({\underline{I}},{\underline{\varphi}})={\mathfrak{T}}_{{\underline{I}}^{*}}({\underline{p}},{\underline{q}})$ performing the initial translation of the actions is of type

[TABLE]

Let us recall that we are constructing an invariant torus with a fixed frequency vector ${\underline{\omega}}$. Following [36], the initial translation can be determined in such a way that, in the integrable approximation, the quasi-periodic motions on the invariant torus $({\underline{p}}={\underline{0}}\,,\,{\underline{q}}\in\mathbb{T}^{n})$ have angular frequencies ${\underline{\omega}}$. Therefore, we determine the vector ${\underline{I}}^{*}$ with positive components (recall the definition of the canonical transformation $({\underline{x}},{\underline{y}})=\mathcal{A}({\underline{I}},{\underline{\varphi}})$) as the nearest to the origin solution of the equations

[TABLE]

We can write the expansion of H^{({\rm IV})}({\underline{p}},{\underline{q}})=H^{({\rm III})}\big{(}{\mathfrak{T}}_{{\underline{I}}^{*}}({\underline{p}},{\underline{q}})\big{)} as

[TABLE]

where, for $l\geq 0$ and $s\geq 0\,$, $f_{l}^{({\rm IV},s)}$ is a homogeneous polynomial of degree $l$ in the actions ${\underline{p}}$ and a trigonometric polynomial of degree either $2s$ or $2s-1$ in the angles ${\underline{q}}\,$. For short, let us introduce the symbol $\mathcal{P}_{l,2s}\,$, which denotes the set of functions that are homogeneous polynomials of degree $l$ in the actions and trigonometric polynomials of degree at most $2s$ in the angles, thus, $f_{l}^{({\rm IV},s)}\in\mathcal{P}_{l,2s}\,$. Moreover, using the Cauchy inequalities, one easily sees that the size (of any suitable norm) of $f_{l}^{({\rm IV},s)}$ can be estimated with an upper bound that is essentially proportional to the $s$-th power of the ratio of $\|{\underline{I}}^{*}\|$ over the analytic radius of convergence of $H^{({\rm III})}$, and it is inversely proportional to the $l$-th power of the minimum component of vector ${\underline{I}}^{*}$. Therefore, ${\underline{I}}^{*}$ plays a major role in the convergence of the expansions, because it is proportional to what is commonly identified as the small parameter of the KAM theory and it rules the radius of convergence for the actions ${\underline{p}}$. At this point, we want emphasize that we have some freedom in the crucial choice of the initial translation vector ${\underline{I}}^{*}$, as it will be discussed in sect. 3.2.

(iii) The standard Kolmogorov normalization algorithm

Let us describe the generic $r$-th step of the Kolmogorov normalization algorithm. We begin with a Hamiltonian of the type

[TABLE]

where $f_{l}^{(r-1,s)}\in\mathcal{P}_{l,2s}\,$, for $l\geq 0$ and $s\geq 0\,$. To fix ideas, we can start with $r=2$ defining $H^{({1})}=H^{({\rm IV})}$. Since we point to a Hamiltonian of type (16), we must remove the main perturbing terms of degree [math] and $1$ in the actions. We will proceed in two separate steps. We first remove part of the unwanted terms via a canonical transformation with generating function $\chi_{1}^{(r)}({\underline{q}})=X^{(r)}({\underline{q}})+{\underline{\xi}}^{(r)}\cdot{\underline{q}}$ (being ${\underline{\xi}}^{(r)}\in\mathbb{R}^{n}$). Thus, we solve with respect to $X^{(r)}({\underline{q}})$ and ${\underline{\xi}}^{(r)}\,$ the equations

[TABLE]

where the $n\times n$ matrix $C^{(r)}$ is defined by the equation $\frac{1}{2}C^{(r)}{\underline{p}}\cdot{\underline{p}}=f_{2}^{(r-1,0)}({\underline{p}})\,$. A unique solution satisfying $\langle X^{(r)}\rangle_{{\underline{q}}}=0$ exists if the frequencies ${\underline{\omega}}$ are non-resonant up to order $2r$, ${\underline{k}}\cdot{\underline{\omega}}\neq 0\,$, with ${\underline{k}}\in\mathbb{Z}^{n}$ such that $0<|{\underline{k}}|\leq 2r$, and if $\det C^{(r)}\neq 0\,$. We must now give the expressions of the functions $\hat{f}_{l}^{(r,s)}$ appearing in the expansion of the new Hamiltonian

[TABLE]

where $\hat{H}^{(r)}=\exp\mathcal{L}_{\chi_{1}^{(r)}}H^{(r-1)}$. To this aim, we will redefine many times the same quantity without changing the symbol. This is made by mimicking the C programming language, by using the notation $a\hookleftarrow b$ which means that the previously defined quantity $a$ is redefined as $a=a+b\,$. Therefore, we initially define

[TABLE]

To take into account the Poisson bracket of the generating function with ${\underline{\omega}}\cdot{\underline{p}}$, we put

[TABLE]

Then, we consider the contribution of the terms generated by the Lie series applied to each function $f_{l}^{(r-1,s)}$ as

[TABLE]

Looking at formulæ (27)–(29), one can easily check that $\hat{f}_{l}^{(r,s)}\in\mathcal{P}_{l,2s}$, for $l\geq 0$ and $s\geq 0\,$. We perform now a “reordering of terms”, by moving the monomials appearing in the expansion of a function $\hat{f}_{l}^{(r,s)}$ to another, in such a way that the so redefined functions $\hat{f}_{l}^{(r,s)}$ are homogeneous polynomials of degree $l$ in the actions and trigonometric polynomials of degree $2s$ or $2s-1$ in the angles, for $l\geq 0$ and $s\geq 0\,$.

In the second part of the $r$-th step of the Kolmogorov’s normalization algorithm, by using another canonical transformation, we remove the part of the perturbation up to the order of magnitude $r$ that actually depends on the angles and it is linear in the actions. Thus, we solve with respect to $\chi_{2}^{(r)}({\underline{p}},{\underline{q}})$ the equation

[TABLE]

where again the solution exists and it is unique if $\langle\chi_{2}^{(r)}\rangle_{{\underline{q}}}=0$ and the frequencies ${\underline{\omega}}$ are non-resonant up to order $2r$. Analogously to what we have done above, we now provide the expressions of the functions $f_{l}^{(r,s)}$ appearing in the expansion of the new Hamiltonian

[TABLE]

where $H^{(r)}=\exp\mathcal{L}_{\chi_{2}^{(r)}}\hat{H}^{(r)}$. We initially define

[TABLE]

In order to take into account the contribution of the terms generated by the Lie series applied to ${\underline{\omega}}\cdot{\underline{p}}\,$, we put

[TABLE]

Then, the contribution of the Lie series applied to the rest of the Hamiltonian $\hat{H}^{(r)}$ implies that

[TABLE]

Finally, we perform a new “reordering of terms”, so that at the end the functions $f_{l}^{(r,s)}\in\mathcal{P}_{l,2s}$ appearing in the expansion (31) of the new Hamiltonian $H^{(r)}$ are again homogeneous polynomials of degree $l$ in the actions and trigonometric polynomials of degree $2s$ or $2s-1$ in the angles, for $l\geq 0$ and $s\geq 0\,$.

Let us recall that the canonical transformation $\mathcal{K}^{(r)}$ inducing the Kolmogorov normalization up to the step $r$ is explicitly given by

[TABLE]

This conclude the $r$-th step of the algorithm that can be further iterated.

3.2 The modified algorithm constructing the Kolmogorov

normal form

As the prediction of the translation vectors ${\underline{\xi}}^{(r)}$ given by (25) can be affected by large errors, we have split the standard Kolmogorov normalization algorithm in two separate steps. First, we iterate for a fixed number of steps the normalization algorithm by setting to zero the translation vectors ${\underline{\xi}}^{(r)}$. This procedure is reminiscent of Arnold’s proof of the KAM theorem ([1] and see also [47] and [13], where such a modification of Kolmogorov normalization algorithm has been recently adapted so as to approximate elliptic lower dimensional tori and to prove their existence). Therefore, under mild theoretical assumptions, such a partial normalization procedure can still converge to a Hamiltonian in Kolmogorov normal form related to a vector of angular frequencies ${\underline{\omega}}^{*}$ different from ${\underline{\omega}}$. In practice, we perform just a finite number of steps of this partial normalization. Then, using this intermediate Hamiltonian H\simeq{\underline{\omega}}^{*}\cdot{\underline{p}}+\mathcal{O}\big{(}\|{\underline{p}}\|^{2}\big{)} as the initial one, we restart a complete standard Kolmogorov normalization algorithm now including the translation vectors ${\underline{\xi}}^{(r)}$ defined in (25). This splitting of the normalization algorithm in two separate steps becomes advantageous if, after the first step, the frequencies ${\underline{\omega}}^{*}$ are sufficiently close to ${\underline{\omega}}$ and, thus, the translation vectors ${\underline{\xi}}^{(r)}$ are small enough. Let us now remark that the frequencies ${\underline{\omega}}^{*}$, related to the intermediate Hamiltonian $H\simeq{\underline{\omega}}^{*}\cdot{\underline{p}}+\mathcal{O}(\|{\underline{p}}\|^{2})$, depend on the initial translation vector ${\underline{I}}^{*}$ of the canonical transformation (21). Therefore, we try to choose ${\underline{I}}^{*}$ in such a way that the frequency vector ${\underline{\omega}}^{*}$ is as close as possible to ${\underline{\omega}}$.

We now provide a detailed description of our algorithm, that has been formulated so as to reproduce the main ideas explained in the previous heuristic discussion. Let us start from a system of the type described by the Hamiltonian $H^{({\rm I})}$ in (14). Then, let us carry out the following steps.

1. (a)

Perform the canonical transformation $\mathcal{B}$ that realizes the Birkhoff normalization up to a finite degree, as described at point (i) of sect. 3.1. 2. (b)

Determine a good initial translation vector ${\underline{\hat{I}}}$ by proceeding as follows.

1. (b1)

Let us refer to the initial conditions as $({\underline{I}}_{0},{\underline{\varphi}}_{0})$; calculate their values in the new coordinates, say $({\underline{I}}_{0}^{*},{\underline{\varphi}}_{0}^{*})=\mathcal{B}({\underline{I}}_{0},{\underline{\varphi}}_{0})$. Then, perform the initial translation of the actions, as at point (ii) of sect. 3.1, by replacing ${\underline{I}}^{*}$ with ${\underline{I}}_{0}^{*}\,$. 2. (b2)

Let us perform the Kolmogorov normalization algorithm up to a fixed $R^{\prime}$-th step, as at point (iii) of sect. 3.1, starting from $H^{({1})}=H^{({\rm IV})}$, but putting the translation vectors ${\underline{\xi}}^{(r)}={\underline{0}}\,$, for $r=2,\,\ldots,\,R^{\prime}\,$. Let us define ${\underline{\omega}}_{0}^{*}$ in such a way that ${\underline{\omega}}_{0}^{*}\cdot{\underline{p}}={\underline{\omega}}\cdot{\underline{p}}+f_{1}^{(0,r+1)}({\underline{p}})$, with $f_{1}^{(0,r+1)}$ as obtained at the end of such procedure. $R^{\prime}$ is a fixed integer parameter that is selected sufficiently large to allow the convergence of the whole algorithm, but also taking into consideration the computational resources available. 3. (b3)

Let us improve our choice of ${\underline{I}}_{0}^{*}$, by approximating numerically (by the finite differences method) the Jacobian matrix $\mathcal{J}_{{\underline{I}}_{0}^{*}}$ of the function ${\underline{\omega}}^{*}({\underline{I}}_{0}^{*})$ and then solving the linear equation $\mathcal{J}_{I_{0}^{*}}({\underline{\hat{I}}}-{\underline{I}}^{*})={\underline{\omega}}-{\underline{\omega}}_{0}^{*}$ in the unknown ${\underline{\hat{I}}}$. 3. (c)

Perform the translation of the actions as at point (ii), sect. 3.1, replacing ${\underline{I}}^{*}$ with ${\underline{\hat{I}}}$. 4. (d)

Let us perform again the Kolmogorov normalization algorithm without translations, as at step (${\rm b}_{2}$). In what follows we will denote as $\{{\mathfrak{H}}^{(r)}\}_{r=1}^{R^{\prime}}\,$, $\{{\mathfrak{X}}_{1}^{(r)},{\mathfrak{X}}_{2}^{(r)}\}_{r=2}^{R^{\prime}}$ and $\{{\mathfrak{K}}^{(r)}\}_{r=2}^{R^{\prime}}$ the obtained finite sequences of the Hamiltonians, the generating functions and the canonical transformations, respectively; so that ${\mathfrak{H}}^{({1})}=H^{({\rm IV})}$, {\mathfrak{H}}^{(r)}=\exp\mathcal{L}_{{\mathfrak{X}}_{2}^{(r)}}\big{(}\exp\mathcal{L}_{{\mathfrak{X}}_{1}^{(r)}}{\mathfrak{H}}^{(r-1)}\big{)} and ${\mathfrak{H}}^{({r})}={\mathfrak{H}}^{({1})}\circ{\mathfrak{K}}^{(r)}$, for $r=2,\,\ldots,\,R^{\prime}\,$. Therefore, ${\mathfrak{K}}^{(r)}$ is explicitly given by a formula analogous to (35), by replacing the symbols $\mathcal{K}$ and $\chi$ with ${\mathfrak{K}}$ and ${\mathfrak{X}}$, respectively. 5. (e)

Let us perform the standard Kolmogorov normalization algorithm, as at point (iii) of sect. 3.1 (with the translation vectors ${\underline{\xi}}^{(r)}$ given by (25)), starting from $H^{({1})}={\mathfrak{H}}^{(R^{\prime})}$.

Let us remark that steps (b1)–(b3) of the previous procedure can be iterated more than once, in order to refine the calculation of the initial translation vector ${\underline{\hat{I}}}\,$; however, from a practical point of view, a single execution of the steps (b1)–(b3) is often enough to successfully perform the subsequent Kolmogorov normalization algorithms, in the sense that one can clearly appreciate their numerical convergence. In particular, this applies also in the case of our model of the planetary problem we are studying.

The knowledge of the normal form (and of the normalizing transformations) allows to explicitly implement a semi-analytic procedure producing the integration of the equations of motion. For instance, if we consider the coordinates $({\underline{x}},{\underline{y}})$ introduced at the very beginning of the present section (where the canonical transformation $\mathcal{D}$ was defined) and the action-angle variables $({\underline{p}},{\underline{q}})$ of the normal form, we have that

[TABLE]

with $\mathcal{C}^{(\infty)}=\lim_{r\to\infty}\mathcal{C}^{(r)}$ and $\mathcal{C}^{(r)}=\mathcal{A}\circ\mathcal{B}\circ{\mathfrak{T}}_{{\underline{\hat{I}}}}\circ{\mathfrak{K}}^{(R^{\prime})}\circ\mathcal{K}^{(r)}$, where the canonical transformation $\mathcal{A}$ is defined at the very beginning of the present section, while $\mathcal{B}\,$, ${\mathfrak{T}}_{{\underline{\hat{I}}}}\,$, ${\mathfrak{K}}^{(R^{\prime})}$ and $\mathcal{K}^{(r)}$ are determined at points (a), (c), (d) and (e) of the algorithm described above, respectively.

3.3 Application to the secular model

We have checked our implementation by comparing the numerical integration of the flow induced by the Hamiltonian $H^{(\mathcal{D})}$ (written in (13)) against the results from the normal form, by means of the scheme (36). The input data needed to start the execution of the algorithm constructing the Kolmogorov normal form are the coefficients appearing in the expansion of $H^{(\mathcal{D})}$ and the frequency vector ${\underline{\omega}}$ identifying the KAM torus we are looking for. In order to determine ${\underline{\omega}}$, we preliminary perform a long-time numerical integration of the motion law t\mapsto\big{(}{\underline{x}}(t)\,,\,{\underline{y}}(t)\big{)} starting from the initial conditions reported in Table 2. Then, we apply the frequency analysis method (see [28] and references therein) to the signals $x_{j}(t)+{\rm i}y_{j}(t)\,$, with $j=1,\,2,\,3\,$. The corresponding fundamental frequencies detected by that numerical method provide us the wanted values

[TABLE]

Of course, formula (36) refers to an ideal scheme, an infinite sequence of canonical transformations cannot be performed on a computer; thus, we limited ourselves to approximate the exact solution by replacing $\mathcal{C}^{(\infty)}$ with $\mathcal{C}^{(25)}$ and putting also $R^{\prime}=25$. Both the (approximate) normalizing transformations $\mathcal{A}\circ\mathcal{B}\circ{\mathfrak{T}}_{{\underline{\hat{I}}}}\circ{\mathfrak{K}}^{(25)}\circ\mathcal{K}^{(25)}$ and its inverse have been computed by using our software for algebraic manipulations. The results of a long-time comparison (lasting for $10^{7}$ years) are summarized in Figure 1: the relative error on the determination of the two-dimensional vectors \big{(}x_{j}(t)\,,\,y_{j}(t)\big{)} is always smaller than $0.003\,$%, for $j=1,\,2,\,3$. We stress that the accuracy of the semi-analytic solution is mainly affected by the unavoidable truncation rules that are applied to the expansions. Let us shortly describe them all together as follows: we have put the same limits on the partial Birkhoff normalization as those holding for $\mathcal{H}^{({\rm sec})}$, i.e., the Hamiltonians $H^{({\rm I})},\,\ldots\,\,H^{({\rm IV})}$ have been truncated up to degree 18 in the square root of the actions ${\underline{I}}\,$; during both the Kolmogorov normalization algorithms the expansions of the Hamiltonians were limited to the fifth degree in actions ${\underline{p}}$ and up to a maximal trigonometric degree in ${\underline{q}}$ equal to 50. The truncation rules applied to the normalizing transformations are analogous to those for the Hamiltonians.

Taking into account all the factors limiting the accuracy of our approximation of the semi-analytic integration scheme (36), we consider that the agreement with the numerical integration is excellent. In our opinion, this is a strong validation of all the procedure constructing the Kolmogorov normal form. We also think that, with a relevant additional effort, our approach could be translated in a computer-assisted proof of existence of such an invariant torus, in a similar way to what was done in [36] and [9]. However, this goes far beyond the scopes of the present work.

4 Effective stability in the neighborhood of the KAM torus

We follow the approach described in [39] and already applied to another planetary problem in [12]. We focus on the Kolmogorov normal form (16), which is related to the locally unique invariant surface, that is travelled by quasi-periodic motion characterized by the angular velocity vector ${\underline{\omega}}$. $H^{(\infty)}$ is now reconsidered as a new starting point for performing a Birkhoff normalization. Such a further constructive procedure finally aims to provide bounds on the diffusion from that torus having angular velocity equal to ${\underline{\omega}}$; this is made by giving a lower estimate on the escape time $T(\varrho)$. Here, $\varrho$ is the distance (as measured in actions) from the initial torus related to ${\underline{\omega}}$ and plays the role of the small parameter in the construction of the Birkhoff normal form.

4.1 Formal algorithm constructing the Birkhoff normal form close to the KAM torus

In order to properly describe the new normalization algorithm, it is convenient to introduce the new symbol $\mathcal{H}^{(2,0)}=H^{(\infty)}$ for the (re)starting Hamiltonian and we explicitly write its expansion as

[TABLE]

being $f_{l}^{(2,0;s)}\in\mathcal{P}_{l,2s}\,$. Unfortunately, in order to obtain final estimates that are good enough, we cannot simply prescribe that, as a rule to perform the $l$-step, all the perturbative terms having the same $l$ degree in actions must be removed at the same time; this would somehow natural but too rough. Therefore, we need to design a normalization procedure based on two indexes. Let us describe the generic $(r,{\mathfrak{r}})$ step, which aims to introduce the new Hamiltonian $\mathcal{H}^{(r,{\mathfrak{r}})}=\exp\mathcal{L}_{{\mathfrak{B}}^{(r,{\mathfrak{r}})}}\>\mathcal{H}^{(r,{\mathfrak{r}}-1)}$, where the expansion of the previous Hamiltonian reads

[TABLE]

being, again, $f_{l}^{(r,{\mathfrak{r}}-1;s)}\in\mathcal{P}_{l,2s}\,$. The generating function is determined so as to remove the main perturbing term among those having $r$ degree with respect to the actions ${\underline{p}}$; in view of the Fourier decay, this is $f_{r}^{(r,{\mathfrak{r}}-1;r)}$, thus we set the homological equation

[TABLE]

that must be obviously solved with respect to ${\mathfrak{B}}^{(r,{\mathfrak{r}})}\in\mathcal{P}_{r,2{\mathfrak{r}}}\,$. We proceed in an analogous way to what has been done in subsect. 3.1, in order to define all the terms appearing in the new expansion

[TABLE]

First, we initially set

[TABLE]

Accordingly to the homological equation (40), we redefine the normal forms term of $r$ degree in ${\underline{p}}$ as

[TABLE]

In order to take into account the contribution of all the terms (but the first two) that are generated by the Lie series applied to ${\underline{\omega}}\cdot{\underline{p}}$, we put

[TABLE]

Then, the contribution of the Lie series applied to the rest of the Hamiltonian $\mathcal{H}^{(r,{\mathfrak{r}}-1)}$ implies that

[TABLE]

where some of the new summands actually have not any influence, because $f_{l}^{(r,{\mathfrak{r}}-1;s)}=0$ if $2\leq l\leq r$ and $s>0$ or $l=r$ and $1\leq s<{\mathfrak{r}}$ (in agreement with equation (39)). All the prescriptions (42)–(45), have been settled so as to preserve the rules about the classes of functions, so that $f_{l}^{(r,{\mathfrak{r}};s)}\in\mathcal{P}_{l,2s}\,$, as it can be easily checked by induction. Since the structure of the expansion (41) is coherent with that in (39), the procedure can be iterated at the next $(r,{\mathfrak{r}}+1)$ step.

In order to start the elimination of the perturbing terms of higher degree in ${\underline{p}}\,$, we simply put $\mathcal{H}^{(r+1,0)}=\lim_{{\mathfrak{r}}\to\infty}\mathcal{H}^{(r,{\mathfrak{r}})}$. Since the expansion (38) of the initial Hamiltonian $\mathcal{H}^{(2,0)}$ is identical to that in (39) when $r=2$ and ${\mathfrak{r}}=1$, then the sequence of the generating functions in our normalization scheme can be ideally represented as follows: ${\mathfrak{B}}^{(2,1)},\>{\mathfrak{B}}^{(2,2)},\>\ldots\>{\mathfrak{B}}^{(2,\infty)},{\mathfrak{B}}^{(3,1)},\>\ldots\>{\mathfrak{B}}^{(3,\infty)},{\mathfrak{B}}^{(4,1)},\>\ldots\>$

4.2 Iterative scheme of estimates

Although explicit computations according to the previous formal algorithm are feasible (see [12]), the expansions of the functions introduced by that procedure become quickly so cumbersome, that it is hard to deal with them for any algebraic manipulator when the truncation orders are increased with respect to both the actions and the angles. Here, we limit ourselves to produce a rather rough lower bound on the diffusion time, in order to drastically reduce the computational difficulty. This is mainly made by iterating a scheme of estimates involving just the norm for each of those functions, instead of computing all their expansions. In practice, for any generic function $g\in\mathcal{P}_{l,2s}\,$, we define its norm as

[TABLE]

being $\{c_{i_{1},\ldots,i_{n},k_{1},\ldots,k_{n}}\}_{{\underline{i}},{\underline{k}}}$ the finite set of coefficients appearing in the corresponding Taylor-Fourier series

[TABLE]

where the notation ${\scriptstyle{{\displaystyle{\sin}}\atop{\displaystyle{\cos}}}}$ means that either the sine or the cosine can appear and, for any fixed harmonic ${\underline{k}}$, the choice is made unique according to the following usual criterion: if there is an index $1\leq j\leq n$ such that $k_{j}<0$ and $k_{1}=\ldots=k_{j-1}=0$ then the sine function appears, otherwise the cosine.

Our iterative scheme of estimates mainly aims to provide a set of computational rules to determine a sequence of majorants $\mathcal{F}_{l}^{(r,{\mathfrak{r}};s)}$ such that

[TABLE]

Let us recall that the induction must be started from $\mathcal{H}^{(2,0)}=H^{(\infty)}$. To fix the ideas, let us refer to the explicit calculation of the Kolmogorov normal form and its truncation rules, as they have been described in subsect. 3.3. Therefore, $\mathcal{F}_{l}^{(2,0;s)}$ can be computed for $2\leq l\leq 5$ and $0\leq s\leq 25$, by simply using the definition (46) of the norm, while for all the remaining initial majorants we put $\mathcal{F}_{l}^{(2,0;s)}=0$ when $l>5$ or $s\geq 25$.

We now describe how the iteration of the estimates works in the case of the generic $(r,{\mathfrak{r}})$ constructive step of the Birkhoff normal form. By induction hypothesis, let us suppose to know the upper bounds related to the previous step, that are $\|f_{l}^{(r,{\mathfrak{r}}-1;s)}\|\leq\mathcal{F}_{l}^{(r,{\mathfrak{r}}-1;s)}$ for $l\geq 2,\ s\geq 0$. First, we estimate the generating function so that

[TABLE]

being222From a practical point of view, for not too large values of ${\mathfrak{r}}$, a lower bound on $\min_{0<|{\underline{k}}|\leq 2{\mathfrak{r}}}|{\underline{k}}\cdot{\underline{\omega}}|$ can be rigorously computed by using interval arithmetics, while the asymptotic behavior can be estimated by using the Diophantine inequality $|{\underline{k}}\cdot{\underline{\omega}}|\geq\gamma/|{\underline{k}}|^{\tau}$, for suitable $\gamma>0$ and $\tau\geq n-1$. For instance, [8] provides a general prescription rule to determine a vector satisfying the Diophantine property with the optimal value for $\tau=n-1$. $\alpha_{\mathfrak{r}}=\min_{0<|{\underline{k}}|\leq 2{\mathfrak{r}}}|{\underline{k}}\cdot{\underline{\omega}}|$. In order to take into account the prescriptions (42)–(43), we initially put

[TABLE]

and we set the following re-definitions

[TABLE]

The new contributions to the perturbing terms are easily evaluated, after having verified the following inequality:

[TABLE]

that is easy to check by induction on $j$. Indeed, in the case $j=1$ one has simply to evaluate the coefficients generated by the derivatives appearing in the Poisson brackets, while noticing that $\mathcal{L}_{{\mathfrak{B}}^{(r,{\mathfrak{r}})}}^{j-1}f_{l}^{(r,{\mathfrak{r}}-1;s)}\in\mathcal{P}_{l+(j-1)(r-1),2(s+(j-1){\mathfrak{r}})}$ is essential to deal with the generic case. Thus, the prescriptions (44)–(45) can be translated into the following rules for the majorants:

[TABLE]

This ends the description of the $(r,{\mathfrak{r}})$ iterative step making part of our scheme of estimates.

Let us restart from our truncated expansions of $\mathcal{H}^{(2,0)}=H^{(\infty)}$, which have been determined as described in subsect. 3.3; the subsequent computation of the majorants $\mathcal{F}_{l}^{(r,{\mathfrak{r}};s)}$ is easy to code and require a small amount of CPU time. In particular, we can provide a simple rule to pass from the generating functions of a fixed degree in the actions ${\underline{p}}$ to the next one. Actually, we define $\mathcal{H}^{(r+1,0)}=\mathcal{H}^{(r,25(r-1))}$ and, correspondingly,

[TABLE]

This is due to the fact that ${\mathfrak{B}}^{(r,{\mathfrak{r}})}=0$ for ${\mathfrak{r}}>25(r-1)$, because of the truncation up to trigonometric degree $50$ of the computation of the Kolmogorov normal form.

In particular, it is immediate to provide an estimate of the remainder terms appearing in the Birkhoff normal form of degree $r$ in the actions ${\underline{p}}$, that is

[TABLE]

being $\mathcal{R}_{l}^{(r+1)}=\sum_{s\geq 0}f_{l}^{(r+1,0;s)}$ for $l>r$. In fact, the following inequality holds true:

[TABLE]

4.3 Evaluation of the stability time for the secular model

In this subsection, we adapt the approach developed in [12] to the present context. It is convenient to consider the domain

[TABLE]

Therefore, for any generic function $g\in\mathcal{P}_{l,2s}\,$, the inequality

[TABLE]

holds true as an immediate consequence of the definition (46) of $\|\cdot\|\,$. Let us consider a particular initial condition \big{(}{\underline{p}}(0),{\underline{q}}(0)\big{)}\in\Delta_{\varrho_{0}}\times\mathbb{T}^{n} with $\varrho_{0}<\varrho$; let $T_{e}$ be the minimum value (for all ${\underline{q}}(0)\in\mathbb{T}^{n}$) for which the corresponding motion laws are such that ${\underline{p}}(t)\in\Delta_{\varrho}$ when $|t|<T_{e}$. Let us refer to $T_{e}$ as the escape time from the domain $\Delta_{\varrho}\,$. This is the quantity that we want to evaluate. To this end we use the elementary estimate

[TABLE]

The latter formula allows us to provide a lower estimate for the escape time from the domain $\Delta_{\varrho}\,$, namely

[TABLE]

where we have kept just the first term of the series $\sum_{l>r}\|\{p_{j},\mathcal{R}_{l}^{(r+1)}\}\varrho^{l+1}\|$ because we restrict ourselves to consider values of $\varrho$ that are safely within its analyticity domain. Moreover, the estimate of $\|\{p_{j},\mathcal{R}_{r+1}^{(r+1)}\}\|$ is similar to (55) for $\mathcal{R}_{l}^{(r+1)}$.

In a practical application, usually $\varrho_{0}$ is fixed by the initial data, while $\varrho$ and $r$ are left arbitrary. Thus, we proceed by looking for the maximal value of $\tau(\varrho_{0},\varrho,r)$ with respect to $\varrho$ and $r$. First we keep $r$ fixed and optimize the function $\varrho\mapsto\tau(\varrho_{0},\varrho,r)$; this allows us to obtain $\varrho=\frac{r+1}{r}\varrho_{0}$ and to introduce the new function

[TABLE]

Next we look for the optimal value $r_{\rm opt}$ of $r$, which maximizes $\tilde{\tau}(\varrho_{0},r)$ when $r$ is allowed to change. This means that we look for the quantity

[TABLE]

which is our best estimate of the escape time, depending only on the initial data. We define the latter quantity as the estimated stability time. After having compared the estimate in (55) with the denominator appearing in (59), it is natural to expect that the function $r\mapsto\tilde{\tau}(\varrho_{0},r)$ asymptotically behaves in the same way as the inverse of the remainder of the Birkhoff normal form, that is $\tau(\varrho_{0},r)\sim C^{r}/[\varrho^{r}(r!)^{n}]$, where $C$ is a suitable constant and we assumed the optimal Diophantine condition for the quasi-periodic motion on the initial KAM torus, i.e. $|{\underline{k}}\cdot{\underline{\omega}}|>\gamma/|{\underline{k}}|^{n-1}$, for a fixed value of $\gamma>0$. Let us recall that this kind of classical estimates is essential to prove that the stability time grows exponentially with respect to the inverse of the distance $\varrho_{0}$ from that torus. On the other hand, for any fixed value of $\varrho_{0}\,$, $\tau(\varrho_{0},r)$ quickly reduces to zero, then there is a finite optimal value $r_{\rm opt}\,$, maximizing $\tilde{\tau}(\varrho_{0},r)$. In the two plots of Figure 2, we reported both the step-wise function $r_{\rm opt}(\varrho)$ and the corresponding estimated stability time $T(\varrho_{0})=\tilde{\tau}(\varrho_{0},r_{\rm opt})$, according to (60). Actually, Figure 2 collects the explicit results for our secular model. This has been possible, because we preliminary computed the upper bounds $\mathcal{F}_{r+1}^{(r+1,0;s)}$ for $0\leq s\leq 25\,r\,$, $r\leq 20\,$. In other words, starting from our truncated expansions of $\mathcal{H}^{(2,0)}=H^{(\infty)}$, we have explicitly calculated the majorants for the first 19 normalization steps constructing the Birkhoff normal form around the initial torus, whose corresponding angular velocity vector is ${\underline{\omega}}$.

In our framework, the best proof for the effective stability of the system would be produced by verifying that the set of possible initial conditions (taking into account the limitations on the knowledge of these values due to the observations) is covered by the domain $\Delta_{\varrho_{0}}\times\mathbb{T}^{n}$. In what follows, we make such a comparison.

For what concerns the set of possible initial conditions, we can refer to the values of $\Delta\xi_{j}$ and $\Delta\eta_{j}$ that are reported for $j=1,\,2$ (corresponding to Jupiter and Saturn, respectively) in Table 3 of [12]: the uncertainties on the determination of the secular actions are given by

[TABLE]

where the secular coordinates $\xi_{j}$ and $\eta_{j}$ can be calculated by using formula (2) and the initial conditions listed in Table 1. Moreover, in (61) we have assumed that the effects induced by the diagonalizing canonical transformation $({\underline{\xi}},{\underline{\eta}})=\mathcal{D}({\underline{x}},{\underline{y}})$ are negligible; this is actually confirmed by visual inspection of the coefficients appearing in the matrix associated with the linear operator $\mathcal{D}$. All the subsequent changes of coordinates defined by the normalization procedure are either near-to-the-identity or rigid translations; thus, it is natural to assume that the values of the actions $\Delta p_{j}=\Delta I_{1}+\Delta I_{2}\simeq 7\,\times 10^{-7}$, for $j=1,\,2,\,3$, approximately gives us the radii of the ball containing all the possible initial conditions, which are coherent with the observations. Let us remark that the value of the secular action of Uranus (that is about $22$ times lighter than Jupiter) is one order of magnitude smaller than that of the major planet of our Solar System; this explains why $\Delta I_{3}$ has not been involved in the previous approximate evaluation of $\Delta p_{j}\,$.

For what concerns our estimates, Figure 2 clearly shows that our model is stable for a time comparable to the estimated age of the Universe (around $10^{10}$ years) in a neighborhood of the initial torus having a radius $\varrho_{0}$ slightly smaller than $10^{-7}$.

Therefore, our estimates about the domain of effective stability could look slightly disappointing: its radius is just one order of magnitude smaller with respect to that related to the possible initial conditions of our model. Nevertheless, we think that our approach could be improved so as to obtain a fully satisfying result. Indeed, the evaluation of the ball radius of the initial conditions is quite pessimistic, as it has been made by following [12] and, therefore, it is based on observations made a few decades ago; the technological progress is continuously reducing the uncertainties of the measures. Moreover, the radius $\varrho_{0}$ of the effective stability domain can be significantly enlarged, by performing explicit calculations for the expansions of the Birkhoff normal form in the neighborhood of the initial KAM torus. This can be preliminary done in such a way to improve the results produced by the iterative scheme of estimates333The effectiveness of an approach complementing explicit expansions of the Birkhoff normal form with rigorous computer-assisted estimates has been recently shown in the classical case-study concerning the Henon-Heiles system. These results are described in C. Caracciolo: “Studio rigoroso della stabilità effettiva di sistemi Hamiltoniani quasi-integrabili: stime computer-assisted”, thesis, Master in Mathematics, Univ. of Roma “Tor Vergata” (2016), that is available on request to the Author..

Acknowledgments

This work has been supported by the research program “Teorie geometriche e analitiche dei sistemi Hamiltoniani in dimensioni finite e infinite”, PRIN 2010JJ4KPA009, financed by MIUR.