Computer-assisted estimates for Birkhoff normal forms

Birkhoff normal forms are commonly used in order to ensure the so called"effective stability"in the neighborhood of elliptic equilibrium points for Hamiltonian systems. From a theoretical point of view, this means that the eventual diffusion can be bounded for time intervals that are exponentially large with respect to the inverse 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 a rather classical scheme of estimates for both the Birkhoff normal forms to any finite order and their remainders. This is made for providing explicit lower bounds of the stability time (that are valid for initial conditions in a fixed open ball), by using a fully rigorous computer-assisted procedure. We apply our approach in two simple contexts that are widely studied in Celestial Mechanics: the H\'enon-Heiles model and the Circular Planar Restricted Three-Body Problem. In the latter case, we adapt our scheme of estimates for covering also the case of resonant Birkhoff normal forms and, in some concrete models about the motion of the Trojan asteroids, we show that it can be more advantageous with respect to the usual non-resonant ones.


Introduction
The Birkhoff normal form represents a milestone in the historical development of the Hamiltonian perturbation theory, because it has been designed in the simplest possible way so as to (at least) partially escape the results obtained by Poincaré for what concerns both the non-existence of the first integrals and the ubiquity of chaos (see [40]). Indeed, it was introduced about one century ago (see [48], [10]- [11] and [5]) in the framework of the dynamics in the neighborhood of an elliptic equilibrium point; in section 2 such a context is recalled by a set of proper definitions. More recently, the scheme of analytical estimates has been further improved in order to ensure that the stability time is exponentially big with respect to (some fractional power of) the inverse of the distance from the equilibrium point (see [34]- [35], [38]- [39] and section 4 of [20] for a brief summary of such a kind of results).
The approach based on suitable estimates for Birkhoff normal forms has allowed to introduce the concept of effective stability for physical systems, that can be stated as follows. Let us focus on the motions starting from initial conditions staying in a ball of radius ̺ 0 , that is centered at an equilibrium point; we say that it is effectively stable, if we can ensure that the solutions of the Hamilton equations stay at a distance not larger than ̺ > ̺ 0 for an interval of time T that is exceeding the (expected) life-time of that system. This kind of approach has been successfully applied to the so called Circular Planar Restricted Three-Body Problem (hereafter, CPRTBP), where the vicinity of the triangular Lagrangian points (L 4 and L 5 ) in a model including Sun and Jupiter as primary bodies is considered. In [24], it is shown that the orbits of the third massless body stay within a ball B ̺ (0) having radius ̺ and centered in L 4 or L 5 for times larger than the age of the universe, where the value of ̺ is approximately the same as ̺ 0 and the ball B ̺ 0 (0) is large enough to cover the initial conditions corresponding to the observations of four Trojan Asteroids. Such a result has been extended in several ways, with applications to problems arising, mainly, in Celestial Mechanics; the following list is far from being exhaustive. In [46], a Trojan Asteroid is shown to be effectively stable in the 3D extension of the model described just above. In [16], the effects induced by Saturn are taken into account indirectly, because the problem is restricted in such a way that the orbit of Jupiter is prescribed so as to be both periodic and a good approximation of a numerical solution for 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 describing a CPRTBP and an Elliptic Restricted Three-Body Problem, respectively. In [44], the effective stability of the rotational motion of Mercury is shown in the framework of a Hamiltonian model. Last but not least, in [43] Birkhoff normal forms are used in order to efficiently calculate counter-terms, with the purpose of controlling the diffusion in a symplectic map of Hénon type. This has been made also in view of further possible extensions to similar models describing the dynamics of particle accelerators.
In spite of its paramount importance, in none of the previously mentioned articles the 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 are based on the explicit calculation of huge numbers of coefficients appearing in the Birkhoff normal forms, made by using software packages that are suitable for performing algebraic manipulations. On the other hand, the computations do not implement interval arithmetic; moreover, for what concerns the series expansions appearing everywhere in the algorithms, the analytic estimates of both their radii of convergence and the sup-norms are missing. The numerical evaluation of such quantities is often made by using simplistic arguments. In this sense, those results are not rigorous proofs. However, completely rigorous computer-assisted estimates have been carried out in a paper that is focused on CPRTBP and is an ancestor of a few of those mentioned above, but its results were strongly 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 allowed to increase very remarkably the threshold of applicability with respect to the small parameter, that is a crucial problem in KAM theory (see [32], [37] and [4]). In fact, while at the very beginning the first purely analytical results were so limited that they were completely inadequate for any realistic application, recent rigorous computer-assisted estimates have been able to get very close to the optimal breakdown threshold on the small parameter, 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 significantly 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 Uranus is considered). For what concerns the very specific case of the CPRTBP, [17] includes rigorous computer-assisted proofs of existence of KAM tori in the vicinity of three Trojan Asteroids. Moreover, in that same work it is shown the numerical evidence that the algorithm constructing the Kolmogorov normal form is successful in 23 cases over the first 34 Trojan Asteroids listed in the IAU Catalogue. Let us recall that more than 7 000 asteroids have been detected about the triangular Lagrangian points of the system having Sun and Jupiter as primary bodies. Such a large number of observed objects, compared with the rather small one for which the previously described results are available, highlights the need for improvements of the mathematical framework, in order to fully explain the orbital stability of these celestial bodies.
Our work aims to contribute in filling some of the gaps that have been discussed above. In particular, we want to settle a completely rigorous scheme of computer-assisted estimates for Birkhoff normal forms. Moreover, we want to test its performances. In view of such a purpose, we decided to focus just on a pair of Hamiltonian systems, that are fundamental examples and have two degrees of freedom. First, we apply our method to the Hénon-Heiles model (see [27] and the very recent review in [12]), that is a very simple example of a perturbed couple of harmonic oscillators. In particular, we study the case with two non-resonant frequencies having opposite signs, because this implies the failure of any trivial proof scheme of the stability. Indeed, the Hamiltonian is not a Lyapunov function 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 in the CPRTBP. Let us recall that the structure of the Hamiltonian for the latter model is very similar to that of the former one; in particular, trivial approaches to the problem of stability are hopeless for both. For what concerns the CPRTBP, we study a variety of situations by considering a few sub-cases where the role of the second body is played by some different planets. This allows us to investigate the problem also with very small mass 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 the expansions of the formal integrals of motion for resonant systems go back at the dawn of the computer age and refer to the Hénon-Heiles model (see, e.g., [26]).
While in KAM theory there are powerful computer-assisted techniques that are not based on the Hamiltonian formalism for canonical transformations (see [7] and [15]), we emphasize that the results about effective stability are intrinsically related with normal forms. 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 of analytical estimates (see [18] for an introduction).
The paper is organised as follows. In section 2, we describe the construction of the Birkhoff normal form by an approach based on Lie series. This is the starting point for the definition of the scheme of estimates, which is explained in section 3. In section 4, we discuss how to use the estimates to prove effective stability for Birkhoff normal forms in both the non-resonant case and the resonant one. In section 5, we describe our results for the Hénon-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 are essential to implement our computer-assisted proofs in a fully rigorous way. In appendix B, we describe in detail an application to the Hénon-Heiles model in a situation, that is tailored in such a way that a relatively small number of computations are needed to obtain a (very limited) result; this small tutorial example is discussed with the aim of making our work more easy to reproduce for a reader that is interested in exporting our approach in other challenging contexts.

Construction of Birkhoff normal form around an elliptic equilibrium point
In this section we are going to describe the construction of the Birkhoff normal form for a particular class of Hamiltonians. This is made by using Lie series operators. Although such a method is rather standard in perturbation theory, a detailed description of the formal algorithm is mandatory, in order to make well definite the discussion in the next sections. The adoption of an approach based on Lie series makes our work further different with respect to those that have been mentioned in the Introduction and deal with effective stability estimates based on Birkhoff normal forms, because they implemented algorithms using 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 (·) = {f, ·}, being {·, ·} the Poisson bracket 1 between two functions. Moreover, the Lie series operator having f as generating function is defined as follows: and it will be used to define near to the identity canonical transformations. Moreover, we denote 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 polynomials with respect to the Poisson brackets.
This simple property is fundamental, because it will be repeatedly used for reorganizing the Hamiltonian terms (appearing in an expansion) with respect to their polynomial degree, after each canonical change of coordinates defined by a Lie series operator.

Hamiltonian framework
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: In the proximity of the origin, the magnitude order of the polynomial terms f ℓ is obviously smaller than that of the main part of the Hamiltonian, which is given by a sum of harmonic oscillators. Moreover, if we introduce the action-angle coordinates for harmonic oscillators, by the following canonical change of variables then the Hamiltonian in (2) transforms to with f ℓ containing terms having total degree in the square root of the actions I 1 , . . . , I n equal to ℓ and trigonometric expansions in ϕ with Fourier harmonics k such that |k| = j |k j | ≤ ℓ. These variables highlight the integrability of the quadratic part, which depends on the actions only. Therefore, we will proceed in a perturbative way, by leading the Hamiltonian to a Birkhoff normal form. This means that, after having performed r canonical changes of coordinates defined by the corresponding Lie series, the Hamiltonian will be such that where the upper index r in every Hamiltonian term refers to the number of normalization steps that have been already performed. In some more words, we want to remove step by step the angular dependence from the perturbative part of the Hamiltonian, in order to decrease the size of the remainder R (r) . This is done to provide a better integrable approximation Z (r) with respect to the problem we are studying.
In order to deal with the construction of the normal form and the subsequent scheme of estimates, it is convenient to work with the complex canonical variables (−iz,z), introduced by the following canonical transformation In these new variables I j = z jzj , then a new normal form term should contain only monomials where z j andz j appear with the same exponent. Moreover, since the polynomial degree is the same if we are using either the canonical variables (y, x) or the complex ones (−iz,z), we will maintain the symbol P ℓ , ∀ ℓ ≥ 0, to denote every class of homogeneous polynomials.

Description of the r-th normalization step
We are going to describe in detail how to perform a single step of the normalization algorithm. Let us suppose that we have been able to perform the first r − 1 steps, in such a way that the Hamiltonian is given the following form: where Z ℓ ∈ P ℓ+2 is depending on the actions only and f (r−1) ℓ ∈ P ℓ+2 , while the upper index r − 1 (appearing also in the symbol denoting the Hamiltonian) refers to the normalization step. When r = 1, this is nothing but Hamiltonian (4) once it has been expressed in complex variables. We want to normalize the main term among the perturbative ones, that is f (r−1) r . For this purpose, we introduce the new r-th Hamiltonian as follows 2 : 2 In formula (8), we make use of the so called exchange theorem (see, e.g., formula (4.6) in [18]) which states that, if χ is a generating function, f (p, q)| (p,q)=exp Lχ(p ′ ,q ′ ) = exp L χ f | (p,q)=(p ′ ,q ′ ) . Therefore, we can apply the Lie series to the Hamiltonian function and, only at the end, rename the variables. For more detailed explanations we defer to the whole section 4.1 of [18]. Here, we do not rename the variables in order to avoid the proliferation of too many symbols.
Since we are asking that the term of polynomial degree equal to r + 2 is depending just on the actions, then we have to determine the generating function χ r in such a way that the following homological equation is satisfied: where Z 0 is the quadratic part, i.e., Z 0 = n j=1 ω j z jzj , and Z r is the new term in normal form we have to define. Hereafter, it is convenient to adopt the standard multi-index notation; for instance, this means that (−iz) ℓ = n j=1 (−i z j ) ℓ j , ∀ z ∈ C n and ℓ ∈ N n . Moreover, for what concerns the complex conjugate variables, the notation has to be read as follows:zl = n j=1zl j j , ∀ z ∈ C n andl ∈ N n . Let us write the generic expansion of the perturbative term f therefore, the generating function that solves the homological equation (9) is such that 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 ω is non-resonant, otherwise ω · (ℓ −l) could vanish. Therefore, we assume that ω satisfies the Diophantine inequality, that is for some fixed values of γ > 0 and τ ≥ n − 1. This condition is not so strict, since for τ > n − 1 it is satisfied by almost all the frequency vectors in R n . The terms with ℓ =l, that cannot be removed by this procedure, depend on the actions only, since z andz appear with the same exponent. As a consequence, we are forced to include them in the new term Z r making part of the normal form, that is 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 change of coordinates induced by exp L χr and to update accordingly all the terms that compose the Hamiltonian. For this purpose, it is convenient to use a notation which mimics a programming code. At first, we put the new terms f ; this simply corresponds to the initial summand with j = 0 in (1). Then, each contribution due to the j-th iteration of the Lie derivative with respect to the generating function χ r (for j ≥ 1) is added to the corresponding class. More precisely, using repeatedly the property described in lemma 2.1, the new summands, which are generated by the application of the Lie series to the terms in normal form (that are denoted with Z s ) and to the perturbative ones (i.e., f (r−1) s ), are gathered in the following way: where the notation a ←֓ b means that a is redefined so as to be equal to the sum given by its previous value plus b. At the end of all these redefinitions, the following expansion of the new Hamiltonian is well defined: Let us emphasize that in this procedure we do not modify the terms that are already in 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 normal form term Z r , while the perturbative terms are different, being redefined as it has been prescribed by formula (14).
The algorithm can be iterated at the next step, by restarting from the Hamiltonian (7) where r − 1 has to be replaced with r.

About the divergence of the Birkhoff normal form
In spite of the fact that we are interested in providing results holding true for a ball of initial conditions B ̺ 0 (0), where ̺ 0 is fixed, now we are going to briefly discuss some asymptotic properties of the Birkhoff normal form.
It is well known that the series introduced by the algorithm constructing the normal form are asymptotically divergent, with the exception of some particular cases. This means 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, because there is an optimal normalization step which minimizes the remainder. The mechanism of divergence is mainly due to the accumulation of the so called "small divisors". In fact, at each step r, new divisors ω · (ℓ −l) (being |ℓ| + |l| = r + 2) are introduced by the definition of the generating function χ r in (11). Because of the redefinitions described in formula (14), the size of the generating functions clearly impacts on the growth of the terms constituting the remainder. If the frequency vector satisfies the Diophantine condition in (12), it is rather easy to see that factorial coefficients O(r!) appear in the estimate of the remainder at the r-th normalization step. This prevents the convergence of the normalization algorithm when it is iterated ad infinitum.
Nevertheless, when the series are estimated on open balls B ̺ (0), with a fixed value of ̺, it is convenient to proceed with the normalization algorithm until the remainder is decreasing. The optimal step r opt is determined, by comparing two subsequent remainders R (r) and R (r−1) . As discussed, e.g., in [18] (see formula (3.30)), this allows to conclude that where C is a positive constant, τ is the exponent appearing in the Diophantine condition (12) and ̺ is the distance from the equilibrium point. As a consequence, an analytical estimate of the remainder at the optimal step can be easily provided; in fact, it is exponentially small with respect to the inverse of a fractional power of the distance from the equilibrium ̺, that in this case assumes the role of small parameter. Near the elliptic equilibrium, the Birkhoff normal form is a good approximation of the real problem and it is convenient to perform a big number of normalization steps. On the other hand, far away from the origin the Birkhoff normal form starts to diverge earlier. Let us mention that in [13] the mechanism of divergence is carefully investigated in a numerical way and it 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 bounds should largely overestimate the effective size of the remainders. In practice, when dealing with explicit expansions, two additional problems have to be tackled: truncations and computational costs. In fact, the Hamiltonian at step r is expanded as in (5) and (15), with Of course, efficient coding and a powerful computing system are more than welcome, because they allow better evaluations of the remainder terms. However, any computer cannot represent all the infinite sequence of terms giving contributions to R (r) . In the next section, we will explicitly provide rigorous upper bounds for the series defining the remainder.

Iterative estimates for Birkhoff normal form
In this section we are going to discuss how to construct a scheme of estimates, in order to keep control of the norms of the Hamiltonian terms, starting from the algorithm described in section 2. In spite of the fact that every purely analytical work based on Birkhoff normal forms exploits a suitable scheme of estimates, our approach significantly differs from all the previously existing ones in the scientific literature for what concerns the following untrivial point, that is clearly highlighted, for instance, in the statements of propositions 3.3 and 3.4. Here, the estimates are designed in such an iterative way, that they take advantage of being applied in the framework of a computer-assisted proof, in order to provide the best possible final results. Before proceeding with the description of the scheme of estimates, the introduction of some new notation is mandatory. For every f ∈ P s whose expansion is the following: we denote by · the functional norm defined as Therefore, in the domain the sup norm of f can be easily estimated as follows: For the sake of simplicity, in the previous definitions we have not introduced more general domains that are polydisks; this has been done with the aim to avoid all the modifications that are necessary to adapt norms and estimates. However, let us recall that in some applications of the Birkhoff normal forms the approach based on polydisks has been useful (see, e.g., [24]).
Hereafter, we will assume to have explicitly built the Birkhoff normal form up to a fixed order r = R I and we will define estimates for the norms of the infinite terms we did not calculate. In order to do that, for every normalization step we will distinguish between two 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 terms with degree greater than R II + 2, that we will control by a sequence of majorants growing in a geometrical way. In the next subsections, we describe how to define iteratively the estimates for the norms of these two kinds of polynomial terms.

Iterative estimates for finite expansions of the power series
Let us suppose that we have already performed the first r − 1 steps of normalization, so that our Hamiltonian is expanded as in (7); moreover, let us assume to know upper bounds for the norm of every term appearing in (7), i.e., We have now to describe how to estimate the norms of the terms that compose the new Hamiltonian H (r) , namely how the estimates F (r−1) s change after a step of Birkhoff normalization. Since the canonical change of coordinates that transforms H (r−1) into H (r) is defined by a Lie series with generating function χ r , we need an estimate for the norm of the Poisson brackets between two functions. In the following lemma, we adapt 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 and suppose that f ≤ F and g ≤ G for some constant G, F ∈ R + ; therefore, In the framework of computer-assisted estimates, sometimes it is more convenient to use the following lemma.
Lemma 3.2 Let χ r and g be polynomial functions with χ r ∈ P r+2 and g ∈ P s+2 , being the expansion of the former function such that Therefore, ∀ j ≥ 1, the following estimate holds true: where Indeed, if we replace D r with (r + 2) χ r , the previous statement easily follows by induction from lemma 3.1. The above definition of D r allows to estimate more carefully the partial derivatives of χ r that appear in the Poisson brackets. This is particularly advantageous 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 is convenient to provide an estimate for the norm χ r ≤ G r and, then, D r can be replaced by (r + 2)G r , for all those r ∈ [R I + 1 , r opt ] for which the coefficients c k,k appearing in the 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 the redefinitions in formula (14) for using inequality (24). Therefore, it is easy to realize that 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 as prescribed 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 those of f (r−1) r divided by ω · (ℓ −l). Let α r be the smallest divisor introduced at the r-th normalization step, i.e., which is well-defined in view of the non-resonance condition assumed in (12). Therefore, we can write χ r ≤ G r , with 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−1) 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) . Proposition 3.3 Let us suppose that the terms of the Hamiltonian introduced at r − 1-th normalization step are bounded as in (21). Therefore, the new Hamiltonian terms of where these new constants are given by the following sequence of redefinitions: being D r = |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), The previous statement is basically a summary of all the discussion explained in the present subsection. In particular, formula (27) can be easily verified, starting from the iterative definitions of the new perturbative Hamiltonian terms in (14) and using the estimates in lemma 3.2.

Power series expansions of Birkhoff normal forms: geometrical increase of majorants for the infinite tails of terms
In this subsection, we are going to describe how to estimate the norm of the terms that cannot be explicitly calculated or iteratively estimated, as it has been described in the previous subsection. In order to do that, it is convenient to bound in a different way the terms constituting the Hamiltonian at the r − 1-th normalization step, i.e., with E, a r−1 ∈ R + . Of course, the estimates above are expected to be less strict than the ones in (21). However, let us imagine to be able to justify the same kind of estimates at the r-th normalization step, holding true for all f (r) s with s ≥ r + 1, with a new constant value a r ; therefore, we could control the sup norm of the remainder (recall equations (5) and (15)) as follows: where we include the upper bounds as iteratively defined by proposition 3.3 when they are available, while for the terms with polynomial degree greater than R II + 2 we use the wanted 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 terms appearing in the expansion (7) of H (r−1) . Therefore, the norms of the terms making part of the new Hamiltonian H (r) satisfy the following inequalities: being where D r = |k|+|k|=r+2 |c k,k | max j {|k j |, |k j |} if the expansion (23) of the generating function χ r is explicitly known, else D r = (r + 2)G r with χ r ≤ G r .
Proof Let us start from the estimate for Z s with 1 ≤ s < r. Since a r−1 < a r , all the estimates in (28) are still valid when a r−1 is replaced by a r . By comparing the polynomial expressions (10) and (13) that are related to the functions f (r−1) r and Z r , respectively, it is easy to realize that the norm of the latter cannot be greater than the former one, because in a normal form term there are just a part of the summands appearing in expansion of the corresponding perturbative function. Therefore, in view of (28) we can write the following chain of inequalities: r ; this ends the justification of the first inequality in formula (30).
Let us now focus on the estimates for f (r) s , with s ≥ r + 1. First, we discuss the case with s = jr + m for m = 1, . . . , r − 1 and j ≥ 1: the new perturbative terms are defined as the sum of f Therefore, inequality (24) allows us to write the following estimate: Using the estimate in hypothesis (30) and the trivial inequality m < r, we obtain Therefore, 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 rewrite formula (14) in the following more extended form: where we have separated the term generated by f (r−1) r . Since χ r solves the homological equation (9), we can rewrite the first two contributions as follows: Therefore, for j ≥ 2 we can write where we used again inequality (24) and F (r−1) r as upper bound for the norm of Z r . The second inequality in formula (30) follows from calculations similar to those we have previously described for the case where s is not a multiple of r.
Let us emphasize that the same scheme of estimates is valid also for a resonant Birkhoff normal 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 the definition of the smallest divisor α r and, as a consequence, the estimate for the generating function χ r .

Computer-assisted proofs about effective stability
Here we are going to use the scheme of estimates described in the previous section for producing results on the effective stability in the neighborhood of an elliptic equilibrium point. We will study separately two different cases: a non-resonant Birkhoff normal form and a resonant one. For what concerns the effective stability in the vicinity of equilibrium points, up to the best of our knowledge other computer-assisted rigorous estimates are available in [8] only, but such a reference has to be considered rather outdated, because of some technical reasons affecting the final quality of the results (see also the comments in the Introduction). Moreover, an application of this kind of approach to a resonant Birkhoff normal form is entirely new.

Estimates on the escape time for non-resonant Birkhoff normal forms
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) depending on the actions only and R (r) = O( I (r+3)/2 ). Therefore, we can calculate the time derivative of the action I j = z jzj as follows: because {I j , Z (r) (I)} = 0. Using the expansion (15), the estimate on the sup norm in (20) and lemma 3.1, we can find the following upper bound for the sup norm of the time derivative of the actions 3 : We can now provide rigorous estimates for the norms of all the terms appearing in the previous series. More precisely, we can use the estimates F (r) s (as recursively defined in (27)) for the first R II terms. Moreover, for what concerns the first R I terms, we can do even better, because we have explicitly performed R I steps of the normalization algorithm; therefore, for such terms the estimates can be replaced by the actual values of the corresponding norms. For all the remaining terms, we limit ourselves to use the uniform estimate in (30). Thus, we just need to compute the value of the upper bound in the r.h.s. of the following inequality: Since the remainder is expected to be exponentially small at the optimal step, it is clear that 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 equilibrium point, where the optimal normalization step is far beyond the number of steps we are able to explicitly perform, by using any kind of software that is specialized for doing computer algebra.
Let us now suppose that the initial value of the action vector belongs to an open domain, e.g., I(0) ∈ ∆ ̺ 0 ; we aim to determine a time T > 0 (as long as possible) such that we can ensure that When an estimate for |İ j | ̺ ∀ j is available, it is convenient to define T as to obtain the wanted confinement. Combining (33) and (34) with (36), we obtain the following lower bound about the escape time from ∆ ̺ : 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 ̺ 0 making a further optimization 4 of the expression above is ̺ 0 = [(r opt + 1)/(r opt + 3)] 1/2 ̺. 4 When r = r opt , for the sake of simplicity, let us assume that the geometrical decrease of the terms in the series appearing at the denominator of equation (37) is so sharp that it can be approximated by its first term, then we can write T (̺, ̺ 0 , r opt ) ≃ g(̺), being g( and C a suitable positive constant. After having remarked that g(̺ 0 ) = lim ̺→∞ g(̺) = 0, one immediately realizes that the function g : [̺ 0 , ∞) → R takes its maximum value in correspondence to the solution of the equation g ′ (̺) = 0, i.e., ̺ = [(r opt + 3)/(r opt + 1)] 1/2 ̺ 0 . Therefore, in applications where the value of ̺ is considered as fixed, we will compute the lower bound for the escape time T ̺ , [(r opt + 1)/(r opt + 3)] 1/2 ̺ , r opt according to formula (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 denominator appearing in formula (37) would be made just by its third summand with the series starting 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 obtain the asymptotic law for the purely analytical estimate about the escape time from ∆ ̺ , i.e., where̺ is a positive constant. For more details we refer to theorem 3.5 in [18], where a complete proof of such a lower bound can be found.

Resonant Birkhoff normal forms: bounds on the escape time
In the case of a resonant Birkhoff normal form, the estimate of the diffusion of the actions is not so straightforward as for the non-resonant one. Indeed, we cannot use directly the remainder 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 Birkhoff normal form up to order r, in such a way that the Hamiltonian of the system has the following structure: where I ∈ R n and ϕ ∈ T n are action-angle canonical coordinates, being ϕ n the resonant angle. Therefore, the estimates in (33) and (35) still hold true for what concerns the diffusion of the actions I 1 , . . . , I n−1 , since we removed the corresponding angles from the normal form. On the other side, it is more difficult to provide good estimates for the diffusion of the action I n , that is conjugated to the resonant angle ϕ n . Since the computation of the stability time T as in formula (37) is valid when all the actions are confined in the domain ∆ ̺ , it is necessary to control the variation of the action I n in order to use such an estimate about the value of T . For this purpose, it is convenient to combine the conservation of the total energy of the system, i.e., the Hamiltonian, with the slow diffusion of the actions I 1 , . . . , I n−1 . More precisely, by omitting the dependency on the index r in formula (39), we can write the following equation: where E is the energy level and we denote withZ 1 the sum of all the normal form terms but the linear ones with respect to the actions (that are still of the form ω · I as in the initial Hamiltonian (4)). By assuming the confinement of all the actions in a ball of radius ̺ 2 , we can evaluate an upper bound of and a lower bound of This allows us to give an estimate of ∆I n = I n; max − I n; min .
Of course, the excursions experienced by the values of the normal form Z mainly depend on the linear terms; this explains why we have written them separately in the equations above. The bounds for the sup-norm and inf-norm appearing in formulae (41)- (42), respectively, can be easily calculated, by using the inequalities (21) and (28). Since the normal form is composed by a finite number of terms, the estimate for ∆Z 1 can be evaluated very strictly when its expansion is explicitly known.
Recalling that the estimate for the variation of I n holds true only if all the actions are confined in a ball of radius ̺ 2 , the maximal radius on I n in order to validate such an estimate is defined as 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)| < ̺ 2 0 ∀ j = 1, . . . , n − 1 and |I n (0)| < (̺ * n ) 2 . Let us remark that the more ̺ 0 is close to ̺, the more the variation of I n is prescribed to be small; at the same time, as it follows immediately from formula (37), the stability time decreases for ̺ 0 going to ̺. In order to balance these two effects, it is convenient to define the parameter ̺ 0 in a different way to what has been done at the end of the previous subsection. In the next section, we will discuss in a practical example how the radius ̺ 0 of the open ball containing the initial conditions can be determined in the case of resonant Birkhoff normal forms.
5 Effective stability of a couple of physical models

An application to the Hénon-Heiles model
The Hénon-Heiles model was introduced in [27] for studying the dynamics of a star in a galaxy. The Hamiltonian that describes the system is composed of a couple of harmonic oscillators perturbed with cubic terms, i.e., Here, we consider the non-resonant case only; in particular, we fix ω 1 = 1 and ω 2 = −( √ 5 − 1)/2, so that the angular velocity vector is Diophantine with τ = 1 according to its definition in (12), because the ratio of its components |ω 1 /ω 2 | is equal to the golden mean 5 . We focus on the study of the dynamics in a neighborhood of the origin. This Hamiltonian belongs to the general class described in (2), for which we have explained how to perform the normalization procedureà la Birkhoff.
Using X̺óνoζ, which is an algebraic manipulator specially designed to implement approaches that are common in the framework of Hamiltonian perturbation theory (see [23]), we explicitly constructed the Birkhoff normal form for this model by representing terms having total polynomial degree less than or equal to 102 and by performing R I = 100 normalization steps. Using the iterative estimates described in section 3, combined with the norms that have been determined for the first R I = 100 terms (for which we have explicit expansions), we provided rigorous estimates for the remainder and the time derivative of the actions, for any fixed value of ̺. Indeed, for the first R II = 1500 terms we iterated the estimates of the norms as described in proposition 3.3, while for the rest of the infinite terms we used the upper bounds reported in formula (28), updating the value of a r as prescribed in (31); finally, we produced an estimate for both the remainder and the escape time T , by using formulae (29) and (37). For each value of ̺ we considered, the computational algorithm was stopped, after having identified the optimal step r opt as the one corresponding to the minimum value of the remainder; thus, we fixed our final estimate for the escape time in such a way that T = T (̺, [(r opt + 1)/(r opt + 3)] 1/2 ̺, 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 parameters R I and R II , allowed to obtain a fully rigorous computer-assisted proof of existence for the KAM torus related to the golden mean ratio of frequencies in the case of the forced pendulum 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. In order to fit with a so successful approach, in our programming codes we have implemented validated numerics exactly in the same way. Therefore, every coefficient appearing in all the polynomial expansions that are explicitly computed (for Hamiltonians or generating functions) are replaced with intervals and all the mathematical operations between intervals are performed by taking into account the round-off errors. For what concerns the iteration of the estimates, validated numerics is used in order to properly provide all the needed bounds. In summary, the whole computational procedure is aiming at ensuring that the final result is fully rigorous. Appendix A is devoted to a short introduction to validated numerics, that is described with the only goal of explaining our implementation of it. Furthermore, in appendix B we have included a pedagogical discussion of a computer-assisted estimate of the stability time for the Hénon-Heiles model in a very simple situation: we have focused on a so small neighborhood of the equilibrium point (i.e., we have fixed ̺ = 0.0001) that a few normalization steps are enough to obtain a final result that is poor but still meaningful. ̺ 0 ̺ r opt a r log 10 |R (ropt) | ̺ log 10 |İ j | ̺ log 10 T 9.96e-04 1.00e-03 232 1. In Table 1 we have included our results for some values of the distance ̺ from the equilibrium point: let us emphasize that each row corresponds to a single computer-assisted proof. Therefore, we think that every item of that list can be useful for comparisons with the results eventually obtained by using other techniques. Moreover, in Figure 1 we have reported the main content of Table 1, that is about the behaviour of both the optimal step and the estimate for the escape time as functions of ̺. In the left box of Figure 1 one can appreciate that the limit value of ̺ for which we are able to produce stability results is around 0.1 . This limitation is due to the rate growth of a r , whose values get close to 10 just after a few normalization steps; this prevents the convergence of the normal form, when ̺ 0.1, because the common ratio of the majorants series is too large: a r ̺ ≥ 1. Looking at the definition in (31), one immediately realizes that having sharp estimates for the generating functions is crucial to obtain good final results. This fact is emphasized also by the occurrence of the plateau in the plot of the function r opt (̺), that is in correspondence with the number of normalization steps R I = 100 for which we explicitly compute the expansions. Indeed, the worsening effect induced by the transition to the mere iteration of the estimate is so remarkable that we have to strongly decrease the ball radius ̺ in order to take a real advantage of performing more normalization steps. With the only noticeable exception of the already mentioned plateau, the plot of the optimal normalization step looks consistent with the expected law, i.e., r opt (̺) ∼ 1/ √ ̺ (see formula (16)). The agreement between the observed behavior and the expectations is even better for what concerns the plot in the right box of Figure 1. In fact, from the asympotic law (38) one can deduce that log T ∼ 1/ √ ̺, that is coherent with the fact that such a semi-logarithmic plot fits rather well with a straight line. Thus, we can conclude that our computer-assisted estimates about (the lower bound of) the stability time preserve its main property: being exponentially large with respect to the inverse of the square root of the distance from the equilibrium point. Let us conclude this subsection with a short discussion about the choice of the parameters R I and R II . Obviously, it is convenient to perform the largest possible number of explicit steps. Therefore, our choice of the parameter R I is mainly due to the computational cost in time and memory that is needed by our algebraic manipulations. The choice 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 to the 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 such a situation, choosing larger values of R II has the only effect to increase the computational time, without any substantial improvement for what concerns the estimate of the remainder. Conversely, for larger values of ̺, the common ratio a r ̺ can approach 1 and the contribution of the infinite tails becomes predominant. Therefore, it is convenient to increase R II in order to estimate more properly the remainder. This explains why, in our calculations, we have usually fixed R II = 1500, but we have suitably increased such a value when a r ̺ 1. However, the computational cost of the iteration of the estimates is negligible; to fix the ideas, when R II = 1500 it took less than a minute, while performing explicitly R I = 100 normalization steps required about one day of CPU-time on a computer equipped with an Intel quad-core I5-6600 -3.3 Ghz and 16 GB of RAM. Let us recall that such an impressive difference of computational cost is due to the fact that the representation of the polynomials defined by the normalization algorithm requires a huge occupation of the memory if compared with a single upper bound on the corresponding norm and, therefore, many more operations are needed to manipulate their explicit Taylor expansions. For planning new applications it is also important to know the scaling law of the computational time T CPU as a function of R I : we have found that for the first 100 normalization steps such a law is not extremely sharp but it can be bounded from top (and conveniently approximated) in such a way that T CPU ∼ R 8 I .

Applications to the asteroidal motions of the Trojans
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 Lagrangian equilibria for the CPRTBP. As a first preliminary step, a rotating frame Oxy with origin O in the centre of mass of the two primary bodies and rescaled physical units is introduced. The Hamiltonian of such a system is usually written as follows (see, e.g., [47]): where µ is the mass ratio between the primaries and (p x , p y ) ∈ R 2 are the kinetic momenta that are canonically conjugate to the positions (x, y) ∈ R 2 , respectively. Therefore, it is convenient 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 where minus and plus refer to the equilibrium points L 4 and L 5 , respectively. In these new variables, the Hamiltonian takes the form After having skipped a constant term, the basic Taylor expansion of the previous Hamiltonian can be written as H(P X , P Y , X, Y ) = Q(P X , P Y , X, Y )+O (P X , P Y , X, Y ) 3 , where the 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 (9 − √ 69)/18, then the equation of motion are linearly stable. This means that there is a class of linear canonical transformations 6 (P X , P Y , X, Y ) = C(y 1 , y 2 , x 1 , x 2 ) conjugating the quadratic approximation to a couple of harmonic oscillators, i.e., Q C(y 1 , y 2 , x 1 , x 2 ) = ν 1 (x 2 1 + y 2 1 )/2 + ν 2 (x 2 2 + y 2 2 )/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 the normalization procedureà la Birkhoff, i.e., being the angular velocities ν 1 and ν 2 such that As an immediate consequence, we have that ν 1 → 1 and ν 2 = O(µ) for µ → 0. The second frequency 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 constructing the Birkhoff normal form (that has been widely discussed in section 2), very small divisors can be introduced because of the fact that |ν 2 | ≪ 1. Therefore, it is natural to expect that a 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 triangular equilibrium point are commonly called "Trojans", according to the tradition started by Max Wolf at the beginning of the XX century. In fact, he decided to choose names from Homer's Iliad for the first bodies which were observed in the vicinity of L 4 or L 5 , in the system having Sun and Jupiter as primary bodies.
In the next subsection we will compare the performances of these two kinds of Birkhoff normal forms, by considering a few realistic values for the parameter µ, which correspond to 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 that is comparable with an overestimate of the residual life of the Sun in the main sequence, i.e., 6 × 10 9 years.

Effective stability of trojan celestial bodies in the Solar system
As it has been discussed in the introduction, an approach merely based on the Birkhoff normal form is not enough to obtain realistic results about the stability of the Jupiter Trojans in the framework of the CPRTBP model. In [24], it has been shown that orbital motions starting from initial conditions contained in the domain (I 1 , I 2 ) ∈ [0, 0.0008] × [0, 0.0005] are effectively stable. This is far from explaining why the regions in the proximity of L 4 and L 5 are so populated; indeed, that domain covers the observational data of a very ̺ 2 0 ̺ 2 T 2.49e-04 2.59e-04 6.36e+08 2.47e-04 2.57e-04 1.01e+09 05e-04 1.83e-04 2.07e-04 5.93e+08 2.02e-04 1.80e-04 2.04e-04 7.23e+08 Table 2: Comparison for the estimates on the stability time between the non-resonant and resonant Birkhoff normal forms. The Jupiter case (µ ≃ 0.000954) with T e.l.t. ≃ 5 × 10 8 .
small fraction of Trojans. The results get even worse when all our rigorous estimates about both the remainder and the stability time are taken into account. As it is clearly shown in the left hand side of Table 2, the construction of the non-resonant Birkhoff normal form allows us to ensure the effective stability for initial conditions such that (I 1 , I 2 ) ∈ [0, 0.00025] × [0, 0.00025]. In that table we reported only a couple of results for two different (and very close) values of ̺. In fact, for this particular problem, it has no 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 stability time is comparable with the so called expected life-time T e.l.t. of the system. This number is rescaled with respect to the revolution period of the celestial body, i.e., we impose that T e.l.t. = (6 Gyrs) · ν 1 2π ≃ 5 × 10 8 (in number of Jupiter revolutions, because ν 1 is expressed in rad/yrs). As it has been highlighted in Figure 2, a small decrease of the value of the radii containing the initial condition ̺ 0 can very significantly change the stability time (because of the exponential dependency of T on the inverse of ̺). Such a behaviour is evident also in Table 2, where with a change of the 1% of ̺, the stability time is almost doubled in the non-resonant case. Let us emphasize that in our estimates, we use the same radius for both the actions, while looking at the initial conditions of the real trojan asteroids for Jupiter in [17], the value of the second action is usually smaller with respect to the first one. Therefore, we believe that considering the initial conditions of the actions in 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 the results 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 to explain how to determine the radius ̺ 0 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 ̺ 0 in order to optimize our evaluation of the lower bound of the escaping time T . Instead, here we focus our attention on the consistency of the procedure, in the sense that we aim to ensure that the motion is confined in the domain ∆ ̺ for all |t| ≤ T . Therefore, we fix ̺ 0 ∈ (0 , ̺) as the unique solution of the equation where ν 1 and ν 2 have opposite signs in view of (47) and β < 1 is another parameter we have the freedom to choose, because it essentially represents the price we are willing to pay in terms of a further restriction concerning the domain of the initial value of the resonant action I 2 (0) . This last concept definitely deserves a more detailed explanation. Indeed, in agreement with (44), we define (̺ * 2 ) 2 = ̺ 2 −∆I 2 , with ∆I 2 given as in (43). Since the main contribution to the maximal variation of the resonant action is due to the linear terms (see formulae (41)-(43)), we have that ∆I 2 ≃ −ν 1 (̺ 2 − ̺ 2 0 )/ν 2 and, therefore, (̺ * 2 ) 2 ≃ β̺ 2 . Let us rephrase a conclusion discussed at the end of section 4.2, by adapting it to the present context: we have that This ends the explanation on the meaning of β as the parameter ruling the restriction on the set of values of I 2 (0). For what concerns the choice of the value of β, we are interested in setting it very close to 1, in order to enlarge the domain of stability in I 2 (0) as much as possible; on the other side, we have to take into account oscillations of the value of the resonant action that are due to the angular dependence of the normal form partZ 1 (I, ϕ 2 ) (see formulae (41)-(42)). Indeed, this so called modulating term for the resonant action I 2 is such thatZ 1 = O( I 2 ). This remark allows to imagine a further optimization procedure on the determination of β, that in principle would not be so difficult, due to the already mentioned fact that the maximal variation ∆I 2 of the resonant action is linear with respect ot the the actions (I 1 , I 2 ) in view of formulae (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 simply set β = 0.9 . Indeed, we consider that an eventual further enlargement of the domain of stability in I 2 (0) is not so crucial, because it cannot be greater than 10 %; therefore, it cannot substantially improve our results, as it is clearly shown in the summary of the comparisons reported in the final Table 6.
The right hand side of Table 2 includes the results based on the construction of a resonant Birkhoff normal form, whose performance is worse with respect to the nonresonant 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 norms of the generating functions χ r for the resonant and non-resonant constructions, in the case of the system having Sun and Jupiter as primary bodies. The increase of the norms proceeds more and less at the same rate up to the first 25 steps; afterwards, the norms of the generating functions related to the resonant normal form start to grow up faster than  (9); on the other hand, the normal form part is larger in the case of the resonant construction and this clearly has an impact on the size of the generating function. In fact, when the r-th normalization step is performed, new perturbative terms are introduced by the application of Lie series to the Hamiltonian H (r−1) ; the main ones are 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 new terms and the main new contribution is expected to be due to L χr Z 2 . This explains why in the non-resonant case the growth of the generating functions has two different trends for odd and even indices, as it is clearly shown by the plots reported in Figure 3. On the contrary, in the resonant case the behaviour of the norms of the generating functions is more regular, at least up to a threshold value of r, for which the worsening effect due to the terms generated by the normal form part becomes predominant with respect to the advantage gained because of the better accumulation of small divisors.
Let us now focus on the comparison between the performances of the resonant and non-resonant constructions, when they are made for different values of µ: the results are summarized in the top-right, bottom-left and bottom-right boxes of Figure 3 and in Tables 3-5, which refer to the systems having Sun-Mars, Sun-Uranus and Saturn-Janus as 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 is remarkably smaller with respect to the corresponding non-resonant ones. Therefore, the series introduced by the resonant normalization procedure are convergent in a bigger neighbourhood of the origin; in other words, in such a situation we are able to iterate the computer-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 µ. This allows us to emphasize that the comparison between the computer-assisted estimates based on the non-resonant Birkhoff normal form and the resonant one are more and more in favour of the latter, when the value of µ tends to zero.

Conclusions and perspectives
Since the very beginning of our research project, one of our main goals was to complement the Birkhoff normal form with a coherent scheme of computer-assisted estimates, in order to 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 such a 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 such an approach).
We have successfully applied our procedure to two simple systems, which are close to elliptic equilibrium points such that the quadratic approximation of the Hamiltonian is not convex (so preventing a trivial proof of stability using the energy as Lyapunov function): the Hénon-Heiles model and the Circular Planar Restricted Three-Body Problem ̺ 2 0 ̺ 2 T 6.00e-07 6.37e-07 3.10e+12 5.89e-07 6.24e-07 5.40e+12 18e-05 1.10e-05 1.18e-05 3.50e+12 1.15e-05 1.08e-05 1.15e-05 6.83e+12 Table 5: As in Table 2 for the Janus case (µ ≃ 3.36 × 10 −9 ) with T e.l.t. ≃ 3 × 10 12 . 6.00 × 10 −7 1.10 × 10 −5 18.33 Table 6: Comparisons between the values of the radii ̺ 2 0 and (̺ * 2 ) 2 which refer to the stability domains for the non-resonant Birkhoff normal form and the resonant one, respectively. 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 the first column.
(CPRTBP). For the latter, we have easily adapted our approach to the construction of the Birkhoff normal form of resonant type and we have shown that the consequent computerassisted estimates are in a better position when the mass ratio between the primary bodies is very small. Nevertheless, the extent of the domain covered by our results is still far from 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, for what concerns the applicability to realistic physical models, that are of interest, in particular, in Celestial Mechanics and Astronomy. We think that in these fields there still are problems that are open and fundamental 7 , where rigorous proofs of stability can be tackled by using an approach purely based on the normalization algorithm for KAM tori, which can ensure a perpetual topological confinement in models with two degrees of freedom. Nevertheless, there is a much wider range of possible applications to Hamiltonian systems defined on phase spaces of higher dimensions. In such a framework, a very promising strategy is based on the local construction of the Birkhoff normal form in the neighborhood of an invariant KAM torus: it has already started to provide some remarkable results (see, e.g., [22]). Here, in the context of the discussion about the perspectives of our computational method, it is natural to describe the scaling properties of our approach with respect to the number n, i.e., the degrees of freedom of the Hamiltonian model. It is well known that the remainder of the Birkhoff normal form is exponentially small with respect to the distance ̺ form the equilibrium point; this allows to deduce the following best possible estimate about the stability time: T ∼ exp[(̺/̺) 1/n ], that is corresponding to the most non-resonant frequencies in the small oscillations limit, being ̺ a positive constant (recall formulae (12) and (38)). Such an asymptotic law is intrinsic to an approach based on Birkhoff normal forms and we have shown that it is in agreement with the results produced by our computational method. This is due to the fact that for small distances ̺ the so called optimal normalization step r opt can largely exceed the maximal polynomial degree R I + 2 of the Hamiltonian terms whose expansions are explicitly 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 degrees of freedom n. Nonetheless, the accuracy of the final estimates about the stability time remarkably 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 for the parameter R I ruling the size of the expansions stored on a computer. Since a polynomial function of maximal degree R I + 2 and depending on 2n variables hosts a number of terms that is O R 2n I , it is natural to expect the following scaling law for the threshold corresponding to the deterioration of our computational results: R I;max ∼ M 1/(2n) , being M a suitable positive constant.
In the very particular case of the PCRTBP, there are other possible sources of further improvements. Indeed, the normal form construction that has been designed in [42] adopted (since the very beginning) canonical coordinates which are much more suitable for such a kind of Celestial Mechanics model. This is the reason why the approximation provided by that normal form describes much more carefully the dynamics, when it is compared with the numerical results based on the Birkhoff normal form. We think that our computer-assisted estimates can be adapted to complement also the constructive algorithm explained in [42]. This should provide results about the effective stability holding true in wider domains, hopefully covering an important fraction of the trojan asteroids.

Acknowledgments
This work was partially supported by the "Mission Sustainability" programme of the Università degli Studi di Roma "Tor Vergata" through the project IDEAS (E81I18000060005) and by the "Progetto Giovani 2019" programme of the National Group of Mathematical Physics (GNFM-INdAM) through the project "Low-dimensional Invariant Tori in FPUlike Lattices via Normal Forms". The authors acknowledge the MIUR Excellence Department Project awarded to the Department of Mathematics of the University of Rome "Tor Vergata" (CUP E83C18000100006), in particular, because of the availability of the computational resources.

A Basics of validated numerics on a computer
In order to give a completely rigorous support to our computer-assisted proofs, in the present work we have implemented everywhere validated numerics in two complementary ways. First, interval arithmetic has been used to calculate the coefficients of the Taylor expansions 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., norms of functions [stability times, respectively]. The present appendix is devoted to summarize our approach to validated numerics, that essentially follows what is described in section 3 of [31].
All our codes are written in C language and all the quantities of interest are computed using the double type in C. The set of floating point numbers of double type that are representable on a computer is defined as follows: where σ is the exponent and the digits d 0 .d 1 d 2 . . . d 52 make up the so called mantissa 8 m; they appear in the binary scientific notation x = ±m 2 σ of each number x belonging to the set R. Therefore, an overflow [underflow] situation can occur when a computational operation attempts to generate a number whose absolute value is greater [smaller] than M = (1 − 2 −53 ) · 2 1024 [than m = 2 −1074 and it is different from 0, respectively]. Let us generically use the symbol * to refer to any elementary arithmetic operations +, −, · and /; moreover, the result provided by a computer when it performs that same operation will be denoted with ⊛. According to such a setting, for instance, a ⊕ b is the computer output for the usual sum a + b, where both a and b belong to R. Therefore, the machine epsilon for double type floating point numbers in C (hereafter, simply ε) is defined as By comparing the definition above with that in (49), one can easily realize that ε = 2 −52 and ε/2 = min η ∈ R : 1 ⊖ η < 1 .
The 64 bit IEEE standard ensures that, if a * b / ∈ R and m < |a * b| < M, then either In words, these prescriptions can be summarized as follows: when a generic 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 9 .
In order to prevent the occurrence of overflows or underflows generated by an elementary arithmetic operation ⊛, we restrict to the "safe range" (according to definition 3.2 in [31]), which is the following set of numbers: (50) 8 In the first row of formula (49), d 52 is fixed to be zero, while d 0 is set to 1 for the numbers appearing in the second row, according to the so called "hidden bit" rule. By taking into account that the integer exponents between −1023 and 1023 require 11 bits to be represented and 1 bit is needed for the sign, it should be evident that the information to be stored for each double type number must be spread over 8 bytes (=64 bits). 9 Also the square root √ enjoys such a peculiar property (see, e.g., section 3.1 of [45]).
We can rephrase the prescriptions of the 64 bit IEEE standard for the arithmetic elementary operations in a very quantitative way, by referring to the relative error, i.e., 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 real numbers 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 square roots 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 by the factor 1 + ε [by 1 − ε] is enough to increase [decrease, resp.] their absolute value by the last significant digit. This last remark joined with the prescriptions of the 64 bit IEEE standard allows us to establish a few simple rules, in order to rigorously perform interval arithmetic. Let us introduce the following binary operators ⊛ + : S × S → R and ⊛ − : S × S → R, that are defined as follows: where, again, * generically denote +, −, · and / (with the only exception of the division by zero). We emphasize that the inequalities in formula (52) are fundamental, because they allow us to implement the rigorous computations of upper bounds (or lower ones) for quantities that have to be estimated. Moreover, the validated numerics can be applied to interval arithmetic by suitably extending the binary operators introduced in (51). In fact, let ⊛ ± : Because of the definition of ⊛ ± and in view of formulae (52)-(53), we can write the following relation that is fundamental for the (rigorous) interval arithmetic: where, once again, the case with 0 ∈ [b − , b + ] and * representing the division must be excluded. Let us stress that two (possibly redundant) list of elements appear after the minimum and the maximum in the definitions (53) of c − and c + , respectively, in order to take into account of the effects induced by the signs, when we are dealing with products and divisions. To fix the ideas, it is convenient to realize that, for instance, the rigorous extension of the sum to intervals, can be defined in a much shorter way, i.e., [a − , Since the algorithm constructing the Birkhoff normal form is based on Lie series (therefore, Poisson brackets), one can immediately realize that the computations of the coefficients appearing in the expansions involve elementary arithmetic operations, only. Thus, their rigorous extension to intervals is enough to perform validated numerics on the truncated expansions of the Hamiltonian as far as they can be explicitly computed. Let us stress that the whole computational algorithm can be iterated provided that each generic arithmetic operation of type [a − , generates a new result such that c − , c + ∈ S. This is the reason why, in our codes, we included tests, in order to verify that such a condition is always satisfied; if it is not so, the running of a program is immediately stopped. The situation concerning the rigorous bounds on the estimates is slightly more complicated, because, very quickly, it occurs a violation of the condition that all the majorants of the norms are given by numbers belonging to the "safe range" set S. This is because the estimates can be iterated for many more normalization steps with respect to the explicit calculation of the expansions, each of them requiring a huge occupation of the memory if compared with a single upper bound on the corresponding norm. Due to the dramatically fast growth of the norms of terms composing the series (that are asymptotically diverging for the number of normalization steps going to infinity), the limit max S is usually trespassed even for values of R II that are relatively small with respect to those considered in the applications described in section 5. Let us recall that we are interested in iterating the estimates of the norms for a large number of normalization steps R II , in order to obtain better results. Therefore, it is convenient to represent the logarithm of the (positive values of the) upper bounds of the norms. For such a purpose, we have to introduce a function log + : S ∩ R + → S such that log x log + x for all positive x ∈ S. This can be done in a rather obvious way, by adapting the approach described in section 3.1 of [45], in order to combine validated numerics with a truncated series expansion of a logarithm. To fix the ideas, let us limit to the case x ∈ S∩[1, ∞), then we define 10 log + x = 2 n ⊙ + log + ξ n , being (ξ n ⊖ + 1) ∈ S ∩ [0 , 0.01] such that ξ 0 = x and ξ j = ξ j−1 + ∀ j = 1, . . . , n, where n is the minimum nonnegative integer such that ξ n ≤ 1.01 and the function √ · + is nothing but (1+ε) √ ·, in agreement with what has been explained in the corresponding footnote 9 . The rigorous upper bound for log 1 + (ξ n ⊖ + 1)) ∀ ξ n ∈ S ∩ [1 , 1.01] is introduced in a similar 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 above with x ≥ 1.
The translation of the iterative estimates in terms of upper bounds on the logarithms of the norms is obvious, when just products and divisions are involved. For instance, let us focus on the last two definitions appearing in the statement of Proposition 3.3, i.e., D r = (r + 2)G r , with G r = F (r−1) r /α r . In our codes we can write the corresponding rigorous estimate as log D r = log + (r + 2) ⊕ + log F (r−1) r ⊕ + log + 1 ⊘ + α r , where we 10 Let us stress that in order to ensure that log x ≤ log + x where log + x = 2 n ⊙ + log + ξ n and provided that log ξ n ≤ log + ξ n ∀ ξ n ∈ S ∩ [1 , 1.01], we exploit the fact that integer powers of two belong to S. mean that log F (r−1) 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 a further example, let us focus on (31): in order to properly define an upper bound for the logarithm of a r = a r r−1 + (r + 1)D r 1/r we have to compute log + (x + y), where the values of log x = r ⊙ + log a r−1 and log y = log + (r + 1) ⊕ + log D r have to be considered as known, because both log a r−1 and log D r have been preliminarily estimated by some rigorous computations; 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 exceeding max S. Without any loss of generality, let us assume that x ≥ y, therefore, it is convenient to set log + (x + y) = log x ⊕ + log + 1 ⊕ + exp + log y ⊖ + log x .
The whole of the elementary operations 11 described in the present appendix represent the bare minimum that is sufficient to implement validated numerics, in order to make fully rigorous our computer-assisted proofs.

B A complete example of application to the Hénon-Heiles model in a case with short expansions
The aim of this appendix is very pedagogical: we explain step-by-step an application of the algorithm described in the sections 2-4. For the sake of clarity, here we focus again on the Hénon-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 described in the present appendix deals with the case where R I = 2 and R II = 5; this means that the expansions of larger polynomial degree to be explicitly computed are quartic while the iteration of the estimates provides upper bounds for terms up to the seventh degree.
In our opinion this choice is a good balance: on the one hand it is non-trivial, on the other hand the representation of the Hamiltonians is still readable because it is not too large. Moreover, we have decided to further simplify our example, by omitting the study of our better estimate for the stability time as function of the radius ̺, whose value is fixed to ̺ = 0.0001. First, let us emphasize that, during our computer-assisted proofs, a Hamiltonian H (r−1) (the Taylor expansion of which is explicitly given in equation (7)) is represented by a set 12 gathering both polynomial functions and numbers: with s > R II . Let us recall that the values appearing in the second and third rows of definition (56) are expressed in logarithmic scale because some of them can eventually exceed the "safe range" set S of representable numbers 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 − 1 has to be replaced with r). As a short reformulation, this means that the 12 Let us stress that in formula (56) some unpleasant lower indexes (that are involving either the minimum or the maximum between nonnegative integers) are somehow unavoidable, because the number of already performed normalization steps, i.e. r − 1, can be greater than R I or not; therefore, either the explicit expansions of f (r−1) min{r,RI} , . . . , f (r−1) RI are missing or the majorants log Z RI+1 , . . . , log Z min{r−1 , RII} disappear in the list of elements making part of the set S (r−1) . 13 It is not easy to code with X̺óνoζ, which is the algebraic manipulator we actually used for all the 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 in the framework provided by Mathematica can be more helpful for a reader that is interested in reproducing our results. We have written such a code (named birkhoff HH.mth) in a version computing exactly the coefficients as algebraic numbers. It makes part of the package freely available at the web address http://www.mat.uniroma2.it/∼locatell/CAPs/BirkCAP AppsA-B.zip and its output is easily converted in rigorous upper bounds on the norms by another program (named estimates.c). Also this further code is included in that same software package and is in charge to compute rigorously all the upper and lower bounds that are needed, in order to conclude the computer-assisted proof of the same result discussed in the present appendix.
pair of positive values for E and log a 0 would be acceptable. However, since the value of the common ratio a r−1 (that is related to the majorants of the infinite tail of terms appearing in the remainder of H (r−1) ) affects the next one, i.e., a r , through the relation (31), we have found reasonable to set E and a 0 in such a way to describe the relative increase of f (0) 1 with respect to Z 0 . 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 0 = log F ) and the definitions of the elementary arithmetic operators allowing to implement validated numerics (like, e.g., ⊘ + ) are given in the previous appendix A.
First normalization step: r = 1. We have to put in Birkhoff normal form the cubic term. In order to do that, we solve the homological equation L χ 1 Z 0 + f (0) 1 = Z 1 and, in view of equation (11), we obtain The new normal form term Z 1 is equal to 0, due to the fact that the normalization step is odd (see (13)). We have now to apply the Lie series of generating function χ 1 to H (0) ; since we are considering terms up to degree 4, the only contributions we have to calculate are 1 2 L 2 χ 1 Z 0 and L χ 1 f 1 ; in fact, according to formula (14), they will produce terms of degree 4, which are collected together in the definition of f 1 . In order to properly describe all the Hamiltonian terms at the end of the first normalization step, we have to update the parameters appearing in the second and in the third row of the generic representation described in (56). Since we know the expansion of the generating function χ 1 , we can use the formulae in Proposition 3.3 to rigorously compute the upper bounds D 1 and log(F (1) s ) for 3 ≤ s ≤ 5; moreover, we can compute the norm of 15 f (1) 2 , that is log + ( f (1) 2 ) = 2.446291. We define log a 1 = 2.599403 by using formula (31). Therefore, the representation of the new Hamiltonian H (1) can be summarized by the following set: 2 , 6.685803, 8.979949, 11.16873, 0, 2.599403 .
Finally, we use formula (29) to determine the rigorous estimate of the sup norm of the remainder after having completed the first normalization step: log(|R (1) | ̺ ) = −34.71383.
Second normalization step: r = 2. For what concerns the second step, we proceed in an analogous way with respect to the first one. The homological equation we have to solve is L χ 2 Z 0 + f (1) 2 = Z 2 . We compute the generating function χ 2 by using formula (11) and the new normal form term Z 2 as defined in (13). All the other terms generated by 15 For brevity reasons, here we have decided to not include the expansion of f 2 . However, its normal form part, i.e., Z 2 , is written within the following description of the second normalization step.
As a new estimate for the remainder, we obtain log(|R (2) | ̺ ) = −39.36487. Since it is smaller with respect to the upper bound for the same quantity at the first normalization step, then it is convenient to further iterate the algorithm.
Third normalization step: r = 3. As a major difference with respect to the second step, now an explicit expansion for f is not available and, therefore, the same holds for the generating function χ 3 . Instead of computing D 3 by using the expansion of χ 3 , we have to set D 3 = 5G 3 (see formula (26) and the definition of D r in Proposition 3.3). Indeed, G 3 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)). After having determined the upper bounds Z 3 and log(F The new remainder can be estimated as log(|R (4) | ̺ ) = −41.66110. Since it has increased, the third normalization step can be considered as the optimal one with respect to these parameters and the algorithm stops here.
End of the algorithm. We reconsider the set S (3) because it represents the Hamiltonian at the third normalization step, which is related to the smallest remainder. The value of ̺ 0 that optimizes the escaping time, as defined at the end of subsection 4.1, is ̺ 0 = 0.00008165. We can finally compute the lower bound for the stability time T by using formula (37) and we obtain log T = 24.92920.
Final comments. It is easy to realize that the explicit computation of the expansions for the terms belonging to the Hamiltonians up to the degree R I + 2 does not depend on the evaluation of the upper bounds for higher order polynomials. This remark explains the reason why it is convenient to separate the computer-assisted proof in two different programs, 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 a version that is adapted to the example discussed in the present appendix. Such a file is conceived as a sort of supplementary material 16 with respect to this work, thus, it can be interesting for readers willing to reproduce the computational algorithm in their own codes (or to adapt our ones). For such a purpose, in BirkCAP AppsA-B.zip we have included a program computing exactly the polynomial expansions in the framework of Mathematica, in order to provide a code that is rather easy to read. On the other hand, the explicit expansions reported in this appendix have been produced by using X̺óνoζ for the algebraic manipulations of polynomials with coefficients that are complex numbers expressed as intervals. We think that this can be useful for the initial comparisons with a new code, eventually designed by a reader interested in much more challenging applications with respect to the very simple example discussed in the present appendix.