Computer-assisted estimates for Birkhoff normal forms
aa r X i v : . [ m a t h - ph ] F e b Computer-assisted estimates for Birkhoffnormal forms ∗ CHIARA CARACCIOLO
Dipartimento di Matematica, Universit`a degli Studi di Roma “Tor Vergata”,via della Ricerca Scientifica 1, 00133 — Roma (Italy).
UGO LOCATELLI
Dipartimento di Matematica, Universit`a degli Studi di Roma “Tor Vergata”,via della Ricerca Scientifica 1, 00133 — Roma (Italy).e-mails: [email protected], [email protected],
Abstract
Birkhoff normal forms are commonly used in order to ensure the so called “ef-fective stability” in the neighborhood of elliptic equilibrium points for Hamiltoniansystems. From a theoretical point of view, this means that the eventual diffusioncan be bounded for time intervals that are exponentially large with respect to theinverse of the distance of the initial conditions from such equilibrium points. Here,we focus on an approach that is suitable for practical applications: we extend arather classical scheme of estimates for both the Birkhoff normal forms to any finiteorder and their remainders. This is made for providing explicit lower bounds of thestability time (that are valid for initial conditions in a fixed open ball), by using afully rigorous computer-assisted procedure. We apply our approach in two simplecontexts that are widely studied in Celestial Mechanics: the H´enon-Heiles modeland the Circular Planar Restricted Three-Body Problem. In the latter case, we ad-apt our scheme of estimates for covering also the case of resonant Birkhoff normalforms and, in some concrete models about the motion of the Trojan asteroids, weshow that it can be more advantageous with respect to the usual non-resonant ones. ∗ Primary: 37J40; Secondary: 37N05, 70–08, 70F07, 70H08.
Key words and phrases: computer-assisted proofs, normal form methods, effective stability, Hamiltoniansystems, Celestial Mechanics.
C. Caracciolo, U. Locatelli
The Birkhoff normal form represents a milestone in the historical development of theHamiltonian perturbation theory, because it has been designed in the simplest possibleway so as to (at least) partially escape the results obtained by Poincar´e for what concernsboth the non-existence of the first integrals and the ubiquity of chaos (see [40]). Indeed, itwas introduced about one century ago (see [48], [10]–[11] and [5]) in the framework of thedynamics in the neighborhood of an elliptic equilibrium point; in section 2 such a contextis recalled by a set of proper definitions. More recently, the scheme of analytical estimateshas been further improved in order to ensure that the stability time is exponentially bigwith respect to (some fractional power of) the inverse of the distance from the equilibriumpoint (see [34]–[35], [38]–[39] and section 4 of [20] for a brief summary of such a kind ofresults).The approach based on suitable estimates for Birkhoff normal forms has allowed tointroduce the concept of effective stability for physical systems, that can be stated asfollows. Let us focus on the motions starting from initial conditions staying in a ball ofradius ̺ , that is centered at an equilibrium point; we say that it is effectively stable , ifwe can ensure that the solutions of the Hamilton equations stay at a distance not largerthan ̺ > ̺ for an interval of time T that is exceeding the (expected) life-time of thatsystem. This kind of approach has been successfully applied to the so called CircularPlanar Restricted Three-Body Problem (hereafter, CPRTBP), where the vicinity of thetriangular Lagrangian points ( L and L ) in a model including Sun and Jupiter as primarybodies is considered. In [24], it is shown that the orbits of the third massless body staywithin a ball B ̺ ( ) having radius ̺ and centered in L or L for times larger than theage of the universe, where the value of ̺ is approximately the same as ̺ and the ball B ̺ ( ) is large enough to cover the initial conditions corresponding to the observations offour Trojan Asteroids. Such a result has been extended in several ways, with applicationsto problems arising, mainly, in Celestial Mechanics; the following list is far from beingexhaustive. In [46], a Trojan Asteroid is shown to be effectively stable in the 3D extensionof the model described just above. In [16], the effects induced by Saturn are taken intoaccount indirectly, because the problem is restricted in such a way that the orbit of Jupiteris prescribed so as to be both periodic and a good approximation of a numerical solutionfor a model including Sun, Jupiter and Saturn. In the models studied in [14] and [33],the continuous Hamiltonian flows have been replaced by symplectic mappings describinga CPRTBP and an Elliptic Restricted Three-Body Problem, respectively. In [44], theeffective stability of the rotational motion of Mercury is shown in the framework of aHamiltonian model. Last but not least, in [43] Birkhoff normal forms are used in orderto efficiently calculate counter-terms, with the purpose of controlling the diffusion in asymplectic map of H´enon type. This has been made also in view of further possibleextensions to similar models describing the dynamics of particle accelerators.In spite of its paramount importance, in none of the previously mentioned articlesthe construction of the Birkhoff normal form has been complemented with any rigorous computer-assisted proof. On the one hand, the results obtained in all these papers arebased on the explicit calculation of huge numbers of coefficients appearing in the Birkhoff omputer-assisted estimates for Birkhoff normal forms rigorous proofs. However, completelyrigorous computer-assisted estimates have been carried out in a paper that is focusedon CPRTBP and is an ancestor of a few of those mentioned above, but its results werestrongly limited by a too naive choice of the initial Hamiltonian expansions (compare [8]with [24]).Computer-assisted proofs have been used in several different fields of Mathematics,for instance, ranging from combinatorial theory (the celebrated Four-Colour theorem,see [1]–[3]) to PDEs (for a recent work, see, e.g., [6]). In particular, they have allowedto increase very remarkably the threshold of applicability with respect to the small para-meter, that is a crucial problem in KAM theory (see [32], [37] and [4]). In fact, whileat the very beginning the first purely analytical results were so limited that they werecompletely inadequate for any realistic application, recent rigorous computer-assisted es-timates have been able to get very close to the optimal breakdown threshold on the smallparameter, at least for some specific problems (see [28] and [15], respectively). Until now,the performances of approaches based on the application of KAM theorem are signific-antly better with respect to those using Birkhoff normal forms (for instance, compare [22]with [41], where a secular planetary model including Sun, Jupiter, Saturn and Uranusis considered). For what concerns the very specific case of the CPRTBP, [17] includesrigorous computer-assisted proofs of existence of KAM tori in the vicinity of three TrojanAsteroids. Moreover, in that same work it is shown the numerical evidence that the al-gorithm constructing the Kolmogorov normal form is successful in 23 cases over the first34 Trojan Asteroids listed in the IAU Catalogue. Let us recall that more than 7 000 aster-oids have been detected about the triangular Lagrangian points of the system having Sunand Jupiter as primary bodies. Such a large number of observed objects, compared withthe rather small one for which the previously described results are available, highlightsthe need for improvements of the mathematical framework, in order to fully explain theorbital stability of these celestial bodies.Our work aims to contribute in filling some of the gaps that have been discussedabove. In particular, we want to settle a completely rigorous scheme of computer-assistedestimates for Birkhoff normal forms. Moreover, we want to test its performances. In viewof such a purpose, we decided to focus just on a pair of Hamiltonian systems, that arefundamental examples and have two degrees of freedom. First, we apply our method tothe H´enon-Heiles model (see [27] and the very recent review in [12]), that is a very simpleexample of a perturbed couple of harmonic oscillators. In particular, we study the casewith two non-resonant frequencies having opposite signs, because this implies the failureof any trivial proof scheme of the stability. Indeed, the Hamiltonian is not a Lyapunovfunction because of the lack of convexity around the origin. As a second application,we reconsider the classical problem of stability of the triangular Lagrangian equilibria inthe CPRTBP. Let us recall that the structure of the Hamiltonian for the latter model is C. Caracciolo, U. Locatelli very similar to that of the former one; in particular, trivial approaches to the problemof stability are hopeless for both. For what concerns the CPRTBP, we study a varietyof situations by considering a few sub-cases where the role of the second body is playedby some different planets. This allows us to investigate the problem also with very smallmass ratios between the primaries and to check when an approach based on resonant
Birkhoff normal forms can be advantageous. Let us recall that the first results about theexpansions of the formal integrals of motion for resonant systems go back at the dawn ofthe computer age and refer to the H´enon-Heiles model (see, e.g., [26]).While in KAM theory there are powerful computer-assisted techniques that are notbased on the Hamiltonian formalism for canonical transformations (see [7] and [15]), weemphasize that the results about effective stability are intrinsically related with normalforms. In order to explicitly construct them, we use the composition of the Lie series,that is a very mature algorithm and can be complemented with a complete scheme ofanalytical estimates (see [18] for an introduction).The paper is organised as follows. In section 2, we describe the construction of theBirkhoff normal form by an approach based on Lie series. This is the starting point forthe definition of the scheme of estimates, which is explained in section 3. In section 4, wediscuss how to use the estimates to prove effective stability for Birkhoff normal forms inboth the non-resonant case and the resonant one. In section 5, we describe our results forthe H´enon-Heiles model and for the CPRTBP. Main conclusions are drawn in section 6.Appendix A is devoted to an introduction to those concepts of validated numerics, that areessential to implement our computer-assisted proofs in a fully rigorous way. In appendix B,we describe in detail an application to the H´enon-Heiles model in a situation, that istailored in such a way that a relatively small number of computations are needed toobtain a (very limited) result; this small tutorial example is discussed with the aim ofmaking our work more easy to reproduce for a reader that is interested in exporting ourapproach in other challenging contexts.
In this section we are going to describe the construction of the Birkhoff normal form fora particular class of Hamiltonians. This is made by using Lie series operators. Althoughsuch a method is rather standard in perturbation theory, a detailed description of theformal algorithm is mandatory, in order to make well definite the discussion in the nextsections. The adoption of an approach based on Lie series makes our work further differentwith respect to those that have been mentioned in the Introduction and deal with effectivestability estimates based on Birkhoff normal forms, because they implemented algorithmsusing Lie transforms.Let us give some definition. Let f be a polynomial function in the canonical variables( y , x ) ∈ R n × R n , then we define the Lie derivative operator with respect to f as L f ( · ) = omputer-assisted estimates for Birkhoff normal forms { f, ·} , being {· , ·} the Poisson bracket between two functions. Moreover, the Lie seriesoperator having f as generating function is defined as follows:exp L f ( · ) = + ∞ X j =0 j ! L jf ( · ) (1)and it will be used to define near to the identity canonical transformations. Moreover, wedenote by P s the class of homogeneous polynomials of degree s in the canonical variables( y , x ). In the following lemma, we describe the behaviour of the homogeneous polynomialswith respect to the Poisson brackets. Lemma 2.1
Let f and g be polynomial functions in P r +2 and P s +2 , respectively. There-fore, { f, g } ∈ P s + r +2 . Moreover, L jf g ∈ P jr + s +2 ∀ j ≥ . This simple property is fundamental, because it will be repeatedly used for reorganizingthe Hamiltonian terms (appearing in an expansion) with respect to their polynomialdegree, after each canonical change of coordinates defined by a Lie series operator.
Since we aim to study the dynamics in the neighbourhood of an elliptic equilibrium point,located in correspondence with the origin, we consider a Hamiltonian of the following type: H ( y , x ) = n X j =1 ω j y j + x j + ∞ X ℓ =1 f ℓ ( y , x ) with f ℓ ∈ P ℓ +2 , ω ∈ R n . (2)In the proximity of the origin, the magnitude order of the polynomial terms f ℓ is obviouslysmaller than that of the main part of the Hamiltonian, which is given by a sum of harmonicoscillators. Moreover, if we introduce the action-angle coordinates for harmonic oscillators,by the following canonical change of variables x j = − p I j cos ( ϕ j ) y j = p I j sin ( ϕ j ) ∀ j = 1 , . . . , n, (3)then the Hamiltonian in (2) transforms to H ( I , ϕ ) = ω · I + + ∞ X ℓ =1 f ℓ ( I , ϕ ) with ( I , ϕ ) ∈ R n × T n , (4)with f ℓ containing terms having total degree in the square root of the actions I , . . . , I n equal to ℓ and trigonometric expansions in ϕ with Fourier harmonics k such that | k | = P j | k j | ≤ ℓ . These variables highlight the integrability of the quadratic part, whichdepends on the actions only. Therefore, we will proceed in a perturbative way, by leadingthe Hamiltonian to a Birkhoff normal form. This means that, after having performed r Let us recall the definition of the Poisson bracket: { f, g } = P nj =1 (cid:16) ∂f∂x j ∂g∂y j − ∂f∂y j ∂g∂x j (cid:17) . C. Caracciolo, U. Locatelli canonical changes of coordinates defined by the corresponding Lie series, the Hamiltonianwill be such that H ( r ) ( I , ϕ ) = Z ( r ) ( I ) + R ( r ) ( I , ϕ ) , (5)where the upper index r in every Hamiltonian term refers to the number of normalizationsteps that have been already performed. In some more words, we want to remove stepby step the angular dependence from the perturbative part of the Hamiltonian, in orderto decrease the size of the remainder R ( r ) . This is done to provide a better integrableapproximation Z ( r ) with respect to the problem we are studying.In order to deal with the construction of the normal form and the subsequent schemeof estimates, it is convenient to work with the complex canonical variables ( − i z , ¯ z ),introduced by the following canonical transformation z j = − p I j e − iϕ j ∀ j = 1 , . . . , n . (6)In these new variables I j = z j ¯ z j , then a new normal form term should contain onlymonomials where z j and ¯ z j appear with the same exponent. Moreover, since the poly-nomial degree is the same if we are using either the canonical variables ( y , x ) or thecomplex ones ( − i z , ¯ z ), we will maintain the symbol P ℓ , ∀ ℓ ≥
0, to denote every class ofhomogeneous polynomials. r -th normalization step We are going to describe in detail how to perform a single step of the normalizationalgorithm. Let us suppose that we have been able to perform the first r − H ( r − ( − i z , ¯ z ) = r − X ℓ =0 Z ℓ ( z ¯ z ) + + ∞ X ℓ = r f ( r − ℓ ( − i z , ¯ z ) , (7)where Z ℓ ∈ P ℓ +2 is depending on the actions only and f ( r − ℓ ∈ P ℓ +2 , while the upper index r − r = 1, this is nothing but Hamiltonian (4) once it has been expressed incomplex variables. We want to normalize the main term among the perturbative ones,that is f ( r − r . For this purpose, we introduce the new r -th Hamiltonian as follows : H ( r ) ( − i z , ¯ z ) = exp L χ r H ( r − ( − i z , ¯ z ) = + ∞ X j =0 j ! L jχ r H ( r − ( − i z , ¯ z ) . (8) In formula (8), we make use of the so called exchange theorem (see, e.g., formula (4.6) in [18]) whichstates that, if χ is a generating function, f ( p , q ) | ( p , q )=exp L χ ( p ′ , q ′ ) = exp L χ f | ( p , q )=( p ′ , q ′ ) . Therefore, wecan apply the Lie series to the Hamiltonian function and, only at the end, rename the variables. For moredetailed explanations we defer to the whole section 4.1 of [18]. Here, we do not rename the variables inorder to avoid the proliferation of too many symbols. omputer-assisted estimates for Birkhoff normal forms r + 2 is depending juston the actions, then we have to determine the generating function χ r in such a way thatthe following homological equation is satisfied: L χ r Z + f ( r − r = Z r , (9)where Z is the quadratic part, i.e., Z = P nj =1 ω j z j ¯ z j , and Z r is the new term in normalform we have to define. Hereafter, it is convenient to adopt the standard multi-indexnotation; for instance, this means that ( − i z ) ℓ = Q nj =1 ( − i z j ) ℓ j , ∀ z ∈ C n and ℓ ∈ N n .Moreover, for what concerns the complex conjugate variables, the notation has to be readas follows: ¯ z ˜ ℓ = Q nj =1 ¯ z ˜ ℓ j j , ∀ z ∈ C n and ˜ ℓ ∈ N n . Let us write the generic expansion ofthe perturbative term f ( r − r so that f ( r − r = X | ℓ | + | ˜ ℓ | = r +2 c ( r − ℓ , ˜ ℓ ¯ z ˜ ℓ ( − i z ) ℓ with c ( r − ℓ , ˜ ℓ ∈ C , ℓ , ˜ ℓ ∈ N n ; (10)therefore, the generating function that solves the homological equation (9) is such that χ r = X | ℓ | + | ˜ ℓ | = r +2 ℓ =˜ ℓ − c ( r − ℓ , ˜ ℓ i ω · ( ℓ − ˜ ℓ ) ¯ z ˜ ℓ ( − i z ) ℓ , (11)where ω is the vector of frequencies of the unperturbed harmonic oscillators in (2). Clearly,this generating function can be properly defined if and only if the frequency vector ω isnon-resonant, otherwise ω · ( ℓ − ˜ ℓ ) could vanish. Therefore, we assume that ω satisfiesthe Diophantine inequality, that is | k · ω | ≥ γ | k | τ ∀ k ∈ Z n \ { } , (12)for some fixed values of γ > τ ≥ n −
1. This condition is not so strict, since for τ > n − R n . The terms with ℓ = ˜ ℓ ,that cannot be removed by this procedure, depend on the actions only, since z and ¯ z appear with the same exponent. As a consequence, we are forced to include them in thenew term Z r making part of the normal form, that is Z r = X | ℓ | = r +2 c ( r − ℓ , ℓ ( − i z ¯ z ) ℓ . (13)From definitions (11) and (13), it is evident that both χ r and Z r belong to P r +2 .In order to conclude the normalization step, we have to perform the canonical changeof coordinates induced by exp L χ r and to update accordingly all the terms that composethe Hamiltonian. For this purpose, it is convenient to use a notation which mimics aprogramming code. At first, we put the new terms f ( r ) ℓ = f ( r − ℓ ; this simply correspondsto the initial summand with j = 0 in (1). Then, each contribution due to the j -thiteration of the Lie derivative with respect to the generating function χ r (for j ≥
1) is
C. Caracciolo, U. Locatelli added to the corresponding class. More precisely, using repeatedly the property describedin lemma 2.1, the new summands, which are generated by the application of the Lie seriesto the terms in normal form (that are denoted with Z s ) and to the perturbative ones (i.e., f ( r − s ), are gathered in the following way: f ( r ) s + jr ← ֓ j ! L jχ r Z s ∀ ≤ s < r, j ≥ ,f ( r ) s + jr ← ֓ j ! L jχ r f ( r − s ∀ s ≥ r, j ≥ , (14)where the notation a ← ֓ b means that a is redefined so as to be equal to the sum givenby its previous value plus b . At the end of all these redefinitions, the following expansionof the new Hamiltonian is well defined: H ( r ) ( − i z , ¯ z ) = r X ℓ =0 Z ℓ ( z ¯ z ) + + ∞ X ℓ = r +1 f ( r ) ℓ ( − i z , ¯ z ) . (15)Let us emphasize that in this procedure we do not modify the terms that are alreadyin normal form. This is why the integrable part has nearly the same expression when (7)is compared with (15), with just one exception due to the occurrence of a new normalform term Z r , while the perturbative terms are different, being redefined as it has beenprescribed by formula (14).The algorithm can be iterated at the next step, by restarting from the Hamiltonian (7)where r − r . In spite of the fact that we are interested in providing results holding true for a ballof initial conditions B ̺ ( ), where ̺ is fixed , now we are going to briefly discuss someasymptotic properties of the Birkhoff normal form.It is well known that the series introduced by the algorithm constructing the normalform are asymptotically divergent, with the exception of some particular cases. Thismeans that the sup-norm of the remainder R ( r ) (appearing in (5)) does not go to zero for r which tends to infinity. However, the Birkhoff normal form is still very useful, becausethere is an optimal normalization step which minimizes the remainder. The mechanismof divergence is mainly due to the accumulation of the so called “small divisors”. In fact,at each step r , new divisors ω · ( ℓ − ˜ ℓ ) (being | ℓ | + | ˜ ℓ | = r + 2) are introduced by thedefinition of the generating function χ r in (11). Because of the redefinitions describedin formula (14), the size of the generating functions clearly impacts on the growth ofthe terms constituting the remainder. If the frequency vector satisfies the Diophantinecondition in (12), it is rather easy to see that factorial coefficients O ( r !) appear in theestimate of the remainder at the r -th normalization step. This prevents the convergenceof the normalization algorithm when it is iterated ad infinitum .Nevertheless, when the series are estimated on open balls B ̺ ( ), with a fixed valueof ̺ , it is convenient to proceed with the normalization algorithm until the remainder is omputer-assisted estimates for Birkhoff normal forms r opt is determined, by comparing two subsequent remainders R ( r ) and R ( r − . As discussed, e.g., in [18] (see formula (3.30)), this allows to concludethat r opt ∼ ( C̺ ) − τ +1 , (16)where C is a positive constant, τ is the exponent appearing in the Diophantine condi-tion (12) and ̺ is the distance from the equilibrium point. As a consequence, an analyticalestimate of the remainder at the optimal step can be easily provided; in fact, it is expo-nentially small with respect to the inverse of a fractional power of the distance from theequilibrium ̺ , that in this case assumes the role of small parameter. Near the ellipticequilibrium, the Birkhoff normal form is a good approximation of the real problem andit is convenient to perform a big number of normalization steps. On the other hand, faraway from the origin the Birkhoff normal form starts to diverge earlier. Let us mentionthat in [13] the mechanism of divergence is carefully investigated in a numerical way andit is shown to be much more subtle with respect to what has been discussed just above.Moreover, this has allowed those authors to conclude that the analytical upper boundsshould largely overestimate the effective size of the remainders.In practice, when dealing with explicit expansions, two additional problems have tobe tackled: truncations and computational costs. In fact, the Hamiltonian at step r isexpanded as in (5) and (15), with R ( r ) = f ( r ) r +1 + f ( r ) r +2 + . . . with f ℓ ∈ P ℓ +2 ∀ ℓ > r. (17)Of course, efficient coding and a powerful computing system are more than welcome,because they allow better evaluations of the remainder terms. However, any computercannot represent all the infinite sequence of terms giving contributions to R ( r ) . In thenext section, we will explicitly provide rigorous upper bounds for the series defining theremainder. In this section we are going to discuss how to construct a scheme of estimates, in order tokeep control of the norms of the Hamiltonian terms, starting from the algorithm describedin section 2. In spite of the fact that every purely analytical work based on Birkhoffnormal forms exploits a suitable scheme of estimates, our approach significantly differsfrom all the previously existing ones in the scientific literature for what concerns thefollowing untrivial point, that is clearly highlighted, for instance, in the statements ofpropositions 3.3 and 3.4. Here, the estimates are designed in such an iterative way, thatthey take advantage of being applied in the framework of a computer-assisted proof, inorder to provide the best possible final results.Before proceeding with the description of the scheme of estimates, the introduction ofsome new notation is mandatory. For every f ∈ P s whose expansion is the following: f = X | ℓ | + | ˜ ℓ | = s c ℓ , ˜ ℓ ( − i z ) ℓ ¯ z ˜ ℓ , (18)0 C. Caracciolo, U. Locatelli we denote by k · k the functional norm defined as k f k = X | ℓ | + | ˜ ℓ | = s | c ℓ , ˜ ℓ | . (19)Therefore, in the domain∆ ̺ = { ( − i z , ¯ z ) : | z j | < ̺ , ≤ j ≤ n } = { ( I , ϕ ) : I j < ̺ , ϕ j ∈ T , ≤ j ≤ n } , the sup norm of f can be easily estimated as follows: | f | ̺ = sup ∆ ̺ | f ( − i z , ¯ z ) | ≤ ̺ s k f k . (20)For the sake of simplicity, in the previous definitions we have not introduced moregeneral domains that are polydisks; this has been done with the aim to avoid all themodifications that are necessary to adapt norms and estimates. However, let us recallthat in some applications of the Birkhoff normal forms the approach based on polydiskshas been useful (see, e.g., [24]).Hereafter, we will assume to have explicitly built the Birkhoff normal form up to afixed order r = R I and we will define estimates for the norms of the infinite terms we didnot calculate. In order to do that, for every normalization step we will distinguish betweentwo classes of terms making part of the power series expansion: those belonging to P ℓ with ℓ ≤ R II , for which we will give explicit estimates of the norms, and the infinite termswith degree greater than R II + 2, that we will control by a sequence of majorants growingin a geometrical way. In the next subsections, we describe how to define iteratively theestimates for the norms of these two kinds of polynomial terms. Let us suppose that we have already performed the first r − k Z s k ≤ Z s ∀ ≤ s ≤ r − , k f ( r − s k ≤ F ( r − s ∀ s ≥ r , (21)with Z s , F ( r − s ∈ R + . We have now to describe how to estimate the norms of the termsthat compose the new Hamiltonian H ( r ) , namely how the estimates F ( r − s change after astep of Birkhoff normalization. Since the canonical change of coordinates that transforms H ( r − into H ( r ) is defined by a Lie series with generating function χ r , we need an estimatefor the norm of the Poisson brackets between two functions. In the following lemma, weadapt to the present context a well known estimate that can be found, e.g., in [19]. Lemma 3.1
Let f and g be polynomial functions such that f ∈ P s +2 , g ∈ P r +2 andsuppose that k f k ≤ F and k g k ≤ G for some constant G , F ∈ R + ; therefore, k { f, g } k ≤ ( s + 2)( r + 2) F G . (22) omputer-assisted estimates for Birkhoff normal forms Lemma 3.2
Let χ r and g be polynomial functions with χ r ∈ P r +2 and g ∈ P s +2 , beingthe expansion of the former function such that χ r = X | k | + | ˜ k | = r +2 c k , ˜ k ( − i z ) k ¯ z ˜ k . (23) Therefore, ∀ j ≥ , the following estimate holds true: (cid:13)(cid:13)(cid:13)(cid:13) j ! L jχ r g (cid:13)(cid:13)(cid:13)(cid:13) ≤ Q j − i =0 ( s + ir + 2) j ! D jr k g k , (24) where D r = P | k | + | ˜ k | = r +2 | c k , ˜ k | max j {| k j | , | ˜ k j |} . Indeed, if we replace D r with ( r + 2) k χ r k , the previous statement easily follows by in-duction from lemma 3.1. The above definition of D r allows to estimate more carefullythe partial derivatives of χ r that appear in the Poisson brackets. This is particularlyadvantageous when the optimal normalization step is not greater than the maximal value R I of the index r for which we compute explicitly the expansion of the generating function χ r . On the other hand, when the number R I is lower than the optimal step r opt , it isconvenient to provide an estimate for the norm k χ r k ≤ G r and, then, D r can be replacedby ( r + 2) G r , for all those r ∈ [ R I + 1 , r opt ] for which the coefficients c k , ˜ k appearing inthe expansion (23) are unknown.In order to determine upper bounds for the norms of the new perturbative terms(introduced by the r -th normalization step of the algorithm), we have to focus on theredefinitions in formula (14) for using inequality (24). Therefore, it is easy to realizethat we substantially need to find an estimate for the norm of the generating function.Obviously, if r ≤ R I , we know exactly the expression of χ r and we can calculate D r asprescribed in lemma 3.2. When the normalization step r > R I , recalling the expansion of χ r in (11), we can say that the coefficients of the generating function are equal to thoseof f ( r − r divided by ω · ( ℓ − ˜ ℓ ). Let α r be the smallest divisor introduced at the r -thnormalization step, i.e., α r = min | ℓ | + | ˜ ℓ | = r +2 ℓ =˜ ℓ | ω · ( ℓ − ˜ ℓ ) | , (25)which is well-defined in view of the non-resonance condition assumed in (12). Therefore,we can write k χ r k ≤ G r , with G r = F ( r − r α r . (26)Analogously, by definition of the new term in normal form Z r in (13), the norm of Z r cannot be greater than Z r = F ( r − r . We are now ready to define the new majorants F ( r ) s ,providing upper bounds for the perturbative terms appearing in the Hamiltonian H ( r ) .2 C. Caracciolo, U. Locatelli
Proposition 3.3
Let us suppose that the terms of the Hamiltonian introduced at r − -thnormalization step are bounded as in (21) . Therefore, the new Hamiltonian terms of H ( r ) are bounded so that k Z r k ≤ Z r = F ( r − r and k f ( r ) s k ≤ F ( r ) s ∀ s > r , where these newconstants are given by the following sequence of redefinitions: F ( r ) ℓ = F ( r − ℓ ∀ ℓ ≥ , F ( r ) s + jr ← ֓ Q j − i =0 ( s + ir + 2) j ! D jr F ( r − s ∀ j ≥ , s ≥ , (27) being D r = P | k | + | ˜ k | = r +2 | c k , ˜ k | max j {| k j | , | ˜ k j |} , if the expansion of the generating function χ r is explicitly known, as it is written in (23) , else D r = ( r + 2) G r , with G r = F ( r − r /α r . The previous statement is basically a summary of all the discussion explained in thepresent subsection. In particular, formula (27) can be easily verified, starting from theiterative definitions of the new perturbative Hamiltonian terms in (14) and using theestimates in lemma 3.2.
In this subsection, we are going to describe how to estimate the norm of the terms thatcannot be explicitly calculated or iteratively estimated, as it has been described in theprevious subsection. In order to do that, it is convenient to bound in a different way theterms constituting the Hamiltonian at the r − k Z s k ≤ E a s +2 r − ∀ ≤ s ≤ r − , k f ( r − s k ≤ E a s +2 r − ∀ s ≥ r , (28)with E , a r − ∈ R + . Of course, the estimates above are expected to be less strict than theones in (21). However, let us imagine to be able to justify the same kind of estimates atthe r -th normalization step, holding true for all f ( r ) s with s ≥ r + 1, with a new constantvalue a r ; therefore, we could control the sup norm of the remainder (recall equations (5)and (15)) as follows: sup ∆ ̺ |R ( r ) | ≤ R II X s = r +1 F ( r ) s ̺ s +2 + E + ∞ X s = R II +1 ( a r ̺ ) s +2 , (29)where we include the upper bounds as iteratively defined by proposition 3.3 when theyare available, while for the terms with polynomial degree greater than R II + 2 we use thewanted geometric estimates. These are obtained by iteration of the rule described below. Proposition 3.4
Let us assume that the estimates in (28) hold true for all the termsappearing in the expansion (7) of H ( r − . Therefore, the norms of the terms making partof the new Hamiltonian H ( r ) satisfy the following inequalities: k Z s k ≤ E a s +2 r ∀ ≤ s ≤ r , k f ( r ) s k ≤ E a s +2 r ∀ s ≥ r + 1 , (30) omputer-assisted estimates for Birkhoff normal forms being a r = a r − (cid:18) r + 1) D r a rr − (cid:19) r , (31) where D r = P | k | + | ˜ k | = r +2 | c k , ˜ k | max j {| k j | , | ˜ k j |} if the expansion (23) of the generating func-tion χ r is explicitly known, else D r = ( r + 2) G r with k χ r k ≤ G r . Proof
Let us start from the estimate for Z s with 1 ≤ s < r . Since a r − < a r , all theestimates in (28) are still valid when a r − is replaced by a r . By comparing the polynomialexpressions (10) and (13) that are related to the functions f ( r − r and Z r , respectively, it iseasy to realize that the norm of the latter cannot be greater than the former one, becausein a normal form term there are just a part of the summands appearing in expansion of thecorresponding perturbative function. Therefore, in view of (28) we can write the followingchain of inequalities: k Z r k ≤ k f ( r − r k ≤ E a r +2 r − ≤ E a r +2 r ; this ends the justification of thefirst inequality in formula (30).Let us now focus on the estimates for f ( r ) s , with s ≥ r + 1. First, we discuss the casewith s = jr + m for m = 1 , . . . , r − j ≥
1: the new perturbative terms are definedas the sum of f ( r − s plus the new contributions due to L jχ r Z m /j ! and L ℓχ r f ( r − j − ℓ ) r + m /ℓ ! for ℓ = 1 , . . . , j −
1. Therefore, inequality (24) allows us to write the following estimate: F ( r ) jr + m = F ( r − jr + m + 1 j ! L jχ r Z m + j − X ℓ =1 ℓ ! L ℓχ r f ( r − j − ℓ ) r + m ≤ F ( r − jr + m + Q j − i =0 ( m + ir + 2) j ! D jr Z m + j − X ℓ =1 Q ℓ − i =0 ( m + ( j − ℓ + i ) r + 2) ℓ ! D ℓr F ( r − j − ℓ ) r + m ! . Using the estimate in hypothesis (30) and the trivial inequality m < r , we obtain F ( r ) jr + m ≤ E a jr + m +2 r − + Q j − i =0 ( i + 1) j ! [( r + 1) D r ] j E a m +2 r − + j − X ℓ =1 Q ℓ − i =0 ( j − ℓ + i + 1) ℓ ! a ℓrr − [( r + 1) D r ] ℓ E a jr + m +2 r − . Therefore, F ( r ) jr + m ≤ E a jr + m +2 r − " r + 1) D r ] j + j − X ℓ =1 (cid:18) jℓ (cid:19) [( r + 1) D r ] ℓ a rℓr − = E a jr + m +2 r − j X ℓ =0 (cid:18) jℓ (cid:19) (cid:18) ( r + 1) D r a rr − (cid:19) ℓ ≤ E a jr + m +2 r − (cid:18) r + 1) D r a rr − (cid:19) j ≤ E " a r − (cid:18) r + 1) D r a rr − (cid:19) r jr + m +2 . C. Caracciolo, U. Locatelli
This concludes the proof of the case s > r , not being a multiple of r .Let us now consider s = jr with j ≥
2. In such a case, it is convenient to rewriteformula (14) in the following more extended form: f ( r ) jr ← ֓ j ! L jχ r Z + 1( j − L j − χ r f ( r − r + j − X ℓ =1 ℓ ! L ℓχ r f ( r − j − ℓ ) r , where we have separated the term generated by f ( r − r . Since χ r solves the homologicalequation (9), we can rewrite the first two contributions as follows:1 j ! L jχ r Z + 1( j − L j − χ r f ( r − r = 1 j ! L j − χ r ( L χ r Z + f ( r − r ) + j − j ! L j − χ r f ( r − r = 1 j ! L j − χ r Z r + j − j ! L j − χ r f ( r − r . Therefore, for j ≥ F ( r ) jr ≤ F ( r − jr + 1 j (cid:13)(cid:13)(cid:13)(cid:13) j − L j − χ r Z r (cid:13)(cid:13)(cid:13)(cid:13) + j − j (cid:13)(cid:13)(cid:13)(cid:13) j − L j − χ r f ( r − r (cid:13)(cid:13)(cid:13)(cid:13) + j − X ℓ =1 (cid:13)(cid:13)(cid:13)(cid:13) ℓ ! L ℓχ r f ( r − j − ℓ ) r (cid:13)(cid:13)(cid:13)(cid:13) ≤ F ( r − jr + 1 j · Q j − i =0 ( r + ir + 2)( j − D j − r k Z r k + ( j − j · Q j − i =0 ( r + ir + 2)( j − D j − r (cid:13)(cid:13) f ( r − r (cid:13)(cid:13) + j − X ℓ =1 Q ℓ − i =0 (( j − ℓ ) r + ir + 2) ℓ ! D ℓr (cid:13)(cid:13)(cid:13) f ( r − j − ℓ ) r (cid:13)(cid:13)(cid:13) ≤ F ( r − jr + 1( j − j − Y i =0 [(1 + i ) r + 2)] D j − r F ( r − r + j − X ℓ =1 ℓ ! ℓ − Y i =0 [( j − ℓ + i ) r + 2] D ℓr F ( r − j − ℓ ) r ≤ j − X ℓ =0 ℓ ! ℓ − Y i =0 [( j − ℓ + i ) r + 2] D ℓr F ( r − j − ℓ ) r , where we used again inequality (24) and F ( r − r as upper bound for the norm of Z r .The second inequality in formula (30) follows from calculations similar to those we havepreviously described for the case where s is not a multiple of r . (cid:4) Let us emphasize that the same scheme of estimates is valid also for a resonant Birkhoffnormal form, since we did not use any explicit definition of the terms Z s appearing in (15).Therefore, in order to cover the extension to such a case, it will be enough to adapt thedefinition of the smallest divisor α r and, as a consequence, the estimate for the generatingfunction χ r . omputer-assisted estimates for Birkhoff normal forms Here we are going to use the scheme of estimates described in the previous section forproducing results on the effective stability in the neighborhood of an elliptic equilibriumpoint. We will study separately two different cases: a non-resonant Birkhoff normal formand a resonant one. For what concerns the effective stability in the vicinity of equilibriumpoints, up to the best of our knowledge other computer-assisted rigorous estimates areavailable in [8] only, but such a reference has to be considered rather outdated, becauseof some technical reasons affecting the final quality of the results (see also the commentsin the Introduction). Moreover, an application of this kind of approach to a resonantBirkhoff normal form is entirely new.
Let us suppose that the Hamiltonian has been already brought in Birkhoff normal form,after having performed r normalization steps, i.e., H ( r ) = Z ( r ) + R ( r ) with Z ( r ) dependingon the actions only and R ( r ) = O ( k I k ( r +3) / ). Therefore, we can calculate the timederivative of the action I j = z j ¯ z j as follows:˙ I j = (cid:8) I j , H ( r ) (cid:9) = (cid:8) I j , Z ( r ) + R ( r ) (cid:9) = (cid:8) I j , R ( r ) (cid:9) ∀ j = 1 , . . . , n , (32)because { I j , Z ( r ) ( I ) } = 0. Using the expansion (15), the estimate on the sup normin (20) and lemma 3.1, we can find the following upper bound for the sup norm of thetime derivative of the actions : | ˙ I j | ̺ = sup ∆ ̺ | ˙ I j | = sup ∆ ̺ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)( I j , + ∞ X s = r +1 f ( r ) s )(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ + ∞ X s = r +1 ̺ s +2 k (cid:8) I j , f ( r ) s (cid:9) k≤ + ∞ X s = r +1 ( s + 2) ̺ s +2 k f ( r ) s k . (33)We can now provide rigorous estimates for the norms of all the terms appearing in theprevious series. More precisely, we can use the estimates F ( r ) s (as recursively definedin (27)) for the first R II terms. Moreover, for what concerns the first R I terms, wecan do even better, because we have explicitly performed R I steps of the normalizationalgorithm; therefore, for such terms the estimates can be replaced by the actual valuesof the corresponding norms. For all the remaining terms, we limit ourselves to use theuniform estimate in (30). Thus, we just need to compute the value of the upper bound The general formula in lemma 3.1 would give a factor 2, but here we take advantage of the fact that I j = z j ¯ z j so we do not have terms with z j . C. Caracciolo, U. Locatelli in the r.h.s. of the following inequality: + ∞ X s = r +1 ( s + 2) ̺ s +2 k f ( r ) s k ≤ R I X s = r +1 ( s + 2) ̺ s +2 k f ( r ) s k + R II X s = R I +1 ( s + 2) ̺ s +2 F ( r ) s + E + ∞ X s = R II +1 ( s + 2) ̺ s +2 a s +2 r . (34)Since the remainder is expected to be exponentially small at the optimal step, it is clearthat the quantity in (33) can be very small as well. From a computational point of view,iterating estimates is much cheaper than doing algebraic calculations on huge expansions;this is why such a procedure can be very convenient when we are near the equilibriumpoint, where the optimal normalization step is far beyond the number of steps we are ableto explicitly perform, by using any kind of software that is specialized for doing computeralgebra.Let us now suppose that the initial value of the action vector belongs to an opendomain, e.g., I (0) ∈ ∆ ̺ ; we aim to determine a time T > I ( t ) ∈ ∆ ̺ ∀ t ∈ [ − T , T ], being the domain ∆ ̺ a little largerthan ∆ ̺ , i.e. ̺ > ̺ . First, let us remark that | I j ( t ) | ≤ | I j (0) | + | I j ( t ) − I j (0) | ≤ ̺ + | ˙ I j | ̺ T ∀ j = 1 , . . . , n. (35)When an estimate for | ˙ I j | ̺ ∀ j is available, it is convenient to define T as ̺ − ̺ max j | ˙ I j | ̺ (36)to obtain the wanted confinement. Combining (33) and (34) with (36), we obtain thefollowing lower bound about the escape time from ∆ ̺ : T ( ̺, ̺ , r ) = ̺ − ̺ R I X s = r +1 ( s + 2) ̺ s +2 k f ( r ) s k + R II X s = R I +1 ( s + 2) ̺ s +2 F ( r ) s + E + ∞ X s = R II +1 ( s + 2) ̺ s +2 a s +2 r , (37)where the dependency on different parameters of such an expression is emphasized;moreover, r = r opt is chosen so as to minimize the estimate (29) of the remainder R ( r ) .It is rather easy to verify that the value of ̺ making a further optimization of theexpression above is ̺ = [( r opt + 1) / ( r opt + 3)] / ̺ . When r = r opt , for the sake of simplicity, let us assume that the geometrical decrease of the termsin the series appearing at the denominator of equation (37) is so sharp that it can be approximated byits first term, then we can write T ( ̺, ̺ , r opt ) ≃ g ( ̺ ), being g ( ̺ ) = C ( ̺ − ̺ ) / [( r opt + 3) ̺ ( r opt +3) ] and C a suitable positive constant. After having remarked that g ( ̺ ) = lim ̺ →∞ g ( ̺ ) = 0, one immediatelyrealizes that the function g : [ ̺ , ∞ ) R takes its maximum value in correspondence to the solution ofthe equation g ′ ( ̺ ) = 0, i.e., ̺ = [( r opt + 3) / ( r opt + 1)] / ̺ . omputer-assisted estimates for Birkhoff normal forms ̺ is considered as fixed, we will computethe lower bound for the escape time T (cid:0) ̺ , [( r opt + 1) / ( r opt + 3)] / ̺ , r opt (cid:1) according toformula (37), where r opt makes optimal the estimate for the remainder R ( r ) .Let us recall that, in the framework of purely analytical estimates, the denominatorappearing in formula (37) would be made just by its third summand with the seriesstarting from r + 1 instead of R II + 1; moreover, one can verify that a r ∼ ( r !) τ +1 C r , being C a suitable positive constant. By choosing r = r opt as in formula (16), one can obtainthe asymptotic law for the purely analytical estimate about the escape time from ∆ ̺ , i.e., T ∼ exp "(cid:18) ˜ ̺̺ (cid:19) τ +1 , (38)where ˜ ̺ is a positive constant. For more details we refer to theorem 3.5 in [18], where acomplete proof of such a lower bound can be found. In the case of a resonant Birkhoff normal form, the estimate of the diffusion of the actionsis not so straightforward as for the non-resonant one. Indeed, we cannot use directly theremainder in order to estimate the variation of the resonant actions as we did in (33).Here, we will discuss the case in which there is only one resonant angle.Let us assume we have already performed the construction of a resonant Birkhoffnormal form up to order r , in such a way that the Hamiltonian of the system has thefollowing structure: H ( r ) ( I , ϕ ) = Z ( r ) ( I , ϕ n ) + R ( r ) ( I , ϕ ) , (39)where I ∈ R n and ϕ ∈ T n are action–angle canonical coordinates, being ϕ n the resonantangle. Therefore, the estimates in (33) and (35) still hold true for what concerns thediffusion of the actions I , . . . , I n − , since we removed the corresponding angles fromthe normal form. On the other side, it is more difficult to provide good estimates forthe diffusion of the action I n , that is conjugated to the resonant angle ϕ n . Since thecomputation of the stability time T as in formula (37) is valid when all the actions areconfined in the domain ∆ ̺ , it is necessary to control the variation of the action I n inorder to use such an estimate about the value of T . For this purpose, it is convenient tocombine the conservation of the total energy of the system, i.e., the Hamiltonian, withthe slow diffusion of the actions I , . . . , I n − . More precisely, by omitting the dependencyon the index r in formula (39), we can write the following equation: I n = E − P n − j =1 ω j I j ω n − ¯ Z ( I , ϕ n ) + R ( I , ϕ ) ω n , (40)where E is the energy level and we denote with ¯ Z the sum of all the normal form termsbut the linear ones with respect to the actions (that are still of the form ω · I as in theinitial Hamiltonian (4)). By assuming the confinement of all the actions in a ball of radius8 C. Caracciolo, U. Locatelli ̺ , we can evaluate an upper bound of I n ; max = max | t |≤ T sup ∆ ̺ ( E − P n − j =1 ω j I j ω n − ¯ Z ( I , ϕ n ) + R ( I , ϕ ) ω n ) , (41)and a lower bound of I n ; min = min | t |≤ T inf ∆ ̺ ( E − P n − j =1 ω j I j ω n − ¯ Z ( I , ϕ n ) + R ( I , ϕ ) ω n ) . (42)This allows us to give an estimate of∆ I n = I n ; max − I n ; min . (43)Of course, the excursions experienced by the values of the normal form Z mainly dependon the linear terms; this explains why we have written them separately in the equationsabove. The bounds for the sup-norm and inf-norm appearing in formulæ (41)–(42), re-spectively, can be easily calculated, by using the inequalities (21) and (28). Since thenormal form is composed by a finite number of terms, the estimate for ∆ ¯ Z can be eval-uated very strictly when its expansion is explicitly known.Recalling that the estimate for the variation of I n holds true only if all the actionsare confined in a ball of radius ̺ , the maximal radius on I n in order to validate such anestimate is defined as ( ̺ ∗ n ) = ̺ − ∆ I n . (44)Therefore, if we define T as in formula (36) for j = 1 , . . . , n −
1, we can assure that I ( t ) ∈ ∆ ̺ ∀ t ∈ [ − T, T ], provided that | I j (0) | < ̺ ∀ j = 1 , . . . , n − | I n (0) | < ( ̺ ∗ n ) .Let us remark that the more ̺ is close to ̺ , the more the variation of I n is prescribedto be small; at the same time, as it follows immediately from formula (37), the stabilitytime decreases for ̺ going to ̺ . In order to balance these two effects, it is convenientto define the parameter ̺ in a different way to what has been done at the end of theprevious subsection. In the next section, we will discuss in a practical example how theradius ̺ of the open ball containing the initial conditions can be determined in the caseof resonant Birkhoff normal forms. The H´enon-Heiles model was introduced in [27] for studying the dynamics of a star in agalaxy. The Hamiltonian that describes the system is composed of a couple of harmonicoscillators perturbed with cubic terms, i.e., H ( y , x ) = ω x + y ω x + y x x − x y , x ) ∈ R n . (45) omputer-assisted estimates for Birkhoff normal forms ω = 1 and ω = − ( √ − /
2, so that the angular velocity vector is Diophantine with τ = 1 according toits definition in (12), because the ratio of its components | ω /ω | is equal to the goldenmean . We focus on the study of the dynamics in a neighborhood of the origin. ThisHamiltonian belongs to the general class described in (2), for which we have explainedhow to perform the normalization procedure `a la Birkhoff.Using X ̺ ´o ν o ζ , which is an algebraic manipulator specially designed to implement ap-proaches that are common in the framework of Hamiltonian perturbation theory (see [23]),we explicitly constructed the Birkhoff normal form for this model by representing termshaving total polynomial degree less than or equal to 102 and by performing R I = 100normalization steps. Using the iterative estimates described in section 3, combined withthe norms that have been determined for the first R I = 100 terms (for which we haveexplicit expansions), we provided rigorous estimates for the remainder and the time de-rivative of the actions, for any fixed value of ̺ . Indeed, for the first R II = 1500 terms weiterated the estimates of the norms as described in proposition 3.3, while for the rest ofthe infinite terms we used the upper bounds reported in formula (28), updating the valueof a r as prescribed in (31); finally, we produced an estimate for both the remainder andthe escape time T , by using formulæ (29) and (37). For each value of ̺ we considered,the computational algorithm was stopped, after having identified the optimal step r opt asthe one corresponding to the minimum value of the remainder; thus, we fixed our finalestimate for the escape time in such a way that T = T ( ̺, [( r opt + 1) / ( r opt + 3)] / ̺, r opt ),in agreement with the discussion at the end of subsection 4.1.Let us recall that a similar technique, with the same meaning of the integer paramet-ers R I and R II , allowed to obtain a fully rigorous computer-assisted proof of existence forthe KAM torus related to the golden mean ratio of frequencies in the case of the forcedpendulum with a value of the small parameter that is ∼
92% of the breakdown threshold(see [9]); as far as we know, this is still the best result of such a very particular kind,for what concerns the applications of KAM theory to a Hamiltonian continuous flow. Inorder to fit with a so successful approach, in our programming codes we have implementedvalidated numerics exactly in the same way. Therefore, every coefficient appearing in allthe polynomial expansions that are explicitly computed (for Hamiltonians or generatingfunctions) are replaced with intervals and all the mathematical operations between inter-vals are performed by taking into account the round-off errors. For what concerns theiteration of the estimates, validated numerics is used in order to properly provide all theneeded bounds. In summary, the whole computational procedure is aiming at ensuringthat the final result is fully rigorous. Appendix A is devoted to a short introduction tovalidated numerics, that is described with the only goal of explaining our implementa-tion of it. Furthermore, in appendix B we have included a pedagogical discussion of acomputer-assisted estimate of the stability time for the H´enon-Heiles model in a verysimple situation: we have focused on a so small neighborhood of the equilibrium point(i.e., we have fixed ̺ = 0 . The prominent role exerted by the so called noble numbers is highlighted in [36]; the golden mean isthe main representative of such a class of numbers, forming a subset of the Diophantine ones with τ = 1. C. Caracciolo, U. Locatelli ̺ ̺ r opt a r log |R ( r opt ) | ̺ log | ˙ I j | ̺ log T ω = 1 and ω = − ( √ − / ̺ from the equi-librium point: let us emphasize that each row corresponds to a single computer-assistedproof. Therefore, we think that every item of that list can be useful for comparisons withthe results eventually obtained by using other techniques. Moreover, in Figure 1 we havereported the main content of Table 1, that is about the behaviour of both the optimalstep and the estimate for the escape time as functions of ̺ . In the left box of Figure 1 onecan appreciate that the limit value of ̺ for which we are able to produce stability resultsis around 0 . a r , whose values get close to 10just after a few normalization steps; this prevents the convergence of the normal form,when ̺ & .
1, because the common ratio of the majorants series is too large: a r ̺ ≥ r opt ( ̺ ), that is incorrespondence with the number of normalization steps R I = 100 for which we explicitlycompute the expansions. Indeed, the worsening effect induced by the transition to themere iteration of the estimate is so remarkable that we have to strongly decrease the ballradius ̺ in order to take a real advantage of performing more normalization steps. With omputer-assisted estimates for Birkhoff normal forms r op t ρ T ρ -1/2 Figure 1: On the left, plot of the optimal normalization step r opt as a function of theball radius ̺ ; on the right, graph of the evaluation of our lower bound about the escapetime T as a function of 1 / √ ̺ . Both the plots refer to results obtained by applyingcomputer-assisted estimates to the H´enon-Heiles model with frequencies ω = 1 and ω = − ( √ − / r opt ( ̺ ) ∼ / √ ̺ (see for-mula (16)). The agreement between the observed behavior and the expectations is evenbetter for what concerns the plot in the right box of Figure 1. In fact, from the asympoticlaw (38) one can deduce that log T ∼ / √ ̺ , that is coherent with the fact that such asemi-logarithmic plot fits rather well with a straight line. Thus, we can conclude thatour computer-assisted estimates about (the lower bound of) the stability time preserveits main property: being exponentially large with respect to the inverse of the square rootof the distance from the equilibrium point.Let us conclude this subsection with a short discussion about the choice of the para-meters R I and R II . Obviously, it is convenient to perform the largest possible numberof explicit steps. Therefore, our choice of the parameter R I is mainly due to the com-putational cost in time and memory that is needed by our algebraic manipulations. Thechoice of the parameter R II is more delicate, because it somehow rules the size of the“infinite tail” of terms in the series. For small values of ̺ , the main contribution tothe remainder (29) is due to the low order terms and the series of common ratio a r ̺ does not affect the result in an appreciable way, even for small values of R II . In sucha situation, choosing larger values of R II has the only effect to increase the computa-tional time, without any substantial improvement for what concerns the estimate of theremainder. Conversely, for larger values of ̺ , the common ratio a r ̺ can approach 1 andthe contribution of the infinite tails becomes predominant. Therefore, it is convenientto increase R II in order to estimate more properly the remainder. This explains why, inour calculations, we have usually fixed R II = 1500, but we have suitably increased sucha value when a r ̺ .
1. However, the computational cost of the iteration of the estim-ates is negligible; to fix the ideas, when R II = 1500 it took less than a minute, whileperforming explicitly R I = 100 normalization steps required about one day of CPU-time2 C. Caracciolo, U. Locatelli on a computer equipped with an
Intel quad-core I5-6600 - 3.3 Ghz and
16 GB ofRAM . Let us recall that such an impressive difference of computational cost is due to thefact that the representation of the polynomials defined by the normalization algorithmrequires a huge occupation of the memory if compared with a single upper bound on thecorresponding norm and, therefore, many more operations are needed to manipulate theirexplicit Taylor expansions. For planning new applications it is also important to know thescaling law of the computational time T CPU as a function of R I : we have found that forthe first 100 normalization steps such a law is not extremely sharp but it can be boundedfrom top (and conveniently approximated) in such a way that T CPU ∼ R . Let us briefly recall the Hamiltonian model we start from (that is described, e.g., in [24]and [17]), in order to study the long term stability in the vicinity of the triangular Lag-rangian equilibria for the CPRTBP. As a first preliminary step, a rotating frame
Oxy withorigin O in the centre of mass of the two primary bodies and rescaled physical units isintroduced. The Hamiltonian of such a system is usually written as follows (see, e.g., [47]): H ( p x , p y , x, y ) = p x + p y yp x − xp y − − µ p ( x − µ ) + y − µ p ( x + 1 − µ ) + y , where µ is the mass ratio between the primaries and ( p x , p y ) ∈ R are the kinetic momentathat are canonically conjugate to the positions ( x, y ) ∈ R , respectively. Therefore, it isconvenient to perform a preliminary transformation to heliocentric polar coordinates;afterwards, we introduce local coordinates ( P X , P Y , X, Y ) centered in a triangular point.They are defined in such a way that x = ( X + 1) cos (cid:16) Y ∓ π (cid:17) + µ , p x = P X cos (cid:16) Y ∓ π (cid:17) − sin (cid:0) Y ∓ π (cid:1) X + 1 ( P Y + 1) ,y = ( X + 1) sin (cid:16) Y ∓ π (cid:17) , p y = P X sin (cid:16) Y ∓ π (cid:17) + cos (cid:0) Y ∓ π (cid:1) X + 1 ( P Y + 1) + µ , where minus and plus refer to the equilibrium points L and L , respectively. In thesenew variables, the Hamiltonian takes the form H ( P X , P Y , X, Y ) = 12 (cid:20) P X + ( P Y + 1) ( X + 1) (cid:21) − P Y − µ ( X + 1) cos (cid:16) Y ∓ π (cid:17) − − µX + 1 − µ r ( X + 1) + 1 + 2( X + 1) cos (cid:16) Y ∓ π (cid:17) . After having skipped a constant term, the basic Taylor expansion of the previous Hamilto-nian can be written as H ( P X , P Y , X, Y ) = Q ( P X , P Y , X, Y )+ O (cid:0)(cid:12)(cid:12) ( P X , P Y , X, Y ) k (cid:1) , wherethe terms that are quadratic with respect to the canonical variables are gathered in Q .It is well known that if µ is smaller than the so called Routh-Gascheau critical value omputer-assisted estimates for Birkhoff normal forms − √ /
18, then the equation of motion are linearly stable. This means that there isa class of linear canonical transformations ( P X , P Y , X, Y ) = C ( y , y , x , x ) conjugatingthe quadratic approximation to a couple of harmonic oscillators, i.e., Q (cid:0) C ( y , y , x , x ) (cid:1) = ν ( x + y ) / ν ( x + y ) /
2. Therefore, the Taylor series expansion of the Hamiltonian H = H ◦ C can be written in the same form described in (2), that is suitable to start thenormalization procedure `a la Birkhoff, i.e., H ( y , y , x , x ) = ν x + y ν x + y + ∞ X ℓ =3 f ℓ ( y , y , x , x ) with f ℓ ∈ P ℓ , (46)being the angular velocities ν and ν such that ν = s p µ − µ + 12 and ν = − s − p µ − µ + 12 . (47)As an immediate consequence, we have that ν → ν = O ( µ ) for µ →
0. The secondfrequency is therefore much slower than the first one, for small values of the mass ratio µ between the primaries. Let us remark that during the standard procedure constructingthe Birkhoff normal form (that has been widely discussed in section 2), very small divisorscan be introduced because of the fact that | ν | ≪
1. Therefore, it is natural to expect thata resonant Birkhoff normal form, aiming to remove just the first angle (i.e., the fastest)can be more advantageous. Let us also recall that the bodies orbiting around a triangularequilibrium point are commonly called “Trojans”, according to the tradition started byMax Wolf at the beginning of the XX century. In fact, he decided to choose names fromHomer’s Iliad for the first bodies which were observed in the vicinity of L or L , in thesystem having Sun and Jupiter as primary bodies.In the next subsection we will compare the performances of these two kinds of Birkhoffnormal forms, by considering a few realistic values for the parameter µ , which correspondto systems having the following couples of primary bodies: Sun–Jupiter, Sun–Uranus,Sun–Mars and Saturn–Janus. In particular, we aim to prove the stability for a time thatis comparable with an overestimate of the residual life of the Sun in the main sequence,i.e., 6 × years. As it has been discussed in the introduction, an approach merely based on the Birkhoff nor-mal form is not enough to obtain realistic results about the stability of the Jupiter Trojansin the framework of the CPRTBP model. In [24], it has been shown that orbital motionsstarting from initial conditions contained in the domain ( I , I ) ∈ [0 , . × [0 , . L and L are so populated; indeed, that domain covers the observational data of a very A procedure determining such a canonical transformation C can be found in section 7 of [19]. Thecode allowing Mathematica to compute the symplectic matrix related to C is freely available at the webaddress ∼ locatell/MCSH/programmi/diag L4oL5.mth C. Caracciolo, U. Locatelli ̺ ̺ T ̺ ( ̺ ∗ ) ̺ T µ ≃ . T e . l . t . ≃ × .small fraction of Trojans. The results get even worse when all our rigorous estimatesabout both the remainder and the stability time are taken into account. As it is clearlyshown in the left hand side of Table 2, the construction of the non-resonant Birkhoffnormal form allows us to ensure the effective stability for initial conditions such that( I , I ) ∈ [0 , . × [0 , . ̺ . In fact, for this particular problem, it hasno physical meaning to prove stability for a time longer than the life-time of the Sun.Therefore, as it has been done in [24], we examine the values of ̺ for which the stabilitytime is comparable with the so called expected life-time T e . l . t . of the system. This numberis rescaled with respect to the revolution period of the celestial body, i.e., we impose that T e . l . t . = (6 Gyrs) · ν π ≃ × (in number of Jupiter revolutions, because ν is expressedin rad / yrs). As it has been highlighted in Figure 2, a small decrease of the value of theradii containing the initial condition ̺ can very significantly change the stability time(because of the exponential dependency of T on the inverse of ̺ ). Such a behaviour isevident also in Table 2, where with a change of the 1% of ̺ , the stability time is almostdoubled in the non–resonant case. Let us emphasize that in our estimates, we use thesame radius for both the actions, while looking at the initial conditions of the real trojanasteroids for Jupiter in [17], the value of the second action is usually smaller with respectto the first one. Therefore, we believe that considering the initial conditions of the actionsin polydisks of different radii for the actions, as in [24], could improve the results.In order to make easier the comparisons, in our opinion it is convenient to report theresults about both the non-resonant normal form and the resonant one next each other,as we have done in Table 2 and also in Figure 2. For the sake of definiteness, we have toexplain how to determine the radius ̺ of the open ball containing the initial conditions.Let us recall that at the end of subsection 4.1 we have chosen a value of ̺ in order tooptimize our evaluation of the lower bound of the escaping time T . Instead, here we focusour attention on the consistency of the procedure, in the sense that we aim to ensure thatthe motion is confined in the domain ∆ ̺ for all | t | ≤ T . Therefore, we fix ̺ ∈ (0 , ̺ ) asthe unique solution of the equation β̺ = ̺ + ν ( ̺ − ̺ ) ν , (48)where ν and ν have opposite signs in view of (47) and β < I (0) . This last concept definitely deserves a more detailed explanation. Indeed, inagreement with (44), we define ( ̺ ∗ ) = ̺ − ∆ I , with ∆ I given as in (43). Since the main omputer-assisted estimates for Birkhoff normal forms T ρ Jupiter: non resonant normal form T (ρ ∗2 ) Jupiter: resonant normal form
Figure 2: Plots of the evaluation of our lower bound of the escape time T (in semi-logscale). On the left, the graph is a function of ̺ , on the right, of ̺ ∗ . The horizontal linecorresponds to T e . l . t . = 5 × . See the text for more details.contribution to the maximal variation of the resonant action is due to the linear terms(see formulæ (41)–(43)), we have that ∆ I ≃ − ν ( ̺ − ̺ ) /ν and, therefore, ( ̺ ∗ ) ≃ β̺ .Let us rephrase a conclusion discussed at the end of section 4.2, by adapting it to thepresent context: we have that I ( t ) ∈ ∆ ̺ ∀ t ∈ [ − T, T ], provided that | I (0) | < ̺ and | I (0) | < ( ̺ ∗ ) ≃ β̺ . This ends the explanation on the meaning of β as the parameterruling the restriction on the set of values of I (0). For what concerns the choice of thevalue of β , we are interested in setting it very close to 1, in order to enlarge the domainof stability in I (0) as much as possible; on the other side, we have to take into accountoscillations of the value of the resonant action that are due to the angular dependence ofthe normal form part ¯ Z ( I , ϕ ) (see formulæ (41)–(42)). Indeed, this so called modulatingterm for the resonant action I is such that ¯ Z = O ( k I k ). This remark allows to imaginea further optimization procedure on the determination of β , that in principle would notbe so difficult, due to the already mentioned fact that the maximal variation ∆ I of theresonant action is linear with respect ot the the actions ( I , I ) in view of formulæ (41)–(43). For the sake of simplicity, we prefer to avoid such a further optimization procedure:in all the systems we have studied in the framework of the CPRTBP model, we have simplyset β = 0 . I (0) is not so crucial, because it cannot be greater than 10 %; therefore, itcannot substantially improve our results, as it is clearly shown in the summary of thecomparisons reported in the final Table 6.The right hand side of Table 2 includes the results based on the construction of aresonant Birkhoff normal form, whose performance is worse with respect to the non-resonant one. The difference of behaviour between the two methods can be explained,by looking at the top-left box in Figure 3, where there is the comparison between normsof the generating functions χ r for the resonant and non-resonant constructions, in thecase of the system having Sun and Jupiter as primary bodies. The increase of the normsproceeds more and less at the same rate up to the first 25 steps; afterwards, the norms ofthe generating functions related to the resonant normal form start to grow up faster than6 C. Caracciolo, U. Locatelli
10 20 30 40 50 60 70 80 || χ r || rJupiter
10 20 30 40 50 60 70 || χ r || rUranus
10 20 30 40 50 60 70 80 || χ r || rMars
10 20 30 40 50 60 || χ r || rJanus Figure 3: Growth of the norms (in semi-log scale) of the generating functions for thenon-resonant Birkhoff normal form (continuous line) and the resonant one (dashed line).From top to down and from left to right, the boxes refer to the cases of the systems havingSun–Jupiter, Sun–Uranus, Sun–Mars and Saturn–Janus as primary bodies, respectively.the other ones. On the one hand, while the algorithm for the resonant normal form isperformed, there are less terms where small divisors are introduced by the solution of thehomological equation (9); on the other hand, the normal form part is larger in the caseof the resonant construction and this clearly has an impact on the size of the generatingfunction. In fact, when the r –th normalization step is performed, new perturbative termsare introduced by the application of Lie series to the Hamiltonian H ( r − ; the main onesare of type L χ r Z k for k < r . Let us recall that in the non–resonant construction we have Z i = 0 for all odd values of index i ; therefore, in that case we have a lower number of newterms and the main new contribution is expected to be due to L χ r Z . This explains whyin the non-resonant case the growth of the generating functions has two different trendsfor odd and even indices, as it is clearly shown by the plots reported in Figure 3. On thecontrary, in the resonant case the behaviour of the norms of the generating functions ismore regular, at least up to a threshold value of r , for which the worsening effect due tothe terms generated by the normal form part becomes predominant with respect to theadvantage gained because of the better accumulation of small divisors. omputer-assisted estimates for Birkhoff normal forms ̺ ̺ T ̺ ( ̺ ∗ ) ̺ T µ ≃ . × − ) with T e . l . t . ≃ × . ̺ ̺ T ̺ ( ̺ ∗ ) ̺ T µ ≃ . × − ) with T e . l . t . ≃ × .Let us now focus on the comparison between the performances of the resonant andnon-resonant constructions, when they are made for different values of µ : the resultsare summarized in the top-right, bottom-left and bottom-right boxes of Figure 3 and inTables 3–5, which refer to the systems having Sun–Mars, Sun–Uranus and Saturn–Janusas primary bodies, respectively. These systems have been chosen because of two reasons:all of them host at least one trojan body and the corresponding values of the mass ratio µ allow a study where such a parameter spans rather regularly several orders of magnitude.In all the plots of Figure 3 (except the one referring to the case of the Jupiter Trojans)the size of the generating functions produced by the resonant construction algorithm isremarkably smaller with respect to the corresponding non-resonant ones. Therefore, theseries introduced by the resonant normalization procedure are convergent in a biggerneighbourhood of the origin; in other words, in such a situation we are able to iterate thecomputer-assisted estimates for larger values of ̺ .The results described in the present subsection are further summarized in Table 6,where the records are listed in decreasing order with respect to the mass ratio µ . Thisallows us to emphasize that the comparison between the computer-assisted estimatesbased on the non-resonant Birkhoff normal form and the resonant one are more and morein favour of the latter, when the value of µ tends to zero. Since the very beginning of our research project, one of our main goals was to complementthe Birkhoff normal form with a coherent scheme of computer–assisted estimates, in orderto provide rigorous evaluations of the effective stability time. This has been accomplished.Moreover, in our opinion we have developed a technique that could be adapted in sucha way to apply also in proximity of equilibrium points that are not purely elliptical (see,e.g., [30] and [29] to find examples of interesting problems that could be studied with suchan approach).We have successfully applied our procedure to two simple systems, which are close toelliptic equilibrium points such that the quadratic approximation of the Hamiltonian isnot convex (so preventing a trivial proof of stability using the energy as Lyapunov func-tion): the H´enon-Heiles model and the Circular Planar Restricted Three-Body Problem8
C. Caracciolo, U. Locatelli ̺ ̺ T ̺ ( ̺ ∗ ) ̺ T µ ≃ . × − ) with T e . l . t . ≃ × . µ ̺ (non − res . ) ( ̺ ∗ ) (reson . ) ( ̺ ∗ /̺ ) Jupiter 9 . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − ̺ and ( ̺ ∗ ) which refer to thestability domains for the non-resonant Birkhoff normal form and the resonant one, re-spectively. The results are reported as a function of different values of the mass ratio µ ,the name of the smaller primary in the corresponding CPRTBP model is reported in thefirst column.(CPRTBP). For the latter, we have easily adapted our approach to the construction of theBirkhoff normal form of resonant type and we have shown that the consequent computer-assisted estimates are in a better position when the mass ratio between the primary bodiesis very small. Nevertheless, the extent of the domain covered by our results is still farfrom being enough to explain the effective stability of most of the Trojans.In the near future, we plan to work for improving our computer-assisted results, forwhat concerns the applicability to realistic physical models, that are of interest, in par-ticular, in Celestial Mechanics and Astronomy. We think that in these fields there stillare problems that are open and fundamental , where rigorous proofs of stability can betackled by using an approach purely based on the normalization algorithm for KAMtori, which can ensure a perpetual topological confinement in models with two degrees offreedom. Nevertheless, there is a much wider range of possible applications to Hamilto-nian systems defined on phase spaces of higher dimensions. In such a framework, a verypromising strategy is based on the local construction of the Birkhoff normal form in theneighborhood of an invariant KAM torus: it has already started to provide some remark-able results (see, e.g., [22]). Here, in the context of the discussion about the perspectivesof our computational method, it is natural to describe the scaling properties of our ap-proach with respect to the number n , i.e., the degrees of freedom of the Hamiltonianmodel. It is well known that the remainder of the Birkhoff normal form is exponentiallysmall with respect to the distance ̺ form the equilibrium point; this allows to deducethe following best possible estimate about the stability time: T ∼ exp[( ˜ ̺/̺ ) /n ], that iscorresponding to the most non-resonant frequencies in the small oscillations limit, being˜ ̺ a positive constant (recall formulæ (12) and (38)). Such an asymptotic law is intrinsic For instance, see C. Caracciolo, U. Locatelli, M. Sansottera, M. Volpi: “Librational KAM tori inthe secular dynamics of the υ –Andromedae planetary system”, in preparation. A preliminary version isavailable on request to the authors. omputer-assisted estimates for Birkhoff normal forms ̺ the so called optimal normalization step r opt can largely exceedthe maximal polynomial degree R I + 2 of the Hamiltonian terms whose expansions areexplicitly computed; in such a condition, our method starts to iterate the estimates, i.e.,a procedure that is very weakly affected by the eventual increase of the number of degreesof freedom n . Nonetheless, the accuracy of the final estimates about the stability timeremarkably deteriorates when r opt becomes greater than R I (see Figure 1). Increasing n strongly reduces the maximal integer value, say R I;max that can be conveniently fixed forthe parameter R I ruling the size of the expansions stored on a computer. Since a polyno-mial function of maximal degree R I + 2 and depending on 2 n variables hosts a number ofterms that is O (cid:0) R n I (cid:1) , it is natural to expect the following scaling law for the thresholdcorresponding to the deterioration of our computational results: R I;max ∼ M / (2 n ) , being M a suitable positive constant.In the very particular case of the PCRTBP, there are other possible sources of furtherimprovements. Indeed, the normal form construction that has been designed in [42]adopted (since the very beginning) canonical coordinates which are much more suitablefor such a kind of Celestial Mechanics model. This is the reason why the approximationprovided by that normal form describes much more carefully the dynamics, when it iscompared with the numerical results based on the Birkhoff normal form. We think thatour computer-assisted estimates can be adapted to complement also the constructivealgorithm explained in [42]. This should provide results about the effective stabilityholding true in wider domains, hopefully covering an important fraction of the trojanasteroids. Acknowledgments
This work was partially supported by the “Mission Sustainability” programme of the Uni-versit`a degli Studi di Roma “Tor Vergata” through the project IDEAS (E81I18000060005)and by the “Progetto Giovani 2019” programme of the National Group of MathematicalPhysics (GNFM–INdAM) through the project “Low-dimensional Invariant Tori in FPU–like Lattices via Normal Forms”. The authors acknowledge the MIUR Excellence De-partment Project awarded to the Department of Mathematics of the University of Rome“Tor Vergata” (CUP E83C18000100006), in particular, because of the availability of thecomputational resources.
A Basics of validated numerics on a computer
In order to give a completely rigorous support to our computer-assisted proofs, in thepresent work we have implemented everywhere validated numerics in two complementaryways. First, interval arithmetic has been used to calculate the coefficients of the Taylorexpansions for the normal forms (as far as possible), while rigorous computations of upper[lower] bounds have been performed to estimate other quantities of interest, e.g., norms0
C. Caracciolo, U. Locatelli of functions [stability times, respectively]. The present appendix is devoted to summarizeour approach to validated numerics, that essentially follows what is described in section 3of [31].All our codes are written in C language and all the quantities of interest are computedusing the double type in C . The set of floating point numbers of double type that arerepresentable on a computer is defined as follows: R = (cid:8) (cid:9) ∪ ( ± X j =0 d j j ! · − : d j = 0 , ∀ j = 0 , . . . ) ∪ ( ± X j =1 d j j ! · σ : d j = 0 , ∀ j = 1 , . . . , σ ∈ Z , − ≤ σ ≤ ) , (49)where σ is the exponent and the digits d .d d . . . d make up the so called mantissa m ;they appear in the binary scientific notation x = ± m σ of each number x belonging tothe set R . Therefore, an overflow [ underflow ] situation can occur when a computationaloperation attempts to generate a number whose absolute value is greater [smaller] than M = (1 − − ) · [than m = 2 − and it is different from 0, respectively]. Letus generically use the symbol ∗ to refer to any elementary arithmetic operations +, − , · and / ; moreover, the result provided by a computer when it performs that same operationwill be denoted with ⊛ . According to such a setting, for instance, a ⊕ b is the computeroutput for the usual sum a + b , where both a and b belong to R . Therefore, the machineepsilon for double type floating point numbers in C (hereafter, simply ε ) is defined as ε = min { η ∈ R : 1 ⊕ η > } . By comparing the definition above with that in (49), one can easily realize that ε = 2 − and ε/ (cid:8) η ∈ R : 1 ⊖ η < (cid:9) .The 64 bit IEEE standard ensures that, if a ∗ b / ∈ R and m < | a ∗ b | < M , then either a ⊛ b = max (cid:8) R ∩ ( −∞ , a ∗ b ) } or a ⊛ b = min (cid:8) R ∩ ( a ∗ b, ∞ ) } ; furthermore, if a ∗ b ∈ R ,then a ⊛ b = a ∗ b . In words, these prescriptions can be summarized as follows: when ageneric elementary operation does not create a situation of either overflow or underflow ,then the result provided by a computer must be correct up to the last significant digit .In order to prevent the occurrence of overflows or underflows generated by an element-ary arithmetic operation ⊛ , we restrict to the “safe range” (according to definition 3.2in [31]), which is the following set of numbers: S = (cid:8) (cid:9) ∪ (cid:8) x ∈ R : 2 − ≤ | x | ≤ (cid:9) . (50) In the first row of formula (49), d is fixed to be zero, while d is set to 1 for the numbers appearingin the second row, according to the so called “hidden bit” rule. By taking into account that the integerexponents between − bits to be represented and 1 bit is needed for the sign, itshould be evident that the information to be stored for each double type number must be spread over8 bytes (=64 bits ). Also the square root √ enjoys such a peculiar property (see, e.g., section 3.1 of [45]). omputer-assisted estimates for Birkhoff normal forms ∀ a, b ∈ S : if a ⊛ b = 0 , then a ∗ b = 0 , else (cid:12)(cid:12)(cid:12)(cid:12) a ∗ ba ⊛ b − (cid:12)(cid:12)(cid:12)(cid:12) ≤ ε . There is an obvious exception to such a general rule: the division by zero. Of course,validated numerics is performed on a very restricted domain with respect to the realnumbers and it is not expected to extend operations that are meaningless in R . Therefore,the computational algorithm must be set in such a way that divisions by zero (and squareroots of negative numbers, etc.) are avoided.Looking at the way the mantissa is represented in the second row of formula (49),one can immediately realize that for that type of floating point numbers multiplying bythe factor 1 + ε [by 1 − ε ] is enough to increase [decrease, resp.] their absolute valueby the last significant digit. This last remark joined with the prescriptions of the 64 bitIEEE standard allows us to establish a few simple rules, in order to rigorously performinterval arithmetic. Let us introduce the following binary operators ⊛ + : S × S R and ⊛ − : S × S R , that are defined as follows: a ⊛ + b = (cid:26) ( a ⊛ b ) ⊙ (1 + ε ) if a ⊛ b ≥ , ( a ⊛ b ) ⊙ (1 − ε ) if a ⊛ b < a ⊛ − b = (cid:26) ( a ⊛ b ) ⊙ (1 − ε ) if a ⊛ b ≥ , ( a ⊛ b ) ⊙ (1 + ε ) if a ⊛ b < . (51)In view of the prescriptions of the 64 bit IEEE standard, for what concerns the true resultsof any arithmetic elementary operation we can conclude that a ⊛ − b ≤ a ∗ b ≤ a ⊛ + b ∀ a, b ∈ S , (52)where, again, ∗ generically denote +, − , · and / (with the only exception of the divisionby zero). We emphasize that the inequalities in formula (52) are fundamental, becausethey allow us to implement the rigorous computations of upper bounds (or lower ones) forquantities that have to be estimated. Moreover, the validated numerics can be applied tointerval arithmetic by suitably extending the binary operators introduced in (51). In fact,let ⊛ ± : S × S R be defined in such a way that [ a − , a + ] ⊛ ± [ b − , b + ] = [ c − , c + ],where c − = min (cid:8) a − ⊛ − b − , a − ⊛ − b + , a + ⊛ − b − , a + ⊛ − b + (cid:9) c + = max (cid:8) a − ⊛ + b − , a − ⊛ + b + , a + ⊛ + b − , a + ⊛ + b + (cid:9) . (53)Because of the definition of ⊛ ± and in view of formulæ (52)–(53), we can write thefollowing relation that is fundamental for the (rigorous) interval arithmetic: a ∗ b ∈ [ a − , a + ] ⊛ ± [ b − , b + ] ∀ a ∈ [ a − , a + ] , b ∈ [ b − , b + ] , (54)where, once again, the case with 0 ∈ [ b − , b + ] and ∗ representing the division must beexcluded. Let us stress that two (possibly redundant) list of elements appear after theminimum and the maximum in the definitions (53) of c − and c + , respectively, in order totake into account of the effects induced by the signs, when we are dealing with products2 C. Caracciolo, U. Locatelli and divisions. To fix the ideas, it is convenient to realize that, for instance, the rigorousextension of the sum to intervals, can be defined in a much shorter way, i.e., [ a − , a + ] ⊕ ± [ b − , b + ] = [ a − ⊕ − b − , a + ⊕ + b + ]. Since the algorithm constructing the Birkhoff normalform is based on Lie series (therefore, Poisson brackets), one can immediately realizethat the computations of the coefficients appearing in the expansions involve elementaryarithmetic operations, only. Thus, their rigorous extension to intervals is enough toperform validated numerics on the truncated expansions of the Hamiltonian as far as theycan be explicitly computed. Let us stress that the whole computational algorithm can beiterated provided that each generic arithmetic operation of type [ a − , a + ] ⊛ ± [ b − , b + ] =[ c − , c + ] generates a new result such that c − , c + ∈ S . This is the reason why, in ourcodes, we included tests, in order to verify that such a condition is always satisfied; if itis not so, the running of a program is immediately stopped.The situation concerning the rigorous bounds on the estimates is slightly more com-plicated, because, very quickly, it occurs a violation of the condition that all the majorantsof the norms are given by numbers belonging to the “safe range” set S . This is because theestimates can be iterated for many more normalization steps with respect to the explicitcalculation of the expansions, each of them requiring a huge occupation of the memory ifcompared with a single upper bound on the corresponding norm. Due to the dramaticallyfast growth of the norms of terms composing the series (that are asymptotically divergingfor the number of normalization steps going to infinity), the limit max S is usually tres-passed even for values of R II that are relatively small with respect to those considered inthe applications described in section 5. Let us recall that we are interested in iterating theestimates of the norms for a large number of normalization steps R II , in order to obtainbetter results. Therefore, it is convenient to represent the logarithm of the (positive valuesof the) upper bounds of the norms. For such a purpose, we have to introduce a functionlog + : S ∩ R + S such that log x . log + x for all positive x ∈ S . This can be done in arather obvious way, by adapting the approach described in section 3.1 of [45], in order tocombine validated numerics with a truncated series expansion of a logarithm. To fix theideas, let us limit to the case x ∈ S ∩ [1 , ∞ ), then we define log + x = 2 n ⊙ + log + ξ n , being( ξ n ⊖ + ∈ S ∩ [0 , .
01] such that ξ = x and ξ j = (cid:0)p ξ j − (cid:1) + ∀ j = 1 , . . . , n , where n isthe minimum nonnegative integer such that ξ n ≤ .
01 and the function (cid:0) √· (cid:1) + is nothingbut (1 + ε ) √· , in agreement with what has been explained in the corresponding footnote 9 .The rigorous upper bound for log (cid:0) ξ n ⊖ + ∀ ξ n ∈ S ∩ [1 , .
01] is introduced in asimilar way to the (simpler) definition of the function exp + ( · ) (see formula (55) below).The definition of log + x for x ∈ S ∩ (0 ,
1) is analogous to the case discussed just abovewith x ≥ D r = ( r + 2) G r , with G r = F ( r − r /α r . In our codes we can write the correspondingrigorous estimate as log D r = log + ( r + 2) ⊕ + log F ( r − r ⊕ + log + (cid:0) ⊘ + α r (cid:1) , where we Let us stress that in order to ensure that log x ≤ log + x where log + x = 2 n ⊙ + log + ξ n and providedthat log ξ n ≤ log + ξ n ∀ ξ n ∈ S ∩ [1 , . S . omputer-assisted estimates for Birkhoff normal forms F ( r − r is a previously defined number belonging to the “safe range” set S . When algebraic sums are involved, the procedure is slightly more complicate. As afurther example, let us focus on (31): in order to properly define an upper bound for thelogarithm of a r = (cid:0) a rr − + ( r + 1) D r (cid:1) /r we have to compute log + ( x + y ), where the valuesof log x = r ⊙ + log a r − and log y = log + ( r + 1) ⊕ + log D r have to be considered as known,because both log a r − and log D r have been preliminarily estimated by some rigorouscomputations; moreover, let us recall that we want to avoid to compute x = exp(log x ) and y = exp(log y ), because they are expected to be too large numbers, eventually exceedingmax S . Without any loss of generality, let us assume that x ≥ y , therefore, it is convenientto set log + ( x + y ) = log x ⊕ + log + (cid:16) ⊕ + exp + (cid:0) log y ⊖ + log x (cid:1)(cid:17) . Thus, we are lead to the problem of defining a function exp + : S ∩ [0 , S suchthat exp( x ) . exp + ( x ) ∀ x ∈ S ∩ [0 , ξ j = x ⊘ + j ∀ j = 0 , . . . , n , where n is the minimum positive integer such that ξ n ∈ S ∩ [0 , . + ( ξ n ) = (cid:18) N L + i =0 ( π i ⊘ + i !) (cid:19) ⊕ + (cid:16)(cid:0) π N ⊘ + N ! (cid:1) ⊘ + (1 ⊖ − ξ n ) (cid:17) , (55)where N L + i =0 ( π i ⊘ + i !) = 1 ⊕ + . . . ⊕ + ( π N ⊘ + N !) with π i = π i − ⊙ + ξ n ∀ i = 1 , . . . , N (being π = 1); moreover, N is the minimum positive integer such that π N ⊘ + N ! ≤ ε ,while P i> N ξ in /i ! ≤ (cid:0) π N ⊘ + N ! (cid:1) ⊘ + (1 ⊖ − ξ n ) is nothing but the estimate of the trun-cated remainder for the expansion in Taylor series of exp( ξ n ). The rigorous computationof exp + ( x ) = exp + ( ξ ) is completed by recursively defining exp + ( ξ j ) = exp + ( ξ j +1 ) ⊙ + exp + ( ξ j +1 ) in such a way to proceed backwards from j = N − j = 0.The whole of the elementary operations described in the present appendix representthe bare minimum that is sufficient to implement validated numerics, in order to makefully rigorous our computer-assisted proofs. B A complete example of application to the H´enon-Heiles model in a case with short expansions
The aim of this appendix is very pedagogical: we explain step-by-step an application ofthe algorithm described in the sections 2–4. For the sake of clarity, here we focus again onthe H´enon–Heiles model (that is the simplest one among those considered in section 5),by performing a very low number of normalization steps. In fact, the example describedin the present appendix deals with the case where R I = 2 and R II = 5; this means thatthe expansions of larger polynomial degree to be explicitly computed are quartic whilethe iteration of the estimates provides upper bounds for terms up to the seventh degree. These functions are included in the lib val num.c library, that makes part of the package freelyavailable at the web address ∼ locatell/CAPs/BirkCAP AppsA-B.zip C. Caracciolo, U. Locatelli
In our opinion this choice is a good balance: on the one hand it is non-trivial, on theother hand the representation of the Hamiltonians is still readable because it is not toolarge. Moreover, we have decided to further simplify our example, by omitting the studyof our better estimate for the stability time as function of the radius ̺ , whose value isfixed to ̺ = 0 . H ( r − (the Taylor expansion of which is explicitly given in equation (7)) is represented by a set gathering both polynomial functions and numbers: S ( r − = (cid:26) Z , . . . , Z min { r − ,R I } , f ( r − { r,R I } , . . . , f ( r − R I , log Z R I +1 , . . . , log Z min { r − , R II } , log F ( r − { r,R I +1 } , . . . , log F ( r − R II , log E , log a r − (cid:27) . (56)The polynomials appearing in the first row of definition (56) are written in terms of theirexpansions, that are Z s = P | ℓ | = s c ( s ) ℓ , ℓ ; ± ( − i z ) ℓ ¯ z ℓ and f ( r − s = P | ℓ | + | ˜ ℓ | = s c ( r − ,s ) ℓ , ˜ ℓ ; ± ( − i z ) ℓ ¯ z ˜ ℓ ,where the coefficients are complex numbers expressed as intervals, that are explicitly com-puted with the help of an algebraic manipulator . Moreover, the values Z s and F ( r ) s arethe rigorous upper bounds for the norms of the corresponding Hamiltonian polynomialswhose expansions are not explicitly computed, while E and a r are the parameters neededto estimate the infinite tail of terms of type f ( r − s with s > R II . Let us recall that the val-ues appearing in the second and third rows of definition (56) are expressed in logarithmicscale because some of them can eventually exceed the “safe range” set S of representablenumbers on a computer, as it has been widely discussed in appendix A.In our opinion, the main concept to be kept in mind can be summarized as follows.The prescriptions included in sections 2–3 allow to perform the r –th normalization step,in such a way to explicitly determine all the elements appearing in S ( r ) , which represents H ( r ) and is a set having finite cardinality with the same structure as that described in (56)(where r − r ). As a short reformulation, this means that the Let us stress that in formula (56) some unpleasant lower indexes (that are involving either the min-imum or the maximum between nonnegative integers) are somehow unavoidable, because the number ofalready performed normalization steps, i.e. r −
1, can be greater than R I or not; therefore, either the ex-plicit expansions of f ( r − { r,R I } , . . . , f ( r − R I are missing or the majorants log Z R I +1 , . . . , log Z min { r − , R II } disappear in the list of elements making part of the set S ( r − . It is not easy to code with X ̺ ´o ν o ζ , which is the algebraic manipulator we actually used for allthe applications discussed in the present paper, because its syntax looks quite difficult for new users.Therefore, we think that an implementation of the algorithm constructing the Birkhoff normal form inthe framework provided by Mathematica can be more helpful for a reader that is interested in repro-ducing our results. We have written such a code (named birkhoff HH.mth ) in a version computingexactly the coefficients as algebraic numbers. It makes part of the package freely available at the webaddress ∼ locatell/CAPs/BirkCAP AppsA-B.zip and its output is eas-ily converted in rigorous upper bounds on the norms by another program (named estimates.c ). Alsothis further code is included in that same software package and is in charge to compute rigorously all theupper and lower bounds that are needed, in order to conclude the computer-assisted proof of the sameresult discussed in the present appendix. omputer-assisted estimates for Birkhoff normal forms S ( r − into S ( r ) and now we are going to do it from S (0) to S (5) . Input.
We focus on the Hamiltonian of the H´enon–Heiles model (45) with ω = 1 and ω = −√ /
2. After having rewritten it in complex canonical variables ( − i z , ¯ z ) (see (6)and (3) for the definition), the set S (0) representing the initial Hamiltonian H (0) can bewritten as follows: S (0) = n Z , f (0)1 , , − , − , − , , . o , (57)where Z = i . ± ( − iz )¯ z − i . ± ( − iz )¯ z and f (0)1 = − i . ± ( − iz ) ( − iz ) − . ± ¯ z − . ± ( − iz ) ¯ z − . ± ( − iz )( − iz )¯ z + i . ± ( − iz )¯ z ¯ z + i . ± ( − iz ) +0 . ± ( − iz ) ¯ z + i . ± ( − iz )¯ z − i . ± ( − iz )¯ z + 0 . ± ¯ z ¯ z . In both the expansions above we adopted the notation c ± = a ± σ a + i b ± σ b for each coeffi-cient, where a and b are the central values referring to the intervals of the real and ima-ginary part, while σ a and σ b allow to determine the half-widths σ a × − e a and σ b × − e b of those same intervals, being e a and e b the number of digits appearing after the float-ing point in the writing of a and b , respectively. For instance, the last summand in theexpansion of Z reads as i [ − . , − . − iz )¯ z . We remark that f (0) s = 0 ∀ s ≥
2, because the H´enon–Heiles model is defined by a cubicHamiltonian. This explains why the third element of S (0) is equal to zero. Since we havedecided to describe the upper bounds of the norms in a logarithmic scale (in order tomaintain the agreement with more challenging applications involving huge expansions),therefore, we have put log( F (0) s ) = − for 3 ≤ s ≤ e − ≃ − . .For what concerns the last two parameters E and log a appearing at the end of therepresentation S ( r − (56) in the case with r = 1, usually they are determined taking intoaccount both the sup norm and the geometrical decay of the power series in a domainwhere the Hamiltonian H (0) is an analytic function. In this particular case, any choice of a We have made such a choice for the values of the angular velocities ω and ω , because all thecoefficients of the expansions defined by the normalization algorithm can be represented in the form( l + m √ /n , being l, m, n ∈ Z . This makes easier the exchange of information between the codes birkhoff HH.mth and estimates.c as it has been described in the previous footnote 13 . C. Caracciolo, U. Locatelli pair of positive values for E and log a would be acceptable. However, since the value of thecommon ratio a r − (that is related to the majorants of the infinite tail of terms appearingin the remainder of H ( r − ) affects the next one, i.e., a r , through the relation (31), we havefound reasonable to set E and a in such a way to describe the relative increase of f (0)1 withrespect to Z . This is the reason why the last two elements making part of S (0) in (57)have been fixed so that log E = 0 and log a = log F (0)1 ⊘ +
3, where log F (0)1 = log + ( k f (0)1 k )and the definitions of the elementary arithmetic operators allowing to implement validatednumerics (like, e.g., ⊘ + ) are given in the previous appendix A. First normalization step: r = 1 . We have to put in Birkhoff normal form the cubicterm. In order to do that, we solve the homological equation L χ Z + f (0)1 = Z and, inview of equation (11), we obtain χ = 0 . ± ( − iz ) ( − iz ) − i . ± ( − iz ) ¯ z + i . ± ( − iz )( − iz )¯ z − i . ± ¯ z − . ± ( − iz )¯ z ¯ z + 0 . ± ( − iz ) − i . ± ( − iz ) ¯ z + 0 . ± ( − iz )¯ z +0 . ± ( − iz )¯ z − i . ± ¯ z ¯ z . The new normal form term Z is equal to 0, due to the fact that the normalization step isodd (see (13)). We have now to apply the Lie series of generating function χ to H (0) ; sincewe are considering terms up to degree 4, the only contributions we have to calculate are L χ Z and L χ f (0)1 ; in fact, according to formula (14), they will produce terms of degree4, which are collected together in the definition of f (1)2 = L χ Z + L χ f (0)1 . In order toproperly describe all the Hamiltonian terms at the end of the first normalization step, wehave to update the parameters appearing in the second and in the third row of the genericrepresentation described in (56). Since we know the expansion of the generating function χ , we can use the formulæ in Proposition 3.3 to rigorously compute the upper bounds D and log( F (1) s ) for 3 ≤ s ≤
5; moreover, we can compute the norm of f (1)2 , that islog + ( k f (1)2 k ) = 2 . a = 2 . H (1) can be summarized by the following set: S (1) = n Z , , f (1)2 , . , . , . , , . o . Finally, we use formula (29) to determine the rigorous estimate of the sup norm of theremainder after having completed the first normalization step: log( |R (1) | ̺ ) = − . Second normalization step: r = 2 . For what concerns the second step, we proceedin an analogous way with respect to the first one. The homological equation we have tosolve is L χ Z + f (1)2 = Z . We compute the generating function χ by using formula (11)and the new normal form term Z as defined in (13). All the other terms generated by For brevity reasons, here we have decided to not include the expansion of f (1)2 . However, its normalform part, i.e., Z , is written within the following description of the second normalization step. omputer-assisted estimates for Birkhoff normal forms S (2) = { Z , , Z , . , . , . , , . } , where Z = − . ± ( − iz ) ¯ z − . ± ( − iz ) ¯ z +1 . ± ( − iz )( − iz )¯ z ¯ z . As a new estimate for the remainder, we obtain log( |R (2) | ̺ ) = − . Third normalization step: r = 3 . As a major difference with respect to the secondstep, now an explicit expansion for f (2)3 is not available and, therefore, the same holdsfor the generating function χ . Instead of computing D by using the expansion of χ ,we have to set D = 5 G (see formula (26) and the definition of D r in Proposition 3.3).Indeed, G can be easily computed because F (2)3 is known and α r can be determined,because it is the smallest divisor which could appear at the third step (see (25)). Afterhaving determined the upper bounds Z and log( F (3) s ) for s = 4 , S (3) = { Z , , Z , . , . , . , , . } . The new estimate of the remainder is log( |R (3) | ̺ ) = − . Fourth normalization step: r = 4 . We can proceed exactly in the same way as forthe third step. The representative set corresponding the new Hamiltonian H (4) is S (4) = { Z , , Z , . , . , . , , . } . The new remainder can be estimated as log( |R (4) | ̺ ) = − . End of the algorithm.
We reconsider the set S (3) because it represents the Hamilto-nian at the third normalization step, which is related to the smallest remainder. Thevalue of ̺ that optimizes the escaping time, as defined at the end of subsection 4.1, is ̺ = 0 . T byusing formula (37) and we obtain log T = 24 . C. Caracciolo, U. Locatelli
Final comments.
It is easy to realize that the explicit computation of the expansionsfor the terms belonging to the Hamiltonians up to the degree R I + 2 does not depend onthe evaluation of the upper bounds for higher order polynomials. This remark explainsthe reason why it is convenient to separate the computer-assisted proof in two differentprograms, that are designed in order to handle with the expansions and the estimates,respectively. These two codes are included in BirkCAP AppsA-B.zip and they are in aversion that is adapted to the example discussed in the present appendix. Such a fileis conceived as a sort of supplementary material with respect to this work, thus, itcan be interesting for readers willing to reproduce the computational algorithm in theirown codes (or to adapt our ones). For such a purpose, in BirkCAP AppsA-B.zip wehave included a program computing exactly the polynomial expansions in the frameworkof
Mathematica , in order to provide a code that is rather easy to read. On the otherhand, the explicit expansions reported in this appendix have been produced by using X ̺ ´o ν o ζ for the algebraic manipulations of polynomials with coefficients that are complexnumbers expressed as intervals. We think that this can be useful for the initial comparisonswith a new code, eventually designed by a reader interested in much more challengingapplications with respect to the very simple example discussed in the present appendix. References [1] K. Appel and W. Haken:
Every planar map is four colorable. Part I. Discharging ,Illinois J. Math., (1977), 429–490.[2] K. Appel and W. Haken: Every Planar Map Is Four Colorable , A.M.S. Contemp.Math., (1989).[3] K. Appel, W. Haken and J. Koch: Every planar map is four colorable. Part II.Reducibility , Illinois J. Math., (1977), 491–567.[4] V.I. Arnold: Proof of a theorem of A. N. Kolmogorov on the invariance of quasi–periodic motions under small perturbations of the Hamiltonian , Usp. Mat. Nauk, (1963), 13; Russ. Math. Surv., (1963), 9.[5] G.D. Birkhoff: Dynamical systems , New York (1927).[6] I. Bal´azs, J. Bouwe van den Berg, J. Courtois, J. Dud´as, J.-P. Lessard, A. V¨or¨os-Kiss, J.F. Williams and X.Y. Yin:
Computer-assisted proofs for radially symmetricsolutions of PDEs , J. Comp. Dyn., (2018), 61–80.[7] A. Celletti and L. Chierchia: KAM stability and Celestial Mechanics , Memoirs ofAMS, (2007), n. 878.[8] A. Celletti and A. Giorgilli:
On the stability of the Lagrangian points in the spatialrestricted problem of three bodies , Cel. Mech. & Dyn. Astr., (1991), 31–38. It is freely available at the web address already reported in the previous footnotes 11 ,
13 . omputer-assisted estimates for Birkhoff normal forms
Improved estimates on the existence ofinvariant tori for Hamiltonian systems , Nonlinearity, (2000), 397–412.[10] T. M. Cherry: On integrals developable about a singular point of a Hamiltoniansystem of differential equations , Proc. Camb. Phil. Soc. (1924), 325–349.[11] T. M. Cherry: On integrals developable about a singular point of a Hamiltoniansystem of differential equations, II , Proc. Camb. Phil. Soc. (1924), 510–533.[12] G. Contopoulos: A Review of the “Third” Integral , Mathematics in Engineering, (2020), 472–511.[13] C. Efthymiopoulos, A. Giorgilli and G. Contopoulos: Nonconvergence of formal in-tegrals: II. Improved estimates for the optimal order of truncation , J. Phys. A: Math.Gen., (2004), 10831–10858.[14] C. Efthymiopoulos and Z. S´andor: Optimized Nekhoroshev stability estimates for theTrojan asteroids with a symplectic mapping model of co-orbital motion , Mon. Not.R. Astron. Soc., (2005), 253-271.[15] J.-Ll. Figueras, A. Haro and A. Luque:
Rigorous computer-assisted application ofKAM theory: a modern approach , Found. Comput. Math., (2017), 1123–1193.[16] F. Gabern and A. Jorba: A restricted four-body model for the dynamics near theLagrangian points of the Sun-Jupiter system , Discrete Contin. Dyn. Syst. Ser. B, (2001), 143–182.[17] F. Gabern, A. Jorba and U. Locatelli: On the construction of the Kolmogorovnormal form for the Trojan asteroids , Nonlinearity, (2005), n.4, 1705–1734.[18] A. Giorgilli: Notes on exponential stability of Hamiltonian systems , “DynamicalSystems, Part I: Hamiltonian systems and Celestial Mechanics”, Pubblicazioni delCentro di Ricerca Matematica Ennio De Giorgi, Pisa (2003), 87–198.[19] A. Giorgilli, A. Delshams, E. Fontich and C. Sim´o:
Effective stability for Hamiltoniansystem near an elliptic equilibrium point, with an application to the restricted threebody problem , Journ. of Diff. Eq., (1989), 167-198.[20] A. Giorgilli and U. Locatelli: Canonical perturbation theory for nearly integrable sys-tems , in: B.A. Steves, A.J. Maciejewski, M. Hendry eds.,
Chaotic Worlds: FromOrder to Disorder in Gravitational N-Body Dynamical Systems , Nato Sc. series,Springer, (2006).[21] A. Giorgilli, U. Locatelli and M. Sansottera:
Kolmogorov and Nekhoroshev theory forthe problem of three bodies , Cel. Mech. & Dyn. Astr., (2009), 159–173.[22] A. Giorgilli, U. Locatelli and M. Sansottera:
Secular dynamics of a planar model ofthe Sun-Jupiter-Saturn-Uranus system; effective stability in the light of Kolmogorovand Nekhoroshev theories , Regular and Chaotic Dynamics, (2017), 54–77.0 C. Caracciolo, U. Locatelli [23] A. Giorgilli and M. Sansottera:
Methods of algebraic manipulation in perturbationtheory , in “Chaos, Diffusion and Non-integrability in Hamiltonian Systems – Applic-ations to Astronomy”, Proceedings of the Third La Plata International School onAstronomy and Geophysics, P.M. Cincotta, C.M. Giordano and C. Efthymiopouloseds., Universidad Nacional de La Plata and Asociaci´on Argentina de Astronom´ıaPublishers, La Plata (2012).[24] A. Giorgilli and Ch. Skokos:
On the stability of the Trojan asteroids , Astron. As-troph., (1997), 254–261.[25] W. Gr¨obner and H. Knapp:
Contributions to the Method of Lie-Series , Bibliograph-isches Institut, Mannheim (1967).[26] F.G. Gustavson:
On constructing formal integrals of a Hamiltonian system near anequilibrium point , Astron. J., (1966), 670–686.[27] M. H´enon and C. Heiles: The applicability of the third integral of motion: somenumerical experiments , The Astronomical Journal, (1964), 73–79[28] M. H´enon: Exploration num´erique du probl`eme restreint IV: Masses ´egales, orbitesnon p´eriodiques , Bulletin Astronomique, (1966), 49–66.[29] T. Johnson and W. Tucker Automated computation of robust normal forms of planaranalytic vector fields , Disc. & Cont. Dyn. Syst. – Series B, (2009), 769–782.[30] `A. Jorba and J. Masdemont: Dynamics in the centre manifold of the collinear pointsof the Restricted Three Body Problem , Physica D, (1999), 189–213.[31] H. Koch, A. Schenkel and P. Wittwer:
Computer-assisted proofs in analysis andprogramming in logic: a case study , SIAM Review, (1996), 565–604.[32] A.N. Kolmogorov: Preservation of conditionally periodic movements with smallchange in the Hamilton function , Dokl. Akad. Nauk SSSR, (1954), 527. Engl.transl. in: Los Alamos Scientific Laboratory translation LA-TR-71-67; reprinted in:Lecture Notes in Physics, .[33] C. Lhotka, C. Efthymiopoulos and R. Dvorak: Nekhoroshev stability at L4 or L5in the elliptic–restricted three–body problem — application to the Trojan asteroids ,Mon. Not. R. Astron. Soc., (2008), 1165–1177.[34] J.E. Littlewood:
On the equilateral configuration in the restricted problem of threebodies , Proc. London Math. Soc. (3) (1959), 343–372.[35] J.E. Littlewood: The Lagrange configuration in celestial mechanics , Proc. LondonMath. Soc. (3) (1959), 525–543.[36] R.S. MacKay and J. Stark: Locally most robust circles and boundary circles for area–preserving maps , Nonlinearity, (1992), 867–888. omputer-assisted estimates for Birkhoff normal forms On invariant curves of area–preserving mappings of an annulus , Nachr.Akad. Wiss. G¨ott., Math. Phys., (1962), 1–20.[38] N.N. Nekhoroshev: Exponential estimates of the stability time of near–integrableHamiltonian systems.
Russ. Math. Surveys, (1977), 1.[39] N.N. Nekhoroshev: Exponential estimates of the stability time of near–integrableHamiltonian systems, 2.
Trudy Sem. Petrovs., (1979), 5.[40] H. Poincar´e: Les m´ethodes nouvelles de la m´ecanique c´eleste , Gauthier–Villars, Paris(1892).[41] M. Sansottera, U. Locatelli and A. Giorgilli:
On the stability of the secular evolutionof the planar Sun-Jupiter-Saturn-Uranus system , Math. Comp. Sim., (2013), 1–14.[42] R.I. P´aez and U. Locatelli: Trojans dynamics well approximated by a new Hamilto-nian normal form , Mon. Not. R. Astron. Soc., (2015), 2177–2188 .[43] M. Sansottera, A. Giorgilli and T. Carletti:
High-order control for symplectic maps ,Physica-D, (2016), 1–15.[44] M. Sansottera, C. Lhotka and A. Lemaˆıtre:
Effective stability around the Cassinistate in the spin-orbit problem , Cel. Mech. & Dyn. Astr., (2014), 75–89.[45] A. Schenkel, J. Wehr and P. Wittwer:
Computer-assisted proofs for fixed point prob-lems in Sobolev spaces , MPEJ, , n. 3 (2000), 1–67.[46] Ch. Skokos and A. Dokoumetzidis: Effective stability of the Trojan asteroids , Astron.Astroph., (2001), 729–736.[47] V. Szebehely:
Theory of Orbits , Academic Press, New York (1967).[48] E. T. Whittaker:
On the adelphic integral of the differential equations of dynamics ,Proc. Roy Soc. Edinburgh, Sect. A,37