This is a preprint of: “Flow map parameterization methods for invariant tori in Hamiltonian systems”, Àlex Haro, Josep Maria Mondelo, Commun. Nonlinear Sci. Numer. Simul., vol. 101, 105859:1–34, 2021. DOI: [10.1016/j.cnsns.2021.105859] FLOW MAP PARAMETERIZATION METHODS FOR INVARIANT TORI IN HAMILTONIAN SYSTEMS ALEX HARO AND J.M. MONDELO Abstract. The goal of this paper is to present a methodology for the computation of invariant tori in Hamiltonian systems com- bining flow map methods, parameterization methods, and sym- plectic geometry. While flow map methods reduce the dimension of the tori to be computed by one (avoiding Poincaré maps), pa- rameterization methods reduce the cost of a single step of the de- rived Newton-like method to be proportional to the cost of a FFT. Symplectic properties lead to some magic cancellations that make the methods work. The multiple shooting version of the methods are applied to the computation of invariant tori and their invari- ant bundles around librational equilibrium points of the Restricted Three Body Problem. The invariant bundles are the first order ap- proximations of the corresponding invariant manifolds, commonly known as the whiskers, which are very important in the dynami- cal organization and have important applications in space mission design. Keywords. Invariant tori; parameterization method; KAM the- ory; RTBP; Lissajous orbits. 1. Introduction Hamiltonian systems are frequently found in physical and engineer- ing applications, from where challenging problems continuously emerge, of both theoretical and practical nature. The developement of efficient methods for computing invariant tori carrying quasi-periodic motion is a driving force in the applications of Hamiltonian systems (see e.g. [52] for early references), in areas such as plasma physics, semiclassical quantum theory, accelerator theory, magnetohydrodynamics, oceanog- raphy and, of course, celestial mechanics. While Lagrangian (maximal dimensional) invariant tori are important in stability studies, partially A.H. is supported by the grants PGC2018-100699-B-I00 (MCIU-AEI-FEDER, UE), 2017 SGR 1374 (AGAUR), MSCA 734557 (EU Horizon 2020), and MDM- 2014-0445 (MINECO), and by the NSF under Grant No. 1440140 to found his residence at MSRI in Berkeley, California, during the Fall 2018 semester. J.M. Mondelo has been supported by the MINECO-AEI grants MTM2014-52209- C2-1-P, MTM2016-80117-P, MTM2017-86795-C3-1-P. We would also like to acknowledge the work of the free software community in providing all the software tools we have used, which are listed at the beginning of Section 4. 1 2 ALEX HARO AND J.M. MONDELO hyperbolic invariant tori are also important in studies of diffusion and chaos in Hamiltonian systems. An important problem in astrodynam- ics is the design of station keeping orbits lying on partially hyperbolic invariant tori around collinear libration points in RTBP approxima- tions, for which the stable manifolds are sort of entry lanes, and the unstable manifolds are the exit lanes (see [21] for a survey of early li- bration points missions, and the web pages of space agencies for many newer ones). This will be the guiding problem of this paper to fix a framework. To date, one of the most succesful approaches to compute invari- ant tori falls in the category of numerical Fourier methods, in which parameterizations of tori are given by (truncated) Fourier expansions, and the arising discretized invariance equations (using e.g. colloca- tion) are solved by numerical methods such as Newton’s method (see e.g. [15, 35, 14, 47, 1] for several variants of this approach in dif- ferents contexts). In spite of the relatively simple formulation of the approach, the main practical drawback is what we refer to as the large matrix problem [31]: the computational bottleneck produced in solv- ing the large dimension of the systems of equations at each iteration step, whose time cost is O(N 3 ) and the memory cost is O(N 2 ), where N is the number of Fourier coefficients of the approximations. We emphasize that the number of Fourier coefficients has to do with the dimension of the tori to be computed and their regularity. Remarkably, the numerical Fourier method introduced in [27] mitigates the curse of dimensionality by reducing the dimension of the tori by 1, in looking for invariant tori for time-T flow maps (where 1/T is one of the fre- quencies of the motion on the torus) instead of looking for invariant tori for Poincaré maps. As of now, it is a well-stablished method that has proven to be among the most adequate in computing partially hyper- bolic invariant tori around collinear points in the RTBP, by reducing the problem to computing invariant curves of time-T flow maps (see [3] for a review). As we see, avoiding the large matrix problem is already important for models such as the RTBP, but is crucial when facing the computation of higher dimensional tori either in higher dimensional problems (e.g. non restricted problems), or in the non-autonomous (periodic and quasi-periodic) improvements of RTBP (such as the el- liptic case, the bicircular case, de quasi-bicircular case, or other models that come from three or more body problems). As of now, this still has not been attempted in a systematic manner. The object of this paper is twofold. First, to overcome the large matrix problem in [27] and second, but not less important, to establish a mathematical framework for the analysis of the derived algorithms. This will be performed by changing the discretization strategy and linking the approach to the so called parameterization method for in- variant manifolds, a general strategy for proving existence of invariant COMPUTATION OF INVARIANT TORI 3 manifolds in a constructive way, so that the methods of proof lead to algorithms of computation, and can be applied to many different contexts (see [5, 6, 7] for the foundational papers for invariant mani- folds attached to fixed points and [29] for a review). More specifically, the algorithms presented here are inspired by non-perturbative KAM strategies [17, 18, 24, 37, 42] and by symplectic geometry [28, 29], that are applied to look for invariant tori of flow maps in the spirit of the methodology introduced in [27]. The algorithms consist in performing Newton-like steps on the invariance equations at the functional level (rather than directly at the numerical level). To do so, the geometrical and dynamical properties of the problem lead to the construction of a frame especially chosen in order to make block triangular the lin- earized invariance equations. Using FFT to switch the representations of invariant tori from samples to Fourier coefficients and vice-versa makes the time cost of each step O(N log N ) and the memory cost O(N ), where N is the number of either Fourier coefficients or samples used to represent the tori. This is a significant improvement of these functional Fourier methods with respect numerical Fourier methods mentioned above. See [29] for some benchmarks comparing parame- terization method-like methodologies in the context of invariant tori in skew-product systems [32, 31, 33], and also [9, 12, 36, 39] for other contexts. We will present several algorithms of computation and continuation of invariant tori, including the isoenergetic case (i.e. invariant tori at a fixed energy level). For the sake of simplicity, we will focus on partially hyperbolic invariant tori of dimension n − 1 of n degrees of freedom Hamiltonians, but many of the ideas can be extrapolated to other cases (including Lagrangian tori and partially elliptic invariant tori), even to non-autonomous Hamiltonian systems (e.g. using the common trick of adding extra degrees of freedom). The methods we present not only compute the invariant torus, but also the invariant stable, unstable and center bundles at the same time. The stable and unstable bundles are the first order approximations of the corresponding stable and unsta- ble invariant manifolds, commonly known as the whiskers. The center bundle provides the tangent directions to the normally hyperbolic in- variant cylinder containing the family of tori being computed. We will come to the a priori unexpected realization that, when invariant tori and invariant bundles are computed at the same time, the algorithms are much more efficient than it they are designed for computing invari- ant tori only. This is another improvement of the standard approach, in which the stable and unstable bundles are computed after the torus is computed [27, 40], and avoids additional O(N 3 ) time and O(N 2 ) memory costs. The algorithms implement multiple shooting, in order to cope with the instability that comes from the hyperbolic part. In 4 ALEX HARO AND J.M. MONDELO summary, the computational bottleneck of the flow map parameter- ization methods presented in this paper is no longer the solution of the invariance equations but the (unavoidable) numerical integration needed in order to evaluate the flow maps. This is a task that can be performed easily in parallel. Last but not least, the convergence of the algorithms presented here could be proved using KAM methods, under appropriate non- degeneracy conditions that we make explicit (related to well-known Kolmogorov and isoenergetic conditions) and smallness of the error of the approximate solutions. We do not pursue this analysis here. We refer to [18, 24, 42, 28, 8, 13, 29] for proofs of several KAM results in a posteriori format, based on the parameterization method. See also [22] for a methodology to perform computer-assisted-proofs based on a posteriori format KAM theorems. Following the standard practice in numerical analysis, we have tested the algorithms with well-known computations, as the ones appearing in [27]. We will see that, already in this case, the algorithms are much faster and let one reach unex- plored regions and compute tori that are about to break. Summary of the paper. Section 2 provides some geometrical back- ground and introduces notation for the invariance equations to be solved and the parameterizations to be computed. Section 3 progres- sively introduces the necessary algorithms for the solution of multiple cohomological equations and the computation of frames, in order to perform Newton steps and finally continuation ones. The section ends with important comments on an actual implementation. Section 4 is devoted to the application of the algorithms to the computation of the family of partially hyperbolic KAM tori born from the equilibrium point L1 of the Earth-Moon RTBP. A performance comparison is made with previous large-matrix methodology. Some dynamical and geomet- ric observables are introduced and graphically represented, in order to discuss global properties of the family. The graphical evolution of a few specific sub-families of tori is also shown in order to illustrate the interaction with other families of objects of the center manifold of L1 . Section 5 presents some conclusions. The paper is ended with two ap- pendices. Appendix A deals with the equivalence of the Poincaré map and the flow map methods, and Appendix B provides the proofs of some cancellations (coming from geometrical properties) that are cru- cial for the design of the algorithms, as well as for eventual proofs of their convergence using KAM methods. 2. Setting 2.1. Notations. In this paper we assume that all objects are real an- alytic. COMPUTATION OF INVARIANT TORI 5 Let T` = R` /[0, 1]` be the standard `-torus. With a slight abuse of notation, we identify a function ξ : T` → R with a funcion ξ : R` → R that is 1-periodic in each of its variables, the components of θ = (θ1 , . . . , θ` ). The average of ξ is Z hξi := ξ(θ)dθ. T` We denote the Fourier coefficients of ξ as {ξˆk }k∈Z` , which are given by Z ˆ ξk = ξ(θ)e−i2πkθ dθ. T` Then X ξ(θ) = ξˆk ei2πkθ , k∈Z` P` where kθ := i=1 ki θi and i denotes the imaginary unit. The Fourier P coefficients go to zero exponentially fast when |k| := `i=1 |ki | goes to infinity. 2.2. Symplectic structures and Hamiltonian systems. We as- sume we are given an open set U ⊂ R2n endowed with an exact symplec- tic structure whose matrix representation is an antisymmetric matrix map Ω : U → R2n×2n , which is invertible, and it is given by Ω(z) = Da(z)> − Da(z), where a : U → Rn , for which the transpose a(z)> is the matrix repre- sentation of the action form at the point z ∈ U . The matrix map Ω induces a symplectic product at each z ∈ U . For the sake of simplicity, we will also assume we are given an almost complex structure compat- ible with the symplectic structure, meaning a map J : U → R2n×2n that is involutive (J(z)2 = −I2n ), symplectic (J(z)> Ω(z)J(z) = Ω(z)) and such that the matrix map G : U → R2n×2n defined by G(z) = −Ω(z)J(z) induces an scalar product at each z ∈ U . The prototypical example is the standard symplectic structure, given by   On −In Ω0 (z) = , In On where In is the n × n identity matrix, and On is the n × n zero matrix. In this case,       1 On In On −In In On a0 (z) = 2 z, J0 (z) = , G0 (z) = . −In On In On On In Given a function H : U → R, we obtain an autonomous Hamiltonian system with n degrees of freedom, (1) ż = XH (z) := Ω(z)−1 DH(z)> . 6 ALEX HARO AND J.M. MONDELO We will denote by ϕ : D ⊂ R×U → U the Hamiltonian flow associated to the Hamiltonian vector field XH . We will often use the notation ϕt = ϕ(t; ·). For t fixed, the time-t map ϕt is exact symplectic (in the appropriate domains), meaning that is symplectic, that is Dϕt (z)> Ω(ϕt (z))Dϕt (z) = Ω(z), and, moreover, a(ϕt (z))> Dϕt (z) − a(z)> = Dpt (z) for a certain primitive function pt : U → R, which in fact is given by Z t  pt (z) = a(ϕs (z))> XH (ϕs (z)) − H(ϕs (z)) ds. 0 The exactness property leads to crucial cancellations that enable the existence of invariant tori. Moreover, it is well-known that invariant tori carrying quasi-periodic dynamics have the special geometrical property of being isotropic. 2.3. Isotropic tori and Calabi vectors. Given a parameterization K : Tm → U of an m-dimensional torus K, we define its Calabi vector C(K) as Z > (2) C(K) = a(K(θ))> DK(θ) dθ, Tm where we emphasize the dependence on K. Its components are the Calabi invariants, Z (3) Ci (K) = a(K(θ))> ∂K ∂θi (θ) dθ, Tm q for i = 1, . . . , m. The corresponding radii of K are ri (K) = |Ci π(K)| , for i = 1, . . . , m, which measure the widths of the torus K (w.r.t. the parameterization K). Notice that Calabi invariants (and the radii) are invariant under the flow of a Hamiltonian vector field XH , as we prove in the following lines: Z > C(ϕt (K(θ))) = a(ϕt (K(θ)))> D(ϕt (K(θ))) dθ ZT m = a(ϕt (K(θ)))> Dϕt (K(θ))DK(θ) dθ ZT m Z > = a(K(θ)) DK(θ) dθ + Dpt (K(θ))DK(θ) dθ Tm Tm > = C(K) , where pt is the primitive function of ϕt , and we apply that pt ◦ K is 1-periodic in all its variables, so its differential has zero average. COMPUTATION OF INVARIANT TORI 7 Remark 2.3.1. Given a torus automorphism A : Tm → Tm , where A ∈ GLm (Z) (i.e. A ∈ Zm×m , det A = ±1), the reparametrization K ◦ A : Tm → U of K has Calabi vector C(K ◦ A) = A> C(K). In the previous constructs we have only used the symplectic proper- ties of phase space, but not the possible geometrical properties of the tori. The torus K is isotropic if its parameterization satisfies DK(θ)> Ω(K(θ))DK(θ) = 0, for any θ ∈ Tm . In such a case, we may define, for i = 1, . . . , m, Z 1 (4) Ci (K) = a(K(θ))> ∂K ∂θi (θ) dθi 0 by taking any fixed (θ1 , . . . , θi−1 , θi+1 , . . . , θm ) (hence, giving a genera- tor of the torus). It is not difficult to check that the definition does not depend on such a choice (just compute the derivatives with respect to θj with j 6= i) and, hence, equals the definition given in (3). Remark 2.3.2. In the numerical computation of invariant tori it is use- ful to monitor these geometrical quantities in order to detect shrinking of the tori. These are geometrical observables we use along the com- putations. 2.4. Invariance equations for invariant tori. A parameterization K̂ : Td −→ U of an invariant d-dimensional torus K̂ with frequency vector ω̂ ∈ Rd satisfies the invariance equation  (5) ϕt K̂(θ̂) = K̂(θ̂ + tω̂), for all θ̂ ∈ Td and t ∈ R. The infinitesimal version of (5) is (6) XH (K̂(θ̂)) = DK̂(θ̂)ω̂. The frequency vector ω̂ is assumed to be (at least) non-resonant or ergodic, that is, k̂ · ω̂ 6= 0 for any k̂ ∈ Zd \ {0}. It is well-known that the Hamiltonian H is constant on an invariant torus: H(K̂(θ̂)) = h for all θ̂ ∈ Td , for a certain energy h. Remark 2.4.1. Equation (6) determines K̂ up to a phase: if K̂(θ̂) is a solution, then, for any α̂ ∈ Rd , K̂α̂ (θ̂) := K̂(θ̂+α̂) is also a solution, that parameterizes the same torus K̂. We then have d degrees of freedom in the determinacy of the parameterization K̂. Remark 2.4.2. The frequency vector of K̂ is defined up to a unimodular matrix  ∈ Zd×d (with determinant ±1), so that to the reparameteri- zation K̂ ◦  : Td → U corresponds the frequency vector Â−1 ω̂. In order to reduce the dimension of the parameterization to be com- puted, and also to avoid the use of Poincaré map, we borrow a trick from [27]. By writing ω̂ = T1 (ω, 1), where ω ∈ Rd−1 and T ∈ R, 8 ALEX HARO AND J.M. MONDELO one looks instead for a parameterization K : Td−1 → U of a (d − 1)- dimensional torus K inside the starting one K̂ satisfying  (7) ϕT K(θ) = K(θ + ω). We will refer to T as the period or flying time of the torus K inside K̂, and ω as its rotation vector. From K satisfying (7) we recover the parameterization K̂ satisfying (6) via the flow through  (8) K̂(θ̂) = ϕθd T K(θ − θd ω) , where θ̂ = (θ, θd ) ∈ Td−1 × T. From K̂ we get a K just defining K(θ) = K̂(θ, 0). The rotation vector ω ∈ Rd−1 is non-resonant, meaning that k·ω ∈ / Z, for any k ∈ Zd−1 \ {0}. Notice that H(K(θ)) = h for all θ∈T . d−1 Remark 2.4.3. Equation (7) determines K up to a phase: if K(θ) is a solution, then, for any α ∈ Rd−1 , K(θ + α) is also a solution that parameterizes the same torus K inside K̂. But Equation (7) determines K up to a time translation of K inside K̂: if K(θ) is a solution then, for any ᾱ ∈ R, ϕᾱ (K(θ)) is also a solution. We then again have d degrees of freedom in the determinacy of the parameterization K. Remark 2.4.4. The rotation vector of K is defined up to a unimod- ular matrix A ∈ Z(d−1)×(d−1) , in such a way that the rotation vector corresponding to the reparameterization K ◦ A : Td−1 → U is A−1 ω. In the problem of existence (and computation) of invariant tori K with rotation vector ω one may consider two different cases: • isochronous case: T is fixed, and the torus K and the energy h are the unknows; • isoenergetic case: h is fixed, and the torus K and the flying time T are the unknowns. In summary, the equations to face are:  (9) ϕT K(θ) −K(θ + ω) = 0, Z (10) H(K(θ)) dθ − h = 0, Td−1 in which we fix either T or h accordingly. Notice that this formulation also makes natural to consider either T or h as continuation parameters. Remark 2.4.5. As it is very well-known, using Poincaré map is another way of reducing the dimension of the problem. Both the Poincaré map and the time-T flow map approaches are equivalent (see Appen- dix A). From the theoretical point of view, the time-T approach has the advantage that one deals with exact symplectic mapings in both the isochronous and isoenergetic cases, while in the Poincaré map ap- proach one produces symplectic mappings once one reduces also to a fixed energy level (so one works in the isoenergetic case). From the COMPUTATION OF INVARIANT TORI 9 numerical point of view, both approaches of course involve numerical integration, but computing Poincaré maps and their differentials is a bit more time consuming task (both in terms of coding and execution) than computing time-T flows and their differentials. Moreover, since we are planning to perform multiple shooting methods, the use of mul- tiple Poincaré maps can increase the difficulty of appropriately locating them and the complexity of the algorithms. 2.5. Partially hyperbolic invariant tori of dimension n − 1. In this paper, we focus in algorithms for computing partially hyperbolic invariant tori of dimension n − 1, i.e. d = n − 1, that is, with stable and unstable bundles of rank 1. The algorithms we present are very easy to adapt to Lagrangian tori, that is d = n, and other lower dimensional partially hyperbolic tori, that is d < n. Invariant tori with elliptic di- rections could be also considered, with the aid of parameters. The case we consider appears very often in applications, as the one presented in this paper. In Section 4 we compute partially hyperbolic invariant tori around the Lagrangian points of the Restricted Three Body Problem, that is n = 3, d = 2. Hence, assume that d = n − 1. A bundle Ŵ of rank 1 (with base K̂) parameterized by a map Ŵ : Tn−1 → R2n satisfying the equation  (11) Dϕt K̂(θ̂) Ŵ (θ̂) = etχ Ŵ (θ̂ + tω̂), with χ ∈ R, is invariant under the linearized flow of XH . If χ < 0 then Ŵ = Ŵ s is the stable bundle, and if χ > 0 then Ŵ = Ŵ u is the unstable bundle. In order to decrease the dimension by one, we proceed again by considering time-T maps. With ω̂ = T1 (ω, 1), we look for a parame- terization W : Tn−2 → R2n of a line bundle W of the torus K such that  (12) DϕT K(θ) W (θ) = eT χ W (θ + ω). Again, we recover Ŵ by taking  Ŵ (θ̂) = e−θd T χ Dϕθd T K(θ − θd ω) W (θ − θd ω), where θ̂ = (θ, θd ). Remark 2.5.1. Notice that with this formulation we assume that the line bundle Ŵ is trivial (and oriented). Non-oriented line bundles can fit this formulation through the double-covering trick (see e.g. [33]). 3. Multiple shooting algorithms As we have seen in previous section, the equations to solve involve numerical integration of orbits from points on tori up to a certain T . Since T cannot be chosen to be small, dynamical instability can make 10 ALEX HARO AND J.M. MONDELO the numerical solution of (9), (10) and (12) difficult. In our expe- rience, with values of eT χ of the order of thousands we are not able to observe quadratic convergence of Newton iterates. Also, continua- tion steps become very small. The effects of dynamical instability can be avoided by reducing integration time through the use of multiple shooting. Multiple shooting is classically introduced in the numerical solution of boundary value problems for ordinary differential equations (see e.g. [50]), but it can also be used in the computation of general invariant objects. 3.1. Multiple shooting invariance equations. Instead of looking for a parameterization K0 of a single (n − 2)-dimensional torus K0 inside the (n − 1)-dimensional torus K̂, we look for parameterizations i=0 of m tori, Ki : T {Ki }m−1 → R2n satisfying the equations n−2  ω (13) ϕT /m Ki (θ) − Ki+1 (θ + m ) = 0, for i = 0, . . . , m − 1. In the above equation, and in the following, we assume the i subindex (here in Ki ) is defined modulo m. In particular, Km = K0 . Note that, if {Ki }m−1 i=0 is solution of (13), any Ki is solution of (9). In the following, we will refer to {Ki }m−1i=0 as a multiple torus, and similar notations will be used for other multiple objects (bundles, frames, etc.) There are several ways to explicit the energy level of the torus K̂, and the one we consider here is that of the average of the energy on the first torus of the chain: Z (14) H(K0 (θ)) dθ − h = 0. Tn−2 We will also use multiple shooting for computing the invariant bundle Ŵ, and thus look for λ and a multiple bundle {Wi }m−1 i=0 , where Wi : T n−2 → R , satisfying 2n  ω (15) DϕT /m Ki (θ) Wi (θ) − λWi+1 (θ + m ) = 0, for i = 0, . . . , m − 1. Again, if {Wi }m−1 i=0 is solution of (15), any Wi is solution of (12), with eT χ = λm . Following the philosophy of the parameterization method, we will look for adapted frames in which Newton’s method for the invariant equations (13) and (15) can be performed through a sequence of co- homological equations that are diagonal in Fourier space. In this way, the large matrix systems that appear when doing direct Fourier dis- cretization of invariance equations are avoided, and all the computa- tional effort goes in the numerical integrations necessary to perform the change to the frame and to obtain the new equations, and in the transformations from Fourier space to grid space and vice versa, which are done through FFT. COMPUTATION OF INVARIANT TORI 11 3.2. Multiple cohomological equations. As we have just mentioned, the use of adapted frames is in the basis of our algorithms, but also is in the core of KAM theory. This section is devoted to their application to the analysis of multiple cohomological equations. In the following, we will fix a non-resonant rotation vector ω ∈ R` . (In the context of this paper, ` = d − 1 = n − 2.) 3.2.1. Small divisors equations. The first equation we consider is the following small divisors cohomological equation, (16) ξ(θ) − ξ(θ + ω) = η(θ), where η : T` → R is given and ξ : T` → R is to be found. A necessary condition to solve this equation is that the average of the right hand side is zero: hηi = 0. The solution of (16) is (formally) straightorward in terms of Fourier coefficients: if X X (17) ξ(θ) = ξˆk ei2πkθ , η(θ) = η̂k ei2πkθ , k∈Z` k∈Z` then the solutions of (16) are (formally) given by η̂k ξˆ0 ∈ R arbitrary, ξˆk = for k 6= 0. 1 − ei2πkω Note that the divisors 1 − ei2πkω become arbitrarily small. For analytic η, the convergence of the series ξ is ensured by stronger non-resonance properties of rotation vector ω, such as the so-called Diophantine con- dition: from now on, we assume there exist γ > 0, τ ≥ d such that |kω − q| ≥ γ|k|−τ 1 for all k ∈ Z \ {0}, q ∈ Z. We will denote by Rη(θ) ` the only ξ(θ) solution of (18) ξ(θ) − ξ(θ + ω) = η(θ) − hηi, with zero average. The multiple version of the small divisors cohomological equation i=0 , ηi : T → R, we want to find functions (16) is: given functions {ηi }m−1 ` i=0 , ξi : T → R, satisfying {ξi }m−1 ` ω (19) ξi (θ) − ξi+1 (θ + m ) = ηi (θ), i = 0 . . . , m − 1. As already mentioned, we assume the subindex i to be defined modulo m, so that, in particular ξm = ξ0 . A telescopic sum turns (19) into a small divisors cohomological equation: m−1 X i (20) ηi (θ + m ω) = ξ0 (θ) − ξ0 (θ + ω). i=0 Hence, the necessary condition to solve the multiple equation is the average condition m−1 1 X hηi := hηi i = 0. m i=0 12 ALEX HARO AND J.M. MONDELO Here, and in the following we will use the notation hgi for the mean of the averages of a set of functions {gi }m−1 i=0 . Notice that any solution of (20) is given by m−1 X (21) ξ0 (θ) = ξˆ0,0 + Rηi (θ + i m ω). i=0 Once the average ξˆ0,0 = hξ0 i is fixed, the remaining ξi can be deter- mined from (19) by the relations 1 ξi (θ + m ω) = ξi−1 (θ) − ηi−1 (θ), i = 1 . . . , m − 1. Roundoff propagation is reduced by computing every ξi independently from m−1 X ηj+i (θ + mi ω) = ξj (θ) − ξj (θ + ω). i=0 Then, m−1 X ξj (θ) = ξˆj,0 + Rηj+i (θ + i m ω), i=0 with ξˆj,0 = ξˆj−1,0 − η̂j−1,0 . The number of translations to be done to η0 , . . . , ηm−1 (and thus the coding and run-time overhead of this second approach) is reduced by using m−1 X (22) ξj (θ + j m ω) = ξˆj,0 + Rηj+i (θ + j+i m ω), i=0 instead. Anyway, the overhead of this second approach is negligible against the cost of numerical integration (an orbit per Fourier coeffi- cient, as we will see). 3.2.2. Non-small divisors equations. The other kind of equation we will consider is a non-small divisors cohomological equation, (23) λξ(θ) − µξ(θ + ω) = η(θ), where λ, µ ∈ R, |λ| 6= |µ|, η : T` → R are given and ξ : T` → R is to be found. Its formal solution is straightforward in terms of Fourier coefficients: using the notation of (17), for all k ∈ Z, η̂k ξˆk = . λ − µei2πkω The multiple version of the non-small divisors cohomological equation (23) is: given λ, µ ∈ R, |λ| 6= |µ|, {ηi }m−1 i=0 , ηi : T → R, find {ξi }i=0 ` m−1 satisfying ω  (24) λξi (θ) − µξi+1 θ + m = ηi (θ) COMPUTATION OF INVARIANT TORI 13 for i = 0, . . . , m−1. As before, and also to reduce roundoff propagation, a telescopic sum turns (24) into a non-small divisors cohomological equation for each ξj : m−1 X m m λ ξj (θ) − µ ξj (θ + ω) = µi λm−1−i ηj+i (θ + i m ω). i=0 These m independent equations (that can be solved in parallel) are all of the same type as (23). Also as before, the number of translations to be done to η0 , . . . , ηm−1 is reduced by subtituting θ by θ + mj ω in the previous equation. 3.3. Computation of adapted frames. In the spirit of the parame- terization method, for a multiple torus {Ki }m−1 i=0 and a multiple bundle i=0 , Pi : T {Wi }i=0 , we look for a multiple frame {Pi }m−1 → R2n×2n , m−1 n−2 such that   ω −1  Λ Si (θ) (25) Pi+1 (θ + m ) DϕT /m Ki (θ) Pi (θ) = , Λ−> for i = 0, . . . , m − 1, with    1  In−1 0 Si (θ) 0 (26) Λ= , Si (θ) = 0 λ 0 0 where In−1 is the (n − 1) × (n − 1) identity matrix, Si1 (θ) is a symmetric (n − 1) × (n − 1) matrix (referred to as torsion matrix), and each 0 (and empty block) stands for a zero matrix of the corresponding dimensions. We will see this is possible if {Ki }m−1 m−1 i=0 and {Wi }i=0 are solutions of equations (13) and (15), respectively. In our algorithm, the key point is that (25) is approximately true if {Ki }m−1 m−1 i=0 and {Wi }i=0 are approximate solutions. i=0 , Li : T First, we define the multiple subframe {Li }m−1 → R2n×n , n−2 by  (27) Li (θ) = DKi (θ) XH (Ki (θ)) Wi (θ) , By differentiating (13) with respect to θ, we obtain  ω (28) DϕT /m Ki (θ) DKi (θ) = DKi+1 (θ + m ).  ω  By differentiating ϕT /m ϕt (Ki (θ) = ϕt Ki+1 (θ + m ) with respect to t and taking t = 0,   ω  (29) DϕT /m Ki (θ) XH Ki (θ) = XH Ki+1 (θ + m ) . From (28), (29), (15), and the definition of Li (θ) in equation (27), we have  ω (30) DϕT /m Ki (θ) Li (θ) = Li+1 (θ + m )Λ. 14 ALEX HARO AND J.M. MONDELO It is well-know that symplecticity properties imply that each sub- frame Li (θ) is Lagrangian, meaning that (31) Li (θ)> Ω(Ki (θ))Li (θ) = 0. Remark 3.3.1. We emphasize conditions (30) and (31) also work for matrices of the form   0 A 0 Li (θ) = Li (θ) , 0 b where A ∈ R(n−1)×(n−1) and b ∈ R are constant and invertible. In par- ticular, one can scale frames in order to mitigate possible degeneracies. The goal is now completing each Lagrangian subframe Li to a sym- plectic frame Pi , by juxtaposing a complementary Lagrangian frame Ni . There are several ways to do so (see [29]). Here, with the aid of the compatible almost complex structure J, we define  (32) N̂i (θ) = J Ki (θ) Li (θ)Gi (θ)−1 , Gi (θ) = Li (θ)> G(Ki (θ))Li (θ). Now, the matrix  (33) P̂i (θ) = Li (θ) N̂i (θ) , is symplectic, that is to say P̂i (θ)> Ω(Ki (θ))P̂i (θ) = Ω0 . Hence,   ω −1 Λ Ŝi (θ) (34) P̂i+1 (θ + ) DϕT /m (Ki (θ))P̂i (θ) = m On Λ−> is symplectic (with respect to Ω0 ), where Ŝi (θ) is a n × n matrix given by ω > ω (35) Ŝi (θ) = N̂i+1 (θ + m ) Ω(Ki+1 (θ + m ))DϕT /m (Ki (θ))N̂i (θ). From symplecticity, it follows that Ŝ i (θ)Λ> = ΛŜ i (θ)> , In order to have (25), we perform a new change of frame by consid- ering, for i = 0, . . . , m − 1, matrices   In Bi (θ) (36) Qi (θ) = , In with Bi (θ)> = Bi (θ), so that they are symplectic. Then,   −1 In −Bi (θ) Qi (θ) = . In For the new frame (37) Pi (θ) = P̂i (θ)Qi (θ), COMPUTATION OF INVARIANT TORI 15 we have   ω −1 Λ Si (θ) Pi+1 (θ + ) DϕT /m (Ki (θ))Pi (θ) = , m Λ−> with (38) Si (θ) = ΛBi (θ) + Ŝi (θ) − Bi+1 (θ + ω m )Λ−> . By splitting the matrix Ŝi (θ), in blocks of sizes (n − 1) × (n − 1), (n − 1) × 1, 1 × (n − 2) and 1 × 1, as in  1  Ŝi (θ) Ŝi2 (θ) (39) Ŝi (θ) = , Ŝi3 (θ) Ŝi4 (θ) and using analogous splittings for Si (θ) and Bi (θ), formula (38) reads (40) 1 Si1 (θ) = Ŝi1 (θ) + Bi1 (θ) − Bi+1 (θ + ω m ), (41) Si2 (θ) = Ŝi2 (θ) + Bi2 (θ) − Bi+1 2 (θ + ω m )λ−1 , (42) Si3 (θ) = Ŝi3 (θ) + λBi3 (θ) − Bi+1 3 (θ + ω m ), (43) Si4 (θ) = Ŝi4 (θ) + λBi4 (θ) − Bi+1 4 (θ + ω m )λ−1 . We then take Bi1 (θ) = In−1 for each i = 0, . . . , m − 1, so that Si1 (θ) = Ŝi1 (θ), and look for Bi3 (θ)> = Bi2 (θ), Bi4 (θ) so that Si3 (θ)> = Si2 (θ)λ = 0 and Si4 (θ) = 0. In summary, we need to find {Bi2 (θ)}m−1 i=0 and {Bi4 (θ)}m−1 i=0 such that (44) Bi2 (θ) − Bi+1 2 (θ + ω m )λ−1 = −Ŝi2 (θ), i = 0, . . . , m − 1, (45) λBi4 (θ) − Bi+1 4 (θ + ω m )λ−1 = −Ŝi4 (θ), i = 0, . . . , m − 1. The solution of these multiple shooting cohomological equations has been discussed in Section 3.2. All the previous developments are summarized in the algorithm that follows. Algorithm 3.3.2. Given {Ki }m−1 m−1 i=0 , {Wi }i=0 satisfying (13), (15), com- pute the multiple frame {Pi }m−1 i=0 and the corresponding reduced expres- m−1 sion of {DϕT /m ◦ Ki }i=0 , as given by (25), by following these steps: (1) Compute {Li }m−1i=0 from (27). (2) Compute {N̂i }m−1 m−1 i=0 from (32), in order to get {P̂i }i=0 from (33). (3) Compute {Ŝi }i=0 from (35), and obtain {Si1 }i=0 as {Ŝi1 }m−1 m−1 m−1 i=0 , according to the matrix splitting in (39). (4) Compute {Bi2 }m−1 3 m−1 3 i=0 from (44), and {Bi }i=0 from Bi = (Bi ) . 2 > 4 m−1 (5) Compute {Bi }i=0 from (45). (6) Define {Bi1 }m−1 1 i=0 as Bi = In−1 for i = 0 . . . m − 1. (7) Compute {Pi }i=0 from (37), with {Qi }m−1 m−1 i=0 given by (36). Remark 3.3.3. Observe that step 3 is the only one that requires numer- ical integration. The other operations are diagonal either in Fourier space or in grid space. 16 ALEX HARO AND J.M. MONDELO Remark 3.3.4. We can extend the arguments and reduce the torsion matrices {Si1 }m−1 i=0 to constant coefficients by considering multiple small divisors cohomological equations (40). To do so, we define m−1 1 X 1 (46) S01 = hS i =1 hŜ (θ)i m i=0 i and solve (47) 1 Bi1 (θ) − Bi+1 (θ + ω m ) = S01 − Ŝi1 (θ) i = 0, . . . , m − 1. as we discussed in Section 3.2. With the choice, the torsion matrices are Si1 (θ) = S01 , for all i = 0, . . . , m − 1. Hence, step (6) of Algorithm 3.3.2 can be replaced by (6’) Compute {Bi1 }m−1 1 1 i=0 from (47), with Si = S0 given from (46). 3.4. Description of a Newton step. The goal of this subsection is to develop the formulation necessary to perform Newton steps in the (multiple shooting) invariance equations for the tori, (13), the energy level, (14), and the bundles, (15). This will be done by solving the lin- earization of these equations around a known approximation expressed in the frame (37). We will consider isochronous and isoenergetic cases. 3.4.1. A Newton step on the torus. Let us consider a multiple torus {Ki }m−1 m−1 i=0 , a multiple bundle {Wi }i=0 and λ 6= 0 satisfying equations (13), (14) and (15) approximately, for a given T (and fixed ω). Let i=0 be the error in the invariant equations, so , Ei : T {Ei }m−1 → R2n n−2 is defined by  ω (48) Ei (θ) := ϕT /m Ki (θ) − Ki+1 (θ + m ), for i = 0, . . . , m − 1. We also consider the energy error for a given h: Eh = hH(K0 )i − h. We plan to give rather explicit formulas for the corrections of tori, {∆Ki }m−1 i=0 , and also of the energy, ∆h, and the flying time ∆T , in the isochronous case (for which ∆T = 0) and the isoenergetic case (for which ∆h = 0). With {Pi }m−1 i=0 the frame defined in (37), we write the correction of the tori in the form ∆Ki (θ) = Pi (θ)ξi (θ). Expanding by Taylor up to first order the invariance equation  ϕ(T +∆T )/m Ki (θ) + Pi (θ)ξi (θ) ω ω ω −Ki+1 (θ + m ) − Pi+1 (θ + )ξ (θ m i+1 + m ) =0 around the approximated tori and flying time, and neglectic second order error terms, we get the equation for one step of Newton’s method:  DϕT /m Ki (θ) Pi (θ)ξi (θ) ω  ω  ω  −Pi+1 θ + m ξi+1 θ + m + XH Ki+1 θ + m ∆τ = −Ei (θ), COMPUTATION OF INVARIANT TORI 17 where ∆τ = ∆T m . ω −1 Multiplying the previous equations by Pi+1 (θ + m ) , if the frame m−1 {Pi }i=0 satisfied (25) exactly, we would obtain   Λ Si (θ) ω (49) ξi (θ) − ξi+1 (θ + m ) + e∆τ = ηi (θ) 0 Λ−> with Λ, Si (θ) defined as in (26), ω −1 (50) ηi (θ) = −Pi+1 (θ + m ) Ei (θ), and       ξi1 (θ) en−1 ηi1 (θ) ξi2 (θ)  0  ηi2 (θ) ξi (θ) =   ξi3 (θ) , e=  0 ,  ηi (θ) =   ηi3 (θ) , ξi4 (θ) 0 ηi4 (θ) being en−1 = (0, . . . , 1)> ∈ Rn−1 . Notice we implicitly consider ξi , ηi : Tn−2 → Rn−1 × R × Rn−1 × R and enumerate the corresponding block components accordingly. Actually, since {Ki }m−1 m−1 i=0 , {Wi }i=0 satisfy equations (13), (15) approximately, the frame also satisfies (25) ap- proximately, so we would need to add an error term to the matrix in (49). This error term can be disregarded, because when multiplied by ξi (θ) becomes of second order. Taking this into accout, we can rewrite (49) as a system of equations in order to obtain (51) ξi1 (θ) + Si1 (θ)ξi3 (θ) − ξi+1 1 (θ + ω m )+ en−1 ∆τ = ηi1 (θ), (52) λξi2 (θ) − ξi+1 2 (θ ω + m) = ηi2 (θ), (53) ξi3 (θ) − ξi+1 3 (θ +m ω ) = ηi3 (θ), (54) λ−1 ξi4 (θ) − ξi+1 4 (θ +m ω ) = ηi4 (θ), for i = 0, . . . , m − 1. Equations (52), (54) can be solved as multiple non-small divisors cohomological equations of the type (24). Equation (53) is a multiple small divisors cohomological equation of the type (19). Solving it as such implies assuming that hη 3 i = 0. In Appendix B we prove that exact symplectic properties of the flow imply that hη 3 i is in fact quadratically small, so the previous assump- tion is coherent with a Newton step. We denote {ξ¯i3 }m−1 i=0 to be the ¯3 solution of (53) with hξ0 i = 0 and define ξ 3 (θ) = ξ 3 + ξ¯3 (θ), i 0,0 i 3 for i = 0, . . . , m − 1 with ξ0,0 free. Notice that {ξi3 }m−1 i=0 is a solution of (53) for any ξ0,0 , which will be fixed later. Substituting {ξi3 }m−1 3 i=0 in (51), we obtain (55) ξ 1 (θ) − ξ 1 (θ + ω ) = η 1 (θ) − S 1 (θ)ξ¯3 (θ) − S 1 (θ)ξ 3 − en−1 ∆τ . i i+1 m i i i i 0,0 18 ALEX HARO AND J.M. MONDELO For this last equation to be solved as a multiple shooting small divisors cohomological equation (19), we need the sum of averages of the right hand sides to be zero, which gives m−1 X m−1 X hSi1 iξ0,0 3 + men−1 ∆τ = hηi1 − Si1 ξ¯i3 i, i=0 i=0 3 which is a linear system for and ∆τ equivalent to ξ0,0 hS 1 iξ0,0 3 + en−1 ∆τ = hη 1 − S 1 ξ¯3 i. R Now, considering the error in energy as Eh = Tn−2 H(K0 (θ))dθ − h, by substitution of the corrected first torus and energy in (14) and linearization, it follows that hDH(K0 (θ))P0 (θ)ξ0 (θ) − ∆hi = −Eh . Moreover, by using that the frame P0 (θ) is approximately symplectic,  DH(K0 (θ))P0 (θ) = −XH (K0 (θ))> Ω(K0 (θ))P0 (θ) ' 0> > n−1 0 en−1 0 . Hence, neglecting second order error terms we get (56) e> 3 n−1 ξ0,0 − ∆h = −Eh . 3 By collecting the linear equations for ξ0,0 , ∆τ and ∆h we get the (n − 1) × n system  3   1  ξ0,0  1  hS i en−1 0   hη − S 1 ξ¯3 i ∆τ = e> n−1 0 −1 −Eh ∆h It suffices system matrix has rank n − 1 to get one-parameter families of solutions. We consider two cases: • Isochronous case: ∆τ = 0, and ∆h = Eh + e> 3 3 n−1 ξ0,0 , where ξ0,0 solves linear equation (57) hS 1 iξ 3 = hη 1 − S 1 ξ¯3 i 0,0 provided that dethS 1 i = 6 0. This is the isochronous twist condition, which corresponds to Kolmogorov condition. (In practice, the energy equation is not considered, so ∆h is not computed.) 3 • Isoenergetic case: ∆h = 0, and ξ0,0 , ∆τ solve the linear equation  1  3   1  hS i en−1 ξ0,0 hη − S 1 ξ¯3 i (58) = , e> n−1 0 ∆τ −Eh provided that  1  hS i en−1 det > 6= 0. en−1 0 This is the isoenergetic twist condition. COMPUTATION OF INVARIANT TORI 19 Once one computes ξ0,0 3 , {ξ3i }m−1 i=0 is fully determined and one can solve (55). The general solution is (59) ξi1 (θ) = ξ0,0 1 + ξ¯i1 (θ), 1 for i = 0, . . . , m − 1, where ξ0,0 is free and {ξ¯i1 }m−1 i=0 is the solution of ¯1 1 (55) with hξ0 i = 0. The freedom to choose ξ0,0 has to do with the phase and time underterminacy of the parameterization of the first (and then 1 all) tori (see Remark 2.4.3). A simple choice is ξ0,0 = 0. All the process described is summarized in the algorithm that follows. Algorithm 3.4.1. (Newton step on a multiple torus, isochronous or m−1 isoenergetic case) Let {Ki }i=0 , {Wi }m−1 i=0 , λ satisfy equations (13), (14), (15) approximately. Obtain the corrected tori and, in the isoenergetic case, also the corrected flying time by following these steps: (1) Compute {Pi }m−1 1 m−1 i=0 , {Si }i=0 , following Algorithm 3.3.2. (2) Compute the multiple error {Ei }m−1 i=0 from (48). (3) Compute the right-hand side of the cohomological equations {ηi }m−1 i=0 from (50). (4) Solve (52), (54) as non-small divisors multiple cohomological equations, in order to obtain {ξi2 }m−1 4 m−1 i=0 , {ξi }i=0 . (5) Solve (53) as small divisors multiple cohomological equation, in order to obtain its zero-average solution {ξ¯i3 }m−1 i=0 . (6) In the isochronous (resp. isoenergetic) case, compute hS 1 i and the right-hand side of the linear system (57) (resp. (58)) and 3 3 solve it in order to obtain ξ0,0 (resp. ξ0,0 , ∆τ ). (7) Solve (51) as small divisors multiple cohomological equation in order to obtain {ξ¯i1 }m−1 1 m−1 i=0 , and obtain {ξi }i=0 from (59) by 1 choosing ξ0,0 = 0. (8) Compute the corrected tori as {Ki (θ) + Pi (θ)ξi (θ)}m−1 i=0 . In the isoenergetic case, obtain also the corrected flying time as T + m∆τ . 3.4.2. A Newton step on the bundle. Consider again parameterizations {Ki }m−1 i=0 of (n − 2)-dimensional tori inside a larger (n − 1)-dimensional torus, and associated parameterizations {Wi }m−1 i=0 of (n−2)-dimensional bundles inside the bundle of the larger torus, that satisfy equations (13) and (15) approximately. (In the implementation, the Ki ’s are the ones we have improved in the previous step.) Denote the error in the invariant equations of the Wi as  (60) EiW (θ) := DϕT /m Ki (θ) Wi (θ) − λWi+1 (θ + m ω ), for i = 0, . . . , m − 1. Recalling the frame {Pi }m−1 i=0 defined in (37), we would like to find corrections Pi (θ)ξi (θ) of the bundles Wi (θ) and a correction ∆λ for the eigenvalue λ such that the corrected bundles Wi (θ) + Pi (θ)ξi (θ) and corrected eigenvalue λ + δ make (15) vanish 20 ALEX HARO AND J.M. MONDELO at first order. The invariance equation on the corrected bundles and eigenvalue is   DϕT /m Ki (θ) W i (θ) + Pi (θ)ξi (θ) ω ω ω  − (λ + ∆λ) Wi+1 (θ + m ) + Pi+1 (θ + m )ξi+1 (θ + m ) = 0. Expanding the parentheses and neglecting errors of second order, as ω the ones with a factor ∆λ ξi+1 (θ + m ), the previous equation becomes  DϕT /m Ki (θ) Pi (θ)ξi (θ) − λPi+1 (θ + ω )ξ (θ + m m i+1 ω ω ) − ∆λ Wi+1 (θ + m ) = −EiW (θ). ω −1 Multiplying by Pi+1 (θ +m ) , using (25) and neglecting second order error terms, we obtain   Λ Si (θ) ω (61) ξi (θ) − λξi+1 (θ + ) − e∆λ = ηi (θ), 0 Λ−> m with Λ, Si (θ) defined as in (26), ω −1 W (62) ηi (θ) = −Pi+1 (θ + m ) Ei (θ), and       ξi1 (θ) 0 ηi1 (θ) ξi2 (θ) 1 ηi2 (θ) ξi (θ) =   ξi3 (θ) , e=  0 , ηi (θ) =   ηi3 (θ) . ξi4 (θ) 0 ηi4 (θ) As before, we implicitly consider ξi , ηi : Tn−2 → Rn−1 ×R×Rn−1 ×R and enumerate the corresponding block components accordingly. Rewriting (61) as a system of equations, we obtain (63) ξi1 (θ) + Si1 (θ)ξi3 (θ) − λξi+1 1 (θ + ω m ) = ηi1 (θ), (64) λξi2 (θ) − λξi+1 2 (θ + ω m ) − ∆λ = ηi2 (θ), (65) ξi3 (θ) − λξi+1 3 (θ + ω m ) = ηi3 (θ), (66) λ−1 ξi4 (θ) − λξi+1 4 (θ + ω m ) = ηi4 (θ), for i = 0, . . . , m − 1. Equations (65) and (66) can be solved as multi- ple non-small divisors cohomological equations of the form (24). Once ξi3 (θ) is known, (63) is also solved as multiple non-small divisors co- homological equation. For (64) to be solved as multiple small divisors cohomological equation of the form (19), we need that m−1 X hηi2 (θ) + ∆λi = 0, i=0 which is achieved by taking ∆λ = −hη 2 i. In doing so, hξ02 i remains free. In this case, the freedom to choose this average is related to the underdeterminacy of selecting the lengths of {Wi }m−1 i=0 (as, say, the undeterminacy of selecting the lenght of an eigenvector of a given eigenvalue of a matrix). A simple choice is to take hξ02 i = 0. COMPUTATION OF INVARIANT TORI 21 The previous steps are summarized in algorithm that follows. Algorithm 3.4.2. (Newton step on a multiple bundle) Let {Ki }m−1 i=0 , {Wi }m−1 i=0 , λ satisfy equations (13), (14), (15) approximately. Obtain the corrected bundle and eigenvalue by following these steps: (1) Compute {Pi }m−1 1 m−1 i=0 , {Si }i=0 , following Algorithm 3.3.2. (2) Compute the multiple error {EiW }m−1 i=0 from (60). (3) Compute the right-hand side of the cohomological equations {ηi }m−1 i=0 from (62). (4) Solve (65), (66), (63) for {ξi3 }m−1 4 m−1 1 m−1 i=0 , {ξi }i=0 , {ξi }i=0 , respec- tively, as multiple non-small divisors cohomological equations. (5) Take ∆λ = −hη 2 i. (6) Take {ξi2 }m−1 2 i=0 as the solution with hξ0 i = 0 of (64) as small divisors multiple cohomological equation. (7) Compute the corrected multiple bundle as {Wi (θ)+Pi (θ)ξi (θ)}m−1 i=0 and the corrected eigenvalue as λ + ∆λ. 3.5. Algorithms for continuation. In this section we explain a me- thodology to continue invariant tori with respect to parameters. We will first consider continuation with respect to time T (the flying time, related with one of the frequencies of the torus) and with respect to the energy h. The methodology can easily be adapted for continuation of external parameters (appearing on the Hamiltonian). In fact, we will only present an algorithm to compute the tangent of the continuation curve with respect to the continuation parameter. This provides a first order aproximation for the seed for the new value, as it is common practice in numerical continuation (see e.g. [2]). Note that both approaches produce the same family of invariant tori, this is, the one that corresponds to the fixed value of the rotation vector ω. The difference is that the first approach is better to aim to a torus with a specific value of return time T , whereas the second is better to aim to a specific energy h. Aiming to a specific energy is useful to produce isoenergetic Poincaré sections, which is a common way to represent the center manifold of the collinear points of the RTBP (see e.g. [41, 26, 27]). In this respect, we also consider at the end of this section the continuation of the objects with respect to the rotation vector ω, in the isoenergetic case. Notice that in this case the family of objects is parameterized by a Cantor set of parameters, the rotation vector, and the derivatives to be computed are in Whitney sense. 3.5.1. Continuation with respect to T . In order to be able to perform continuation with respect to T , our goal now is to compute the deriva- tives with respect to T of the parameterizations of tori and bundles that solve equations (13) and (15). Assume that these equations define implicitly {Ki }m−1 m−1 i=0 , {Wi }i=0 as functions of T . In order not to burden the notation, we do not write the dependence on T of Ki , Wi , λ, but 22 ALEX HARO AND J.M. MONDELO we denote by ∂T Ki (θ), ∂T Wi (θ), ∂T λ their corresponding derivatives with respect to T . By differentiating (13) with respect to T we obtain  DϕT /m Ki (θ) ∂T Ki (θ) ω  − ∂T Ki+1 (θ + m ) + m1 XH Ki+1 (θ + m ω ) = 0. Considering the derivatives ∂T Ki (θ) in the frame, this is, assuming ∂T Ki (θ) = Pi (θ)ξi (θ), the previous equation is rewritten as  DϕT /m Ki (θ) Pi (θ)ξi (θ) ω ω  − Pi+1 (θ + m )ξi+1 (θ + m ) = − m1 XH Ki+1 (θ + m ω ) . ω −1 Multiplying by Pi+1 (θ + m ) , we obtain equation (49) but without the e∆τ term and with a different definition for ηi , namely   en−1   0  ω −1 ηi (θ) = − m1 Pi+1 (θ + m ) XH Ki+1 (θ + mω ) = − m1   0n−1  , 0 for i = 0, . . . , m − 1. It can be solved as described in subsection 3.4.1, under isochronous nondegeneracy condition. By differentiating (15) with respect to T , we obtain     ∂ ∂T Dϕ T /m K i (θ) W i (θ) + DϕT /m K i (θ) ∂T Wi (θ) ω ω − (∂T λ)Wi+1 (θ + m ) − λ ∂T Wi+1 (θ + m ) = 0. The above expression is rewritten as  DϕT /m Ki (θ) ∂T Wi (θ) − λ ∂T Wi+1 (θ + ω m ) − ∂T λ Wi+1 (θ + ω m ) = −EiW (θ) with     EiW (θ) = ∂ ∂T DϕT /m Ki (θ) Wi (θ)  (67) =m λ DXH Ki+1 (θ + m ω ) Wi+1 (θ + mω ) h i + D2 ϕT /m Ki (θ) Wi (θ), ∂T Ki (θ) .   Here D2 ϕT Ki (θ) ·, · is the bilinear form given by the second dif- ferential of ϕT evaluated at Ki (θ). These equations are of the same type we have considered in Section 3.4.2, by taking frames and writing ∂T Wi (θ) = Pi (θ)ξi (θ), ∂T λ = ∆λ. The steps to follow in order to solve the two systems of multiple cohomological equations from which {∂T Ki }m−1 m−1 i=0 , {∂T Wi }i=0 , ∂T λ can be obtained are summarized in the algorithm that follows. In order not to burden the notation, {ξi }m−1i=0 is used to denote the solution of both systems. COMPUTATION OF INVARIANT TORI 23 Algorithm 3.5.1. (Continuation step with respect to T ) Let {Ki }m−1 i=0 , m−1 {Wi }i=0 , λ be implicit functions of T through equations (13), (14), (15). Find {∂T Ki }m−1 m−1 i=0 , {∂T Wi }i=0 , ∂T λ through the following steps: (1) Compute {Pi }m−1 1 m−1 i=0 , {Si }i=0 , following Algorithm 3.3.2. (2) Take ξi2 (θ) = ξi4 (θ) = 0, i = 0, . . . , m − 1. (3) Compute hS 1 i and ξ0,0 3 = − m1 hS 1 i−1 en−1 . (4) Find {ξi1 }m−1 1 i=0 as the solution with hξ0 i = 0 of the small divisors multiple cohomological equation 1 ξi1 (θ) − ξi+1 (θ + ω m ) 3 = − m1 en−1 − Si1 (θ)ξ0,0 . (5) Obtain {∂T Ki }m−1 i=0 as ∂T Ki (θ) = Pi (θ)ξi (θ). (6) Compute {EiW }m−1i=0 from (67). m−1 (7) Compute {ηi }i=0 from (62). (8) Solve (65), (66), (63) for {ξi3 }m−1 4 m−1 1 m−1 i=0 , {ξi }i=0 , {ξi }i=0 , respec- tively, as multiple non-small divisors cohomological equations. (9) Take ∂T λ = −hη 2 i. (10) Take {ξi2 }m−1 2 i=0 as the solution with hξ0 i = 0 of (64) as small divisors multiple cohomological equation. (11) Compute {∂T Wi }m−1i=0 as ∂T Wi (θ) = Pi (θ)ξi (θ). Remark 3.5.2. Step (2) of Algorithm 3.5.1 (see also later Algorithm 3.5.3 and Algorithm 3.5.4), set to zero the components of the derivatives of the parameterizations of the tori in the hyperbolic directions. This is geometrically very natural, since partially hyperbolic tori are presented in families (for instance contained in a center manifold of an equilib- rium point, or in a normally hyperbolic cylinder), and the derivatives are tangent to such families (and that invariant objects). 3.5.2. Continuation with respect to h. In this section we will consider the continuation of invariant tori with respect to the energy h, so we have to compute the derivatives with respect to h of the parameteriza- tions of tori, bundles and flying time T that solve equations (13), (14), (15). Assume that these equations define implicitly T , {Ki (θ)}m−1 i=0 , m−1 {Wi (θ)}i=0 and λ as functions of h. We will follow the same criterion in the notation as in the previous section, omiting explicitly the depen- dence on h of these objects, but we denote by ∂h T , ∂h Ki (θ), ∂h Wi (θ), ∂h λ their corresponding derivatives with respect to h. By differentiating (13) and (14) with respect to h we obtain  DϕT /m Ki (θ) ∂h Ki (θ) ω  − ∂h Ki+1 (θ + m ) + m1 XH Ki+1 (θ + m ω ) ∂h T = 0 and hDH(K0 (θ))∂h K0 (θ)i − 1 = 0. 24 ALEX HARO AND J.M. MONDELO Considering the derivatives ∂h Ki (θ) in the frame, this is, assuming ∂h Ki (θ) = Pi (θ)ξi (θ), and denoting ∂h τ = m1 ∂h T , the previous equa- tions are rewritten as  DϕT /m Ki (θ) Pi (θ)ξi (θ) ω ω ω  − Pi+1 (θ + m )ξi+1 (θ + m ) + XH Ki+1 (θ + m ) ∂h τ = 0. and hDH(K0 (θ))P0 (θ)ξ0 (θ) − 1i = 0. ω −1 Multiplying by Pi+1 (θ + m ) , we obtain equation (49) with ∆τ = ∂h τ and ηi = 0. We also obtain equation (56) but with 1 in place of ∆h and 0 in place of Eh , that is e> 3 n−1 ξ0,0 = 1. This can be solved as described in subsection 3.4.1, under isoenergetic nondegeneracy condition. By differentiating (15) with respect to h, we obtain  DϕT /m Ki (θ) ∂h Wi (θ) − λ ∂h Wi+1 (θ + ω m ) − ∂h λ Wi+1 (θ + ω m ) = −EiW (θ) with     EiW (θ) = ∂ DϕT /m Ki (θ) ∂h Wi (θ) (68) λ ω  ω =m DXH Ki+1 (θ + m ) Wi+1 (θ + m )∂h T 2   + D ϕT /m Ki (θ) Wi (θ), ∂h Ki (θ) These equations are of the same type we have considered in Sec- tion 3.4.2, by taking frames and writing ∂h Wi (θ) = Pi (θ)ξi (θ), ∂h λ = ∆λ. The steps to follow in order to solve the two systems of multiple cohomological equations from which {∂h Ki }m−1 m−1 i=0 , ∂h T , {∂h Wi }i=0 , ∂h λ can be obtained are summarized in the algorithm that follows. As in Section 3.5.1, we use {ξi }m−1 i=0 to denote the solution of both systems. Algorithm 3.5.3. (Continuation step with respect to h) Let {Ki }m−1 i=0 , m−1 T , {Wi }i=0 , λ be implicit functions of h through equations (13), (14), (15). Find {∂h Ki }m−1 m−1 i=0 , ∂h T , {∂h Wi }i=0 , ∂h λ through the following steps: (1) Compute {Pi }m−1 1 m−1 i=0 , {Si }i=0 , following Algorithm 3.3.2. (2) Set ξi2 (θ) = ξi4 (θ) = 0 and ξi3 (θ) = ξ0,0 3 for i = 0, . . . , m − 1, where  3   1 −1   ξ0,0 hS i en−1 0n−1 = > ∂h τ en−1 0 1 (3) Take {ξi1 }m−1 i=0 as the solution of 1 ξi1 (θ) − ξi+1 (θ + ω m ) = −en−1 ∂h τ − Si1 (θ)ξ0,0 3 , with hξ01 i = 0, and set ∂h T = m ∂h τ . (4) Obtain {∂h Ki }m−1 i=0 as ∂h Ki (θ) = Pi (θ)ξi (θ). (5) Compute {EiW }m−1 i=0 from (68). COMPUTATION OF INVARIANT TORI 25 (6) Compute {ηi }m−1 i=0 from (62). (7) Solve (65), (66), (63) for {ξi3 }m−1 4 m−1 1 m−1 i=0 , {ξi }i=0 , {ξi }i=0 , respec- tively, as multiple non-small divisors cohomological equations. (8) Take ∂h λ = −hη 2 i. (9) Take {ξi2 }m−1 2 i=0 as the solution with hξ0 i = 0 of (64) as small divisors multiple cohomological equation. (10) Compute {∂h Wi }m−1i=0 as ∂h Wi (θ) = Pi (θ)ξi (θ). 3.5.3. Continuation with respect to ω, in the isoenergetic case. Another natural continuation problem is, given a certain energy level h, continue the invariant tori on such an energy level and their invariant bundles with respect to the frequencies. Assume that equations (13), (14), (15) define implicitly T , {Ki (θ)}m−1 m−1 i=0 , {Wi (θ)}i=0 and λ as functions of ω. We will again omit explicitly the dependence on the rotation vector ω of these objects, but we denote by ∂ω T , ∂ω Ki (θ), ∂ω Wi (θ), ∂ω λ their corresponding derivatives with respect to ω. Since the domain of the rotation vector is a Cantor set, these derivatives are understood formally (in fact, these are derivatives in the sense of Whitney, but we will avoid technicalities here). By differentiating (13) and (14) with respect to ω we obtain  DϕT /m Ki (θ) ∂ω Ki (θ) ω  − ∂ω Ki+1 (θ + m ) + m1 XH Ki+1 (θ + mω ) ∂ω T = m1 DKi+1 (θ + m ω ) and hDH(K0 (θ))∂ω K0 (θ)i = 0. Considering ∂ω Ki (θ) = Pi (θ)ξi (θ), denoting ∂ω τ = m1 ∂ω T , and mul- ω −1 tiplying by Pi+1 (θ + m ) , we obtain equation (49) with ∆τ = ∂ω τ and   1 In−2 ηi (θ) = m . O Moreover, equation (56) reads e> 3 n−1 ξ0,0 = 0. This can be solved as described in subsection 3.4.1, under isoenergetic nondegeneracy condi- tion. By differentiating (15) with respect to ω, we obtain  DϕT /m Ki (θ) ∂ω Wi (θ) − λ ∂ω Wi+1 (θ + ω m ) − ∂ω λ Wi+1 (θ + ω m ) = −EiW (θ) with     EiW (θ) = ∂ DϕT /m Ki (θ) ∂ω Wi (θ) λ ω  ω (69) =m DXH Ki+1 (θ + m ) Wi+1 (θ + m )∂ω T   + D2 ϕT /m Ki (θ) Wi (θ), ∂ω Ki (θ) λ ω − m DWi+1 (θ + m ). 26 ALEX HARO AND J.M. MONDELO These equations are of the same type we have considered in Section 3.4.2, by taking frames and writing ∂ω Wi (θ) = Pi (θ)ξi (θ), ∂ω λ = ∆λ. Similarly as we proceed in previous sections, we present in an al- gorithm the steps to solve the two systems of multiple cohomological equations to compute {∂ω Ki }m−1 m−1 i=0 , ∂ω T , {∂ω Wi }i=0 , ∂ω λ. Notice that these are equations for the partial derivatives with respect to the com- ponents of ω, that we formulate in paralell. Algorithm 3.5.4. (Continuation step with respect to ω, for energy h fixed) Let {Ki }m−1 m−1 i=0 , T , {Wi }i=0 , λ be implicit functions of h through equations (13), (14), (15). Find {∂ω Ki }m−1 m−1 i=0 , ∂ω T , {∂h Wi }i=0 , ∂ω λ through the following steps: (1) Compute {Pi }m−1 1 m−1 i=0 , {Si }i=0 , following Algorithm 3.3.2. (2) Set ξi2 (θ) = ξi4 (θ) = 0 and ξi3 (θ) = ξ0,0 3 for i = 0, . . . , m − 1, where    3   1 −1 1 In−2 ξ0,0 hS i en−1 m  0  = > ∂ω τ en−1 0 0 (3) Take {ξi1 }m−1 i=0 as the solution of   1 1 In−2 ω 1 ξi (θ) − ξi+1 (θ + m ) = m − en−1 ∂ω τ − Si1 (θ)ξ0,0 3 , 0 with hξ01 i = 0, and set ∂ω T = m ∂ω τ . (4) Obtain {∂ω Ki }m−1 i=0 as ∂ω Ki (θ) = Pi (θ)ξi (θ). (5) Compute {EiW }m−1 i=0 from (69). m−1 (6) Compute {ηi }i=0 from (62). (7) Solve (65), (66), (63) for {ξi3 }m−1 4 m−1 1 m−1 i=0 , {ξi }i=0 , {ξi }i=0 , respec- tively, as multiple non-small divisors cohomological equations. (8) Take ∂ω λ = −hη 2 i. (9) Take {ξi2 }m−1 2 i=0 as the solution with hξ0 i = 0 of (64) as small divisors multiple cohomological equation. (10) Compute {∂ω Wi }m−1 i=0 as ∂ω Wi (θ) = Pi (θ)ξi (θ). Remark 3.5.5. In the implementation of the continuation with respect to rotation vector, one have to select continuation steps for which the rotation vectors are Diophantine. In practice, since in the computer the rotation vectors are rational, these has to be selected as resonant but of very high order. The continuation can run into troubles when finding strong resonances, since these are more difficult to jump. 3.6. Some comments about the implementation. As it is com- mon in implementation of the parameterization method in KAM-like contexts (see [29] for an overview), the implementation of all the previ- ous algorithms relies on two numerical representations of all the func- tions ζ : Td → R involved. In the grid representation, the function is represented as a set of its values in a uniform grid of Td . In the COMPUTATION OF INVARIANT TORI 27 Fourier representation, the function is represented as a set of approx- imate Fourier coefficients. Through the Discrete Fourier Transform (DFT), that provides approximations of the Fourier coefficients, one representation can be converted into the other. For an easier exposi- tion, we will assume d = 1 in this section (which is the case in all the numerical exploration of Section 4). All the arguments generalize to d > 1, although an actual implementation is subtle (see [29] for com- ments). For detailed expositions on the DFT and its applications, see e.g. [25, 4, 46]. Let ζ : T → R be a function. Choose N > 0. Its grid representation N −1 is given by {ζj }j=0 , with ζj = ζ(j/N ). Its Fourier representation is given by a finite set of (complex) approximate Fourier coefficients, [N/2] {ζ̃k }k=0 , where [·] denotes integer part. The two representations are related by the DFT, N −1 1 X −i2πk j ζ̃k = ζj e N. N j=0 −1 N −1 The DFT is a linear, one-to-one map between {ζj }N j=0 and {ζ̃k }k=0 . Namely, the expression for the inverse DFT is N X −1 j ζj = ζ̃k ei2πk N . k=0 The DFT and its inverse are efficiently evaluated through a family of algorithms known as Fast Fourier Transform (FFT), which allow −1 N −1 computing {ζ̃k }N k=0 from {ζj }j=0 in O(N log N ) operations. The DFT is N -periodic in k, and, for real ζ (which is our case), satisifes the Hermitian symmetry. This is, for k ∈ Z, ζ̃k = ζ̃k+N , ζ̃−k = (ζ̃k )∗ , where ∗ denotes complex conjugate. This gives rise to redundancy in N −1 {ζ̃k }k=0 , which is eliminated by truncating the Fourier representation at k = [N/2], as we have done. Another consequence of the Hermitian symmetry is that ζ̃0 is real and, if N is even, ζ̃N/2 is also real. The DFT coefficients and the Fourier coefficients are related by ∞ X ζ̃k = ζ̂k + (ζ̂k−lN + ζ̂k+lN ). l=1 Actual bounds of the difference ζ̃k − ζ̂k can be obtained from the Cauchy estimates, that ensure exponential decay in |k| of |ζ̂k | for analytic ζ (see e.g. [34]). This is the basis for a detailed error analysis of the algorithms, that we have not pursued here (see e.g. [22]), but a direct consequence of this approximation is that, for the algorithms to work, N has to be large enough for |ζ̃k − ζ̂k | to be small. In the computations of Section 4, the minimum value of N used is 32. Another consequence 28 ALEX HARO AND J.M. MONDELO is that the relative error of the DFT approximation of the Fourier coefficients increases with k. Even with (large) DFT queues of the order of the machine epsilon, we have observed instability in Newton iterates (divergence after apparent convergence). In our implementation, we prevent it by setting to zero the upper half of the DFT coefficients after each Newton iteration. All the steps of algorithms 3.3.2, 3.4.1, 3.4.2, 3.5.1, 3.5.3 can be done in O(N ) operations in either grid or Fourier representation. Consider, for instance, Algorithm 3.3.2. In the evaluation of (27), Wi (θ) is a function we already have, XH (Ki (θ)) is obtained from Ki (θ) in O(N ) operations in grid form, as −1 {XH (Ki (j/N ))}N j=0 , and DKi (θ) obtained in O(N ) operations in Fourier form from Ki (θ) as gi )k = i2πk(K (DK fi )k . Note that this last expression is actually an approximation: the equality is satisfied by the Fourier coefficients of Ki and DKi , not the DFT ones. As another example, the evaluation of (35) is done in O(N ) operations ω ω in grid form if N̂i+1 (θ + m ), Ki+1 (θ + m ) and DϕT /m (Ki (θ))N̂i (θ) are ω known in grid form. The computation of Ki+1 (θ + m ) from Ki+1 (θ) is done in O(N ) operations in Fourier form (again using approximate identities). The computation of DϕT /m (Ki (θ))N̂i (θ) from Ki (θ) and N̂i (θ) in grid form requires the numerical integration of the differential equations (1) with his first variationals applied to several vectors on N trajectories. The different steps of the algorithms stated consist of evaluating equations like the ones just mentioned and solving multiple cohomo- logical equations. The solution of multiple cohomological equations discussed in Section 3.2 is also done in O(N ) operations in Fourier form. As mentioned, converting from grid to Fourier and vice-versa us- ing FFT requires O(N log N ) operations. For the values of N for which numerical integration is feasible, log N is small enough for the cost of FFT to be considered O(N ). Numerical integration is present in step 3 of Algorithm 3.3.2 (Eq. (35)), step 6 of Algorithm 3.5.1 (Eq. (67)), step 5 of Algorithm 3.5.3 (Eq. (68)) and step 5 of Algorithm 3.5.4 (Eq.(69)). Note that Eqs. (67), (68), (69) require numerical integra- tion of the second variational equations. Note also that all algorithms perform numerical integration when computing the frame through Al- gorithm 3.3.2. The cost of numerical integration is formally also O(N ), but for re- alistic estimates it needs to be considered separatedly. Consider for instance Equation (35) as discussed above. In the computations of Section 4, n = 3, so N̂i (θ) has 3 columns, and the system of differential equations to be numerically integrated consists of 24 equations. Taking COMPUTATION OF INVARIANT TORI 29 for example the Runge-Kutta-Felhberg method of orders 7 and 8 for numerical integration, that evaluates the vector field 13 times, since orbits in Section 4 take over 50 integration steps, the factor multiply- ing N is already 15600 times the number of operations required by the evaluation of one differential equation. It turns out that, compared to numerical integration, the remaining computational cost is almost neg- ligible. Numerical integration can be parallelized in a straightforward manner by distributing the N trajectories to be integrated among the threads or processes available. An actual implementation of a continuation procedure in order to follow a family of tori by keeping ω constant requires an strategy in order to choose the number of samples N and control the continua- tion step size as we go along the family. We end this section with a proposal of such a strategy in Algorithm 3.6.1, that has worked well in the numerical computations presented in Section 4. For shortness, we will represent a multiple torus, its multiple bundle, flying time and eigenvalue as T = (K0 , . . . , Km−1 , W0 , . . . , Wm−1 , T, λ). For step size control, the following norm of this compound object is considered:  m−1 X 1/2 2 2 2 2 kT k = T + λ + hkKi k2 i + hkWi k2 i , i=0 where kKi k2 stands for the function θ 7→ kKi (θ)k2 , kWi k2 is interpteted analogously, and the averages are approximated as discrete averages of the grid values (i.e as the 0-th DFT coefficient of the function averaged). The errors in the torus and the bundle are estimated as err(T ) = max kEi (j/N )k∞ , errW (T ) = max kEiW (j/N )k∞ , 0≤i ε1 , double N and perform Newton steps (Algo- rithms 3.4.1, 3.4.2) until err(T 0 ) ≤ ε1 . If unsuccessful, half α and restart the algorithm. (6) Set T ← T 0 (i.e., accept the new torus). Perform continuation step size control as α ← αndes /nit . In step 1, Algorithm 3.5.3 can be used instead, by setting δ ← (∂h K0 , . . . , ∂h Km−1 , ∂h W0 , . . . , ∂h Wm−1 , ∂h T, ∂h λ). 4. An application: computation of the Lissajous family of tori in the Restricted Three Body Problem In this section we apply the algorithms described to the computa- tion of the partially hyperbolic tori that emerge from the L1 equilibrium point of the circular, Spatial Restricted Three-body Problem (RTBP), in the Earth-Moon case. This family of tori is known as the Lissajous family by the astrodynamics community, and plays a fundamental role in libration point dynamics. An outgrowth of our algorithms is the simultaneous computation of the stable, unstable and center bundles of the invariant tori and other observables, such as the Lyapunov mul- tipliers and Calabi invariant, that provide geometrical and dynamical information during the computation of the families of tori. In par- ticular, we will obtain information about the quality of hyperbolicity properties, that is useful to detect possible bifurcations and breakdown phenomena [10, 12, 31]. Apart from the intrinsic importance of the example, our choice is motivated by the fact that the Lissajous family has been extensively described in [27], and has been used as a testbed for different algorithms by other authors (see e.g. [3]). Hence, one of our goals is to compare the performance of the algorithms described here with the ones used in [27], on a more thorough numerical exploration of this family. All the numerical explorations presented here have been done by a program that follows a family of tori with constant ω by performing continuation steps through Algorithm 3.6.1 up to a maximum given number, plus additional stopping criteria that will be specified below. The explorations have been carried out on a Fujitsu Celsius R940 work- station, with two 8-core Intel Xeon E5-2630v3 processors at 2.40GHz, running Debian GNU/Linux 9.11 with the Xfce 4.12 desktop. The source code has been written in C, compiled with GCC 6.3.0 and linked against the Glibc 2.24, LAPACK 3.7.0, FFTW 3.3.5 and PGPLOT 5.2.2 libraries. The code uses OpenMP 4.0 extensions in order to per- form simultaneous numerical integrations in the continuation of one COMPUTATION OF INVARIANT TORI 31 family and also to perform numerical continuation of several families at once. The figures have been generated with gnuplot 5.2. 4.1. The Lissajous family of tori of the RTBP. The circular, spatial Restricted Three-Body Problem (RTBP) describes the motion of a particle of infinitesimal mass under the attraction of two massive bodies known as primaries, with masses m1 > m2 > 0. The primaries are assumed to revolve uniformly in circles around their common center of mass. In a rotating system of reference with the primaries in the horizontal coordinate plane, known in astronomical terms as synodic, the primaries can be made to lie at fixed positions in the x1 axis. After a rescaling in space and time, and defining the mass ratio µ = m1m+m 2 2 , the coordinates of the primaries m1 , m2 become (µ, 0, 0), (µ − 1, 0, 0), their masses become 1−µ, µ respectively, and their period of revolution becomes 2π. The motion of the infinitesimal mass is then described by the autonomous Hamiltonian system with Hamiltonian 1 1−µ µ H(x1 , x2 , x3 , p1 , p2 , p3 ) = (p21 + p22 + p23 ) − x1 p2 + x2 p1 − − , 2 r1 r2 where r12 = (x1 − µ)2 + x22 + x23 , r22 = (x1 − µ + 1)2 + x22 + x23 . The value of the hamiltonian will be denoted as “the energy” from now on. The RTBP is shown to have 5 fixed points: the collinear ones, L1 , L2 , L3 , due to Euler, and the triangular ones, L4 , L5 , due to La- grange (see e.g. [51]). Following the astrodynamical convention, we will consider L1 to be the point located between the primaries. The x coordinate of this point is xL1 = µ − 1 + γ1 , with γ1 the positive root of one of Euler’s quintic equations, γ15 − (3 − µ)γ14 + (3 − 2µ)γ13 − µγ12 + 2µγ1 − µ = 0. The linear behaviour around L1 is of the type center×center×saddle. Namely, for the value of µ we use, Spec DXH (L1 ) = {i2πωp0 , −i2πωp0 , i2πωv0 , −i2πωv0 , λ0 , −λ0 }. From now on we will focus our attention on this point, and we will consider the primaries to be the Earth and the Moon, with mass pa- rameter µ = 1.215058560962404 × 10−2 , for which xL1 ≈ −0.83692, ωp0 ≈ 0.371529, ωv0 ≈ 0.361096, λ0 ≈ 2.932056. Lyapunov’s center theorem (see e.g. [43, 48]) ensures the existence of a family of periodic orbits (p.o.), known as the planar (resp. verti- cal) Lyapunov family, that fills a 2D manifold tangent to the ±i2πωp0 (resp. ±i2πωv0 ) eigenplane. The planar (resp. vertical) denomination comes from the fact that the eigenvectors of eigenvalues ±i2πωp0 (resp. ±i2πωv0 ) have zero x3 , p3 (resp. x1 , x2 , p1 , p2 ) coordinates. Both fami- lies start at the energy of L1 , that will be denoted as h0 , and evolve through higher energies. Denote by Tph (resp. Tvh ) the period of the pla- h nar (resp. vertical) Lyapunov p.o. of energy h. Denote also as e±i2πνp 32 ALEX HARO AND J.M. MONDELO h (resp. e±i2πνv ) the multipliers of modulus one of the monodromy ma- trix of the planar (resp. vertical) Lyapunov p.o. of energy h, with νph (resp. νvh ) chosen in [0, 1/2], as found when computing numerically. This is, h e±i2πνj ∈ Spec DϕTjh (xhj ), νjh ∈ [0, 1/2], j = p, v, where xhp (resp. xhv ) is an initial condition in the planar (resp. vertical) periodic orbit of energy h. Lyapunov’s center theorem also ensures that h→h h→h Tph −→0 1/ωp0 , Tvh −→0 1/ωv0 , and h h→h 0 0 h h→h 0 0 {e±i2πνp } −→0 {e±i2πωv /ωp }, {e±i2πνv } −→0 {e±i2πωp /ωv }. From the numerical values of ωp0 , ωv0 , we have h→h h→h νph −→0 1 − ωv0 /ωp0 , νvh −→0 ωp0 /ωv0 − 1. The Lissajous family of tori mentioned above is made of the KAM tori generated by the 4D central part of L1 , that are contained inside the 4D center manifold of this point. Denote by ωp , ωv the frequencies of any torus in the family, chosen as to have ωp → ωp0 and ωv → ωv0 as the torus collapses to L1 . We will refer to ω̂ = (ωp , ωv ) as the vector of natural frequencies. Denote by K̂ a parameterization of a torus K̂ of the family satisfying the invariance equation (5). Define (70) νp (ωp , ωv ) = 1 − ωv /ωp , νv (ωp , ωv ) = ωp /ωv − 1. As stated previously, we will not compute a parameterization K̂ of the whole torus, but of an invariant curve parameterized by K inside it (recall that we actually compute a collection of invariant curves). Following [27], we will use the energy h and the vertical rotation number ρ = νv as parameters in order to represent the tori of the Lissajous family. It was numerically found that, when varying h, νv in the region enclosed by the α, β, γ curves of Fig. 1, they uniquely determine a torus in the family. The α curve, with coordinates (h, νvh ), represents the vertical Lyapunov family, from its birth at L1 at energy h = h0 (point A) to its first 1:1 bifurcation, at energy h = −1.49590 (point D). The β curve represents the planar Lyapunov family from its birth to its first 1:1 bifurcation, at energy h = hB = −1.58718 (point B, in which the Halo family of p.o. appears). In order to have continuity at the point A, the vertical coordinate of the points of the β curve is not νph but the νv (ωp , ωv ) value of limiting nearby tori. From (70), it is found to be 1 − 1. 1 − νph The γ curve, which is the segment from point B to point D, corresponds to the separatrix between the Lissajous family of tori and other families COMPUTATION OF INVARIANT TORI 33 of quasi-periodic motion in the center manifold of L1 , that are described in [27]. As it has been mentioned in the introduction, the large matrix1 ap- proach in [27] was to write K(θ) as a truncated Fourier series, K(θ) = PNf  A0 + k=1 Ak cos(2πkθ) + Bk sin(2πkθ) , an then turn Eq. (7) into a finite non-linear system of equations by imposing it at 1 + 2Nf equally spaced values of θ. Multiple shooting was also implemented: all the computations were done with m = 2. The computational bottleneck of this procedure is that large values of Nf give rise to large systems of equations. The different sub-regions inside the α, β, γ curves in Fig. 1, that are not disjoint but nested, are labeled according to the value of Nf obtained in the computations of [27] for the tori inside them. A global upper limit of 100 was chosen for Nf , so tori in the ”> 100” sub-region were actually not computed. 0.09 α 0.08 <25 0.07 0.06 <50 0.05 ρ 0.04 <100 0.03 A 0.02 β >100 0.01 0 γ B C D -1.6 -1.59 -1.58 -1.57 -1.56 -1.55 -1.54 -1.53 -1.52 -1.51 -1.5 -1.49 h Figure 1. (Adapted from [27]) Energy-rotation number representation of the Lissajous family of invariant tori of the RTBP around L1 for the Earth-Moon mass parame- ter. 4.2. On the generators of tori in the Lissajous family. As stated previously, we will not compute a parameterization K̂ of the whole torus, but of an invariant curve parameterized by K inside it (recall that we actually compute a collection of invariant curves), to which we will refer to as a generator of the torus. The choice has consequences in the determination of the geometrical observables and the selection of the frequencies. We consider two cases, that we will distinguish as vertical generators and planar generators. We will denote as vertical generator of the invariant torus a param- eterized curve of the form Kv (θ) = K̂(θ, θv∗ ), for some fixed θv∗ ∈ [0, 1]. 1According to the nomenclature of [31]. 34 ALEX HARO AND J.M. MONDELO A calculation shows that  ϕ1/ωv Kv (θ)) = Kv θ + νv , which is Eq. (7) with T = 1/ωv , ω = νv , with νv defined as in (70). Close to vertical Lyapunov p.o., an invariant curve of the linearized flow around the Lyapunov p.o. satisfies the linearized version of the previous equation and thus provides an approximate solution, in order to obtain a first torus and start continuation. When globalizing the invariant curve Kv (for the time-T flow) to the invariant torus (for the vector field) via Eq. (8) we get K̂v (θ̂) = K̂v (θ1 , θ2 ) = ϕθ2 /ωv (Kv (θ1 − θ2 νv )) = K̂(Âv θ̂ + θ̂v∗ ), where     1 1 ∗ 0 (71) Âv = , θ̂v = ∗ . 0 1 θv This is a reparameterization of the invariant torus K, for which the frequencies are ω̂v = Â−1 v ω̂ = (ωp − ωv , ωv ). The geometrical observables provided by the Calabi invariants of the two parameterizations of the torus K are related by the identities     1 0 1 0 C(K̂v ) = C(K̂), C(K̂) = C(K̂v ), 1 1 −1 1 see Remark 2.3.1, and, hence, C1 (K̂) = C1 (K̂v ) = C(Kv ) (which also follows from the definition of Kv ) and C2 (K̂) = −C1 (K̂v ) + C2 (K̂v ). Notice that C2 (K̂v ) = C(ϕθ/ωv (Kv (−θνv ))) = C(ϕT θ (Kv (−θω))). In summary, we can compute the Calabi invariants of K̂ from K̂v . We will denote as planar generator of the invariant torus a parame- terized curve of the form Kp (θ) = K̂(θp∗ , −θ), for some fixed θp∗ ∈ [0, 1]. A calculation shows then that  ϕ1/ωp Kp (θ)) = Kp θ + νp , wich is Eq. (7) with T = 1/ωp , ω = νp , with νp defined as in (70). Close to a planar Lyapunov orbit, an invariant curve of the linearized flow around the Lyapunov p.o. satisfies the linearized version of the previous equation, and thus provides an approximate solution, in order to obtain a first torus and start continuation. When globalzing the invariant curve Kp (for the time-T flow) to the invariant torus (for the vector field) via Eq. (8) we get in this case K̂p (θ̂) = K̂p (θ1 , θ2 ) = ϕθ2 /ωp (Kp (θ1 − θ2 νp )) = K̂(Âp θ̂ + θ̂p∗ ), where    ∗ 0 1 ∗ θ (72) Âp = , θ̂p = p . −1 1 0 COMPUTATION OF INVARIANT TORI 35 This is another reparameterization of the invariant torus K, for which the frequencies are ω̂p = Âp−1 ω̂ = (ωp − ωv , ωp ). The Calabi invariants of the two parameterizations of the torus K are related by the identities     0 −1 1 1 C(K̂p ) = C(K̂), C(K̂) = C(K̂p ), 1 1 −1 0 see Remark 2.3.1, and, hence, C1 (K̂) = C1 (K̂p ) + C2 (K̂p ) and C2 (K̂) = −C1 (K̂p ) = −C(Kp ) (as it follows from the definition of Kp ). Notice that C2 (K̂p ) = C(ϕθ/ωp (Kp (−θνp ))) = C(ϕT θ (Kp (−θω))). In sum- mary, we can compute the Calabi invariants of K̂ from K̂p . 4.3. The numerical explorations. We have performed two numeri- cal explorations of the Lissajous family. In the first one we compute tori with ρ > νvh0 , whereas in the second one we compute tori with ρ < νvh0 . Since the tori with ρ > νvh0 are included in the “< 100” sub-region of Fig. 1, in the first exploration we are able to compare the performance of this paper’s parameterization procedure against the one of the large matrix approach of [27]. In our first exploration, we have chosen 74 values of ρ > νvh0 , equally spaced between 0.02944 and 0.08754, and “nobilized”2 with an abso- lute tolerance of 1.6 × 10−4 . For each of these values of ρ, we have performed continuation of invariant curves given by vertical generators Kv , for constant ω = ρ and increasing T , starting from a curve on a narrow torus around a vertical Lyapunov p.o. and finishing by collaps- ing to another vertical Lyapunov p.o. of a higher energy. We have also simultaneously computed their invariant bundles. To do so, we have performed continuation with respect to T using Algorithm 3.6.1, mak- ing predictions through Algorithm 3.5.1, and refining each prediction through Algorithm 3.4.1 (isochronous case). The parameters used in Algorithm 3.6.1 have been: m = 4, ε = 10−7 , εW = 10−5 , ε1 = 10−8 , ε2 = 10−12 , ndes = 4, nα = 5. The stopping criterion has been that, when approaching the final vertical p.o., |C1 (K̂)| = |C(Kv )| < 0.001 (see Section 4.2). We have repeated this first exploration using the large-matrix approach of [27], selecting the parameters accordingly for a fair comparison. The results of this numerical exploration are shown in Fig. 2. The left plot corresponds to the large matrix approach, whereas the right plot corresponds to this paper’s parameterization one. Both plots show the total number of Fourier coefficients used in the computation. For the left plot this is 1 + 2Nf . For the right plot, it is considered to be N/2, because of the DFT queue cleaning strategy mentioned in 2A noble number is one whose continued fraction expansion coefficients are equal to one from a position on. 36 ALEX HARO AND J.M. MONDELO Section 3.6. The total computing time3 of the large-matrix approach is 68308 seconds, of which 22793 are spent in the computation of the tori, whereas the rest are used in the computation of the stable and unstable bundles through a slight modification of the method presented in [40]. The total computing time of the parameterization approach of this paper (that includes tori and bundles) is 5992 seconds. The total number of tori computed is 4141 with the large-matrix method vs. 7008 with the parameterization one. This is due to the fact that, as can be appreciated in Fig. 2, the continuation strategy of the large-matrix procedure is able to use larger step sizes for tori with small values of Nf . 1+2Nf N/2 0.09 300 0.09 256 0.08 250 0.08 0.07 200 0.07 128 ρ 0.06 150 ρ 0.06 0.05 0.05 64 100 0.04 0.04 50 32 0.03 0.03 0 -1.59 -1.56 -1.53 -1.5 -1.59 -1.56 -1.53 -1.5 h h Figure 2. Tori with ρ > νvh0 of the Lissajous family of the Earth-Moon RTBP around L1 , computed with the large matrix procedure (left) and the parameterization one (right). In a second exploration, we have computed invariant tori that are born from planar Lyapunov orbits that have rotation number ρ < νvh0 . Namely, we have chosen 31 values of ρ, equally spaced between 0.02785 and 0.00317, nobilized with an absolute tolerance of 1.6 × 10−4 and for each of the values we have performed continuation of invariant curves given by planar generators Kp , for constant ω = 1 − 1/(1 + ρ) and increasing T , starting from a narrow torus around a planar Lyapunov p.o. The parameters used in the continuation algorithms have been the same as before. The continuations have stopped either by reaching the computational limit or, as in the first exploration, when |C1 (K̂)| = |C(Kv )| < 0.001 (see Section 4.2). The motivation of this second exploration is twofold. On the one hand, to perform a “stress test” of our procedure by exploring phase space beyond the computations of [27]. On the other, to relate the 3Thecomputing times given will always be qualified as “total”, meaning the sum of all the times used by all the threads. The actual wall-clock time is roughly this time divided by the number of cores (16 in our case). This rule is not followed exactly because of uneven load balancing: the continuation of some constant-ω families takes longer than others. COMPUTATION OF INVARIANT TORI 37 N/2 0.09 8192 0.08 4096 0.07 2048 0.06 1024 0.05 ρ 512 0.04 256 0.03 128 0.02 64 0.01 32 0 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 h Figure 3. Number of Fourier coefficients of all the tori and invariant bundles computed of the Lissajous family of the Earth-Moon RTBP around L1 , using this paper’s parameterization strategy. behavior of dynamical and geometrical observables to the destruction of invariant tori (see the next section). In this exploration, we have always achieved convergence of Newton’s method provided that the continuation step is small enough and the number of points N large enough. The computational limits on these two quantities, chosen as 10−5 and 8192, respectively, have been set in order to obtain reasonable run time and storage requirements in a single workstation. Notice that working with such a large number of Fourier coefficients is unfeasible with the large matrix approach, since it would require the solution of non-linear sytems of equations with a size of nearly 100000 × 100000. A limit, equal to 10000, has also been put for the maximum number of tori computed in each constant ρ family of this second exploration. With all these limits, this exploration has run for a total time of 27.0816 days and has generated a total of 130574 tori that, each compressed as a bz2 file, take up 127.51 GiB of disk space.4 The results of this second exploration are shown in Fig. 3, that is analogous to Fig. 2 right but including both explorations. Many of the figures that follow will refer to the two explorations as a whole. 4.4. Dynamical and geometric observables. In this section we will describe the behaviour of different dynamical and geometric observ- ables of the invariant tori of the Lissajous family that are obtained 4Recall that the “wall-clock” time is roughly the total time divided by 16. On the other hand, several strategies, that we have not pursued here, can be used to reduce greatly these storage requirements, like using binary files with single precision floating-point numbers, and not storing all the tori but a grid of them fine enough in order to recover tori not in the grid by interpolation (see e.g. [44]). 38 ALEX HARO AND J.M. MONDELO during their computation. These observables provide insight on both numerical and dynamical vicissitudes faced by the method. We will also display the evolution of three constant ρ families through the 3D representation of some of their tori in configuration space. The results presented here complement those in [27] that, by using iso-energetic Poincaré sections, provides a detailed account of the evolution of the Lissajous family of invariant tori, together with its interaction with other families of invariant tori and periodic orbits. For a parameterization K = K0 : T → R6 of the invariant curve of the time T -flow with rotation number ω, obtained with multiple shooting with m = 4 steps, that is a generator of the 2D invariant torus given by (8), the dynamical observables we consider here are: • the flying time T of the generator (notice that we get frequencies for the 2D-torus by the formula ω̂ = T1 (ω, 1)); • the unstable Floquet multiplier Λu = λ−m , that provides in- formation about hyperbolity properties of the invariant curve, and from which one can obtain the Floquet exponent of the 2D-torus by χ = T1 log Λu . Notice we select m so that λ−1 is not too big, in order to mitigate numerical unstabilities. In our computations, λ−1 runs in the interval [4, 8]. The geometric observables we consider are: • The Calabi invariants C1 (K̂), C2 (K̂) of the parameterization K̂(θ1 , θ2 ) of the torus with natural frequencies ω̂ = (ωp , ωv ). These invariants give insight on the size of the generators in area units. Section 4.2 provides formulae for their computation from C(K̂p ), C(K̂v ). • The (minimum) distances between several pairings of bundles on the generator curve: T K, the tangent bundle of the generator curve (generated by K 0 ); X, the bundle generated by the vector field on the curve; the stable and the unstable bundles E s and E u , respectively; the central bundle, E c , that has rank 4 and contains the tangent bundle to the generator, the vector field and, hence, the tangent bundle to the 2D torus. The distances we consider are: d(T K, X), to measure the transversality of the flow to the generator, and d(E s , E u ), d(E s , E c ), and d(E u , E c ), to measure the quality of hyperbolicity geometrical properties. The bundles are generated by selected columns of the matrix map P = P0 : T → R6×6 , that we write as  P (θ) = K 0 (θ) XH (K(θ)) W s (θ) N 1 (θ) N 2 (θ) W u (θ) for reference. Hence, at a point K(θ) of the invariant curve, the fiber of the stable bundle E s is generated by W s (θ), the fiber of the unsta- ble bundle E u is generated by W u (θ), and fiber of the center bundle COMPUTATION OF INVARIANT TORI 39 E c is generated by K 0 (θ), XH (K(θ)), N 1 (θ), N 2 (θ). Notice that the tangent bundle to the generator is generated by K 0 (θ), and the tangent bundle of the 2D torus is generated by K 0 (θ), XH (K(θ)). There are several ways of defining distances or angles between vector subspaces of a given normed vector space. Here, the vector space is R6 , with the norm induced by the standard scalar product, and the distance we con- sider between a vector subspace E1 of dimension 1 and another vector subspace E2 is the length of the projection onto E2⊥ of a unit vector in E1 (the angle between E1 and E2 is the arcsin of this lenght). Finally, we define the distance between two bundles as the minimum distance between corresponding fibres of the bundles. Remark 4.4.1. We emphasize that the quality of hyperbolicity proper- ties of the invariant torus are not only given by the Floquet multipliers in the stable and unstable directions, that have to be away from 1, but also by the positivity of the angles between the stable, unstable and center directions. There are mechanisms of breakdown of invariant tori that involve the degeneration of some of these angles, that go to zero, while the stable and unstable Floquet multipliers remain far from 1. See [10, 11, 12, 23, 30, 33]. Remark 4.4.2. Reversibility properties of the RTBP imply that stable and unstable bundles can be obtained from each other using reversors, and that they have same angles with the center manifold of the torus. These properties could also be used to reduce the cost of the algorithms presented here (reducing, for instace, the cost of generating the frame). We prefer not doing so for the sake of generality. By monitoring these observables during the continuation we can get insight about dynamical and geometric properties of the torus (and its invariant bundles), and detect numerical unstabilities caused by degen- eracies of these properties (such as the hyperbolicity, regularity of the frame, size of the generator). For instance, we recall that one stopping criterion is that |C1 (K̂)| = |C(Kv )| < 0.001, revealing that the torus is approaching a periodic orbit. We have collected these observables from the two numerical explorations exposed in the previous section, and the results are summarized in Figure 4. In this massive computation we observe that: • The unstable multiplier ranges from 413.205 to 3344.26, from which the spectral condition of hyperbolicity of the torus is satisfied; • The distance between the stable and unstable bundles is bigger than 0.492489, and the distante between the stable and cen- ter bundle, and the unstable and center bundle, is bigger than 0.0615721, from which the geometrical conditions of hyperbol- icy is also satisfied. 40 ALEX HARO AND J.M. MONDELO The continuations of families of Lissajous tori with smaller rotation numbers stop because of the computational limit of 10000 tori for each constant-ρ family. But, as we see from the behavior of the observables, this phenomenon is not apparently due to the fact that hyperbolicity breaks down. However, the fact that the step size of the continuation becomes smaller and the number of Fourier coefficients becomes larger reveals that the torus is losing regularity (the analyticy strip of the complex domain of the parameterization of the torus goes to zero), indicating an obstruction for the existence of the torus and that it is breaking down. There is another possible mechanism of breakdown, it is what we call KAM breakdown. These tori lie on the center manifold of the L1 point, W c (L1 ), which is a 4D symplectic manifold. So, inside W c (L1 ), these tori are KAM tori, and the basic mechanism of break- down is the collision with resonances (the overlap criterion [16]), which can be more geometrically described as the obstruction produced by homoclinic and heteroclinic webs produced by the invariant manifolds of unstable periodic orbits inside the center manifold (the obstruction criterion in [45, 19]). So, it is very likely that in this case the break- down is produced by this phenomenon inside the center manifold. We will come back to this issue later. In the following, we will particularize the results for three families of Lissajous tori, with rotation numbers ρ = 0.031865, ρ = 0.019091 and ρ = 0.013584. Figs. 5, 6 and 7 show several samples of tori of these three families, projected on the configuration space in different forms and views: (left) as grids on the parameterized surfaces {K̂(θ1 , θ2 )}(θ1 ,θ2 )∈T2 (see Eq. (8)), including two generators of the homothopy group of the torus, given by {K̂(θ, 0)}θ∈T , in blue, and {K̂(0, θ)}θ∈T , in red; (right) as opaque surfaces, with the same scale and range on all axes, with colors corresponding to the different sides of the surface, revealing self-intersections of the projections of tori on configu- ration space. Actually, instead of Eq. (8), the expression  j  K̂(θ1 , θ2 ) = ϕ(θ2 − j )T Kj θ1 − (θ2 − )ω , m m with j = [mθ2 ], has been used, in order to take advantage of multiple shooting. In addition to these views, for each of these families we have plotted in Figs. 9, 10 and 11 the dynamical and geometrical observables as functions of the energy h. We will describe our findings below. The results for the family ρ = 0.031865 are summarized in Figs. 5 and 9. From Fig. 5 we appreciate how the family begins with a small torus around a vertical Lyapunov p.o., that grows up to approximately the size of the planar Lyapunov orbit of the same energy, and then COMPUTATION OF INVARIANT TORI 41 T Λu 0.09 4.2 0.09 0.08 4 0.08 3000 0.07 3.8 0.07 0.06 0.06 2500 3.6 0.05 0.05 2000 ρ 3.4 ρ 0.04 0.04 3.2 0.03 1500 0.03 0.02 3 0.02 1000 0.01 2.8 0.01 0 500 0 2.6 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 h h ^ )| ^ )| |C1(K |C2(K 0.09 1E+00 0.09 1E+00 0.08 0.08 0.07 0.07 0.06 1E-01 0.06 1E-01 0.05 0.05 ρ ρ 0.04 0.04 0.03 1E-02 0.03 1E-02 0.02 0.02 0.01 0.01 0 1E-03 0 1E-03 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 h h s u d(TK,X) d(E ,E ) 0.09 1E+00 0.09 1 0.08 0.08 0.07 0.07 0.8 1E-01 0.06 0.06 0.6 0.05 1E-02 0.05 ρ ρ 0.04 0.04 0.4 0.03 0.03 1E-03 0.02 0.02 0.2 0.01 0.01 0 1E-04 0 0 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 h h s c u c d(E ,E ) d(E ,E ) 0.09 1 0.09 1 0.08 0.08 0.07 0.8 0.07 0.8 0.06 0.06 0.6 0.6 0.05 0.05 ρ ρ 0.04 0.04 0.4 0.4 0.03 0.03 0.02 0.2 0.02 0.2 0.01 0.01 0 0 0 0 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 h h Figure 4. Dynamical and geometrical observables starts bending until it “is about to close”. Then it opens again and shrinks until it collapses to a vertical Lyapunov p.o. of higher energy. All this is done while increasing in size, since energy also increases. Fig. 9 displays the observables for this family. Since the family is born in a vertical Lyapunov p.o. and dies in another vertical Lyapunov p.o. of a higher energy, the Calabi invariant of the generator starts being 0 and finishes being 0. Notice that, in both cases, close to the p.o. the Calabi invariant goes to zero assymptotically as a linear function of the difference of the energy with the one of the p.o. 42 ALEX HARO AND J.M. MONDELO 0.2 0.04 0.1 z 0 0.02 -0.1 z 0 -0.2 -0.02 0.0008 0.2 -0.04 0.0004 0.1 0 0 -0.8385 y -0.1 -0.8382 -0.0004 -0.2 -0.9 -0.8 y -0.8379 -0.8376 -0.0008 x x -0.8373 0.2 0.06 0.1 0.04 z 0 0.02 -0.1 z 0 -0.2 -0.02 -0.04 0.06 0.2 -0.06 0.03 0.1 0 0 -0.86 y -0.1 -0.85 -0.03 -0.2 -0.9 -0.8 y -0.84 -0.83 -0.06 x x -0.82 0.2 0.12 0.1 0.09 z 0 0.06 0.03 -0.1 z 0 -0.03 -0.2 -0.06 -0.09 0.12 0.2 -0.12 0.06 0.1 0 0 -0.88 y -0.1 -0.86 -0.06 -0.2 -0.9 -0.8 y -0.84 -0.82 -0.12 x x -0.8 0.2 0.2 0.1 0.15 z 0 0.1 0.05 -0.1 z 0 -0.05 -0.2 -0.1 -0.15 0.21 0.2 0.14 0.1 -0.2 0.07 0 0 -0.93 -0.07 y -0.1 -0.9 -0.2 -0.9 -0.8 -0.87 -0.14 y -0.84 -0.21 -0.81 x x -0.78 0.2 0.25 0.1 0.2 z 0 0.15 0.1 0.05 -0.1 z 0 -0.05 -0.2 -0.1 -0.15 0.27 0.2 -0.2 0.18 0.1 -0.25 0.09 0 0 -0.96 -0.93 -0.09 y -0.1 -0.9 -0.87 -0.2 -0.9 -0.8 -0.18 y -0.84 -0.81 -0.78 -0.27 x x 0.2 0.25 0.1 0.2 z 0 0.15 0.1 0.05 -0.1 z 0 -0.05 -0.2 -0.1 -0.15 0.06 0.2 -0.2 0.1 -0.25 0.03 0 0 -0.94 y -0.1 -0.92 -0.03 -0.2 -0.9 -0.8 -0.9 y -0.88 -0.06 -0.86 x x -0.84 Figure 5. Two views of a sample of invariant tori with ρ = 0.031865. Letf: Parameterized surfaces, including generators (Blue: invariant curve {K̂(θ, 0)}θ∈T . Red: in- variant curve {K̂(0, θ)}θ∈T ). Right: with the same scale on all axes, and with opaque surfaces. COMPUTATION OF INVARIANT TORI 43 The results for the family ρ = 0.019091 are summarized in Figs. 6 and 10. Fig. 6 shows a sample tori starting in a planar Lyapunov p.o. and, hence, their homothopy group generators are exchanged with respect to Fig. 5. The evolution with energy is similar to the one of Fig. 6, with three main differences: the torus “seems to close after bending” (a zoom of the fourth torus reveals that it does not actually close), there is a more important accumulation of wireframe lines at the “boundary that closes and opens”, and the family starts from a planar Lyapunov p.o. instead of a vertical one. This is in fact the reason the Calabi invariant of the computed generators starts being 0 (and, again, with an asymptotic linear behaviour) and increases till the end of the continuation (converging to the Calabi invariant of the vertical Lyapunov p.o.), as it is observed in Fig. 10. In the iso-energetic Poincaré section plots of [27] (and, for lower en- ergy levels, in the ones of previous references like [26, 41]), it is numer- ically seen that, from energy hB to hC , a double homoclinic connection inside W c (L1 ) of the planar Lyapunov family of p.o. acts as separa- trix from the Lissajous family and the quasi-Halo family of tori. From energy hC to hD , this role is taken by heteroclinic connections (also in- side W c (L1 )) between the vertically symmetric families of p.o. that are born at the second 1:1 bifurcation of the planar Lyapunov p.o. Plots of orbits of these last families of p.o. can be found in Fig. 8 of [27] (they are known as “axial” by other autors, e.g. [20]). In Fig. 6, it is observed how the tori approach these connections. This fact is better appreciated in Fig. 8, that shows magnifications of projections of the third and fifth tori of Fig. 6. The view of the third torus has been chosen in order to stress the fact that the torus represented approaches two different, vertically simmetric quasi-Halo tori, as can be inferred from the Poincaré representations of the center manifold in references [27, 26, 41]. When approaching these connections, dynamics becomes slow and, since the parameterization K̂ is tied to the dynamics through the invariance equation (5), it produces the accumulation of wireframe lines of Fig. 6 and the small values of d(T K, X) of Fig. 4. We believe this kind of stiffness to be responsible for the drastic reduction of step length of the second exploration and, to a lesser extent, of the first one. As it has been commented, these connections are responsible for the destruction of the families of invariant tori in the center manifold. The phenomenon of breakdown is illustrated with the family with ρ = 0.013584. The results for this family are summarized in Figs. 7 and 11. In Fig. 7 we observe that the torus is increasingly pinched, while dynamical and geometrical observables in Fig. 11 do not sug- gest the torus is being destroyed. However, the fact the continuation step is getting very small and the number of Fourier coefficients of the approximations is getting larger, related to the mentioned “pinching 44 ALEX HARO AND J.M. MONDELO 0.2 0.0008 0.1 0.0006 z 0 0.0004 0.0002 -0.1 z 0 -0.0002 -0.2 -0.0004 -0.0006 0.06 0.2 -0.0008 0.04 0.1 0.02 0 -0.85 -0.845 0 y -0.1 -0.02 -0.9 -0.8 -0.84 -0.835 -0.04 y -0.2 -0.83 -0.825 -0.06 x x 0.2 0.06 0.1 0.04 z 0 0.02 -0.1 z 0 -0.02 -0.2 -0.04 0.08 0.2 -0.06 0.04 0.1 0 0 -0.87 -0.86 y -0.1 -0.04 -0.9 -0.8 -0.85 -0.84 y -0.2 -0.83 -0.08 x -0.82 x 0.2 0.08 0.1 0.06 z 0 0.04 0.02 -0.1 z 0 -0.02 -0.2 -0.04 -0.06 0.1 0.2 -0.08 0.05 0.1 0 0 -0.87 -0.86 y -0.1 -0.85 -0.84 -0.05 -0.9 -0.8 y -0.2 -0.83 -0.82 -0.1 x -0.81 x 0.2 0.15 0.1 0.1 z 0 0.05 -0.1 z 0 -0.05 -0.2 -0.1 0.2 0.2 -0.15 0.1 0.1 0 0 -0.92 -0.9 y -0.1 -0.88 -0.86 -0.1 -0.9 -0.8 y -0.2 -0.84 -0.82 -0.2 x -0.8 x 0.2 0.25 0.1 0.2 0.15 z 0 0.1 0.05 -0.1 z 0 -0.05 -0.2 -0.1 -0.15 0.3 0.2 -0.2 0.2 -0.25 0.1 0.1 0 0 -1 -0.96 -0.1 y -0.1 -0.92 -0.88 -0.9 -0.8 -0.2 y -0.2 -0.84 -0.8 -0.3 x -0.76 x 0.2 0.25 0.1 0.2 0.15 z 0 0.1 0.05 -0.1 z 0 -0.05 -0.2 -0.1 -0.15 0.2 0.2 -0.2 -0.25 0.1 0.1 0 0 -0.96 -0.93 y -0.1 -0.9 -0.87 -0.1 -0.9 -0.8 y -0.2 -0.84 -0.81 -0.2 x -0.78 x Figure 6. Two views of a sample of invariant tori with ρ = 0.019091. Left: parameterized surfaces, including generators (Blue: invariant curve {K̂(θ, 0)}θ∈T . Red: in- variant curve {K̂(0, θ)}θ∈T ). Right: with the same scale on all axes, and with opaque surfaces. COMPUTATION OF INVARIANT TORI 45 0.15 0.1 0.0009 0.05 0.0006 z 0 0.0003 -0.05 -0.1 z 0 -0.15 -0.0003 0.2 -0.0006 0.06 0.1 0.03 0 -0.0009 y -0.1 -0.86 0 -0.85 -0.03 y -0.2-0.9 -0.8 -0.84 x -0.83 -0.06 x -0.82 0.15 0.1 0.06 0.05 0.04 z 0 0.02 -0.05 -0.1 z 0 -0.15 -0.02 0.2 -0.04 0.08 0.1 0.04 0 -0.06 y -0.1 -0.87 -0.86 0 -0.04 y -0.2-0.9 -0.8 -0.85 -0.84 x -0.83 -0.08 x -0.82 0.15 0.1 0.08 0.05 0.06 z 0 0.04 -0.05 0.02 -0.1 z 0 -0.15 -0.02 0.2 -0.04 0.1 -0.06 0.1 0.05 0 -0.08 y -0.1 -0.87 -0.86 0 -0.05 -0.2-0.9 -0.8 -0.85 -0.84 y -0.83 -0.82 -0.1 x x -0.81 0.15 0.1 0.15 0.05 0.1 z 0 0.05 -0.05 -0.1 z 0 -0.15 -0.05 0.2 -0.1 0.2 0.1 0.1 0 -0.15 y -0.1 -0.92 -0.9 0 -0.1 -0.2-0.9 -0.8 -0.88 -0.86 y -0.84 -0.82 -0.2 x x -0.8 Figure 7. Two views of a sample of invariant tori with ρ = 0.013584. Left: parameterized surfaces, including generators (Blue: invariant curve {K̂(θ, 0)}θ∈T . Red: in- variant curve {K̂(0, θ)}θ∈T ). Right: with the same scale on all axes, and with opaque surfaces. phenomenon”, envisages the breakdown of the torus inside the center manifold. 5. Conclusion We have presented in this paper a very efficient method to com- pute invariant tori in Hamiltonian systems. To do so, we have first reduced the dimensionaly of the objects, by considering invariant tori for flow maps, and then taken advantage of the geometrical and dynam- ical properties of invariant tori in Hamiltonian systems. The method also provides online information on the linearized dynamics around the 46 ALEX HARO AND J.M. MONDELO 0.08 0.06 0.04 0.02 z 0 -0.02 -0.04 -0.06 -0.08 -0.87 -0.86 -0.85 -0.84 -0.83 -0.82 -0.81 x 0.25 0.2 0.15 0.1 0.05 z 0 -0.05 -0.1 -0.15 0.25 -0.2 0.2 0.15 -0.25 0.1 0.05 -0.98 0 -0.96 -0.94 -0.05 y -0.92 -0.9 -0.1 -0.88 -0.15 -0.86 x -0.84 -0.2 -0.82 -0.8 -0.78 -0.25 Figure 8. Magnifications of projections of the third and fifth tori of Fig. 6, corresponding to ρ = 0.019091. torus, as well as on other geometrical properties. He have focused our attention on partially hyperbolic invariant tori with rank-one stable and unstable manifolds, and in this case the method provides not only parameterizations of the tori but also of the linear approximations of those manifolds. Tests have been performed for the computation of invariant tori around libration points of the Circular Restricted Three Body Problem, for which there is an extensive literature, so any one could easily compare the performances. COMPUTATION OF INVARIANT TORI 47 0.1 1 0.09 0.08 0.8 0.07 0.06 0.6 C 0.05 d 0.04 0.4 d(Es,Eu) 0.03 ^) d(Eus,TK 0.02 0.2 ^) d(E s,TK 0.01 d(Eu,Ecc) d(E ,E ) 0 0 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 h h Figure 9. Calabi invariant C(Kv ) = C1 (K̂), and dis- tances between several pairings of bundles, for the family ρ = 0.031865. 0 1 -0.1 0.8 -0.2 -0.3 0.6 C d -0.4 0.4 d(Es,Eu) -0.5 s ^ d(Eu,TK ) 0.2 ^) d(E s,TK -0.6 d(Eu,Ecc) d(E ,E ) -0.7 0 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 h h Figure 10. Calabi invariant C(Kp ) = −C2 (K̂), and distances between several pairings of bundles, for the family ρ = 0.019091. 0 1 -0.1 0.8 -0.2 -0.3 0.6 C d -0.4 0.4 s u d(E ,E ) -0.5 d(Es,TK^) 0.2 ^) d(Eus,TK -0.6 d(Eu,Ecc) d(E ,E ) -0.7 0 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 -1.6 -1.58 -1.56 -1.54 -1.52 -1.5 h h Figure 11. Calabi invariant C(Kp ) = −C2 (K̂), and distances between several pairings of bundles, for the family ρ = 0.013584. As we have already mentioned, we only present algorithms of compu- tation, based on Newton’s method. Eventually, a proof of convergence could be completed using KAM techniques for obtaining results in a posteriori format. Although we have not attempted this, we do provide information about the magical cancellations appearing in the linearized 48 ALEX HARO AND J.M. MONDELO equations, that are key in the KAM proofs, and, with much more effort, could be implemented as computer assisted proofs [22]. The increasing complexity of problems and applications has spurred the research of this paper. The algorithms presented here are the first ones in a new generation of algorithms to compute invariant tori and their manifolds in Hamiltonian systems. We plan to extend the methodologies to more complex problems in the future, with an eye in the applications. Appendix A. From Poincaré map to time-T map In this section we will see how to obtain, from an invariant torus of a the Poincaré map, an invariant torus of a time-T map. For instance, the invariant torus of the Poincaré map could have been computed from a center manifold reduction around an equilibrium point at a certain fixed energy level (see e.g. [49, 38] for normal form methods and [29] for direct parameterization methods, applied to the computation of the center manifold of a colinear fixed point in the RTBP). Assume we are given a parameterization KP : Td−1 → R2n of a (d − 1)-dimensional torus KP inside a d-dimensional torus K̂, produced by a Poincaré map P associated to a transversal section to the vector field XH . We also assume that the rotation vector of KP is ω (which is assumed to be Diophantine). That is, we assume for all θ ∈ Td−1 ϕTP (θ) (KP (θ)) = KP (θ + ω), where TP : Td−1 → R gives for each θ ∈ Td−1 the time for a point KP (θ) to return to the transversal section. The flying time depends then on the point on the torus. We want to find a parameterization K : Td−1 → R2n of a (d − 1)- dimensional torus K for which the flying time is constant. To do so, we look for τ : Td−1 → R and T such that K(θ) = ϕτ (θ) (KP (θ)), ϕT (K(θ)) = K(θ + ω). Hence, since ϕT (K(θ)) = ϕT +τ (θ) (KP (θ)) and K(θ + ω) = ϕτ (θ+ω) (KP (θ + ω)) = ϕτ (θ+ω)+TP (θ) (KP (θ)), we impose that τ (θ) − τ (θ + ω) = TP (θ) − T. This is the well-known small divisors equation (discussed here in Sec- tion 3.2). We adjust then T to be the average of TP , T = hTP i, and solve for τ . Notice that τ is defined up to a constant, the average that we take as 0. This is natural since K is determined up to a time translation (see Remark 2.4.3). COMPUTATION OF INVARIANT TORI 49 Appendix B. Quadratically small averages In this section we will prove that, given a multiple torus {Ki (θ)}m−1 i=0 , approximately invariant with errors ω Ei (θ) = ϕT /m (Ki (θ)) − Ki+1 (θ + m ), then the averages of {ηi3 (θ)}m−1 i=0 given by  ω > ω   31  3 −DKi+1 (θ + m ) Ω(Ki+1 (θ + m ))Ei (θ) η (θ) ηi (θ) = ω ω =: i32 −XH (Ki+1 (θ + m ))Ω(Ki+1 (θ + m ))Ei (θ) ηi (θ) are quadratically small with respect to the errors (and their deriva- tives). This is a crucial step in our algorithms, and also in a eventual proof of their convergence using KAM methods. First, we start by proving that hη 32 i is quadratically small. In fact m−1 X m−1 X hηi32 (θ)i = hDH(Ki+1 (θ + ω m ))Ei (θ)i i=0 i=0 m−1 X ω = hH(ϕT /m (Ki (θ)) − H(Ki+1 (θ + m ))i i=0 m−1 X Z 1 − h (1 − s) D2 H(Ki+1 (θ + ω m ) + sEi (θ))[Ei (θ), Ei (θ)] dsi i=0 0 m−1 X Z 1 =− h (1 − s) D2 H(Ki+1 (θ + ω m ) + sEi (θ))[Ei (θ), Ei (θ)] dsi, i=0 0 which is quadratically small in the errors. Second, we prove that hη 31 i is quadratically small. We start using the exactness of the symplectic form: m−1 X m−1 X ω > ω > hηi31 (θ)i =− hDKi+1 (θ + m ) Da(Ki+1 (θ + m )) Ei (θ)i i=0 i=0 m−1 X ω > ω + hDKi+1 (θ + m ) Da(Ki+1 (θ + m ))Ei (θ)i i=0 m−1 X = hDEi (θ)> a(Ki+1 (θ + ω m )) + DKi+1 (θ + ω > m ) (∆1 ai (θ) − ∆2 ai (θ))i, i=0 where we use ω >  0 =hD a(Ki+1 (θ + m )) Ei (θ) i =hEi (θ)> D(a(Ki+1 (θ + m ω ))) + a(Ki+1 (θ + ω > m )) DEi (θ)i, 50 ALEX HARO AND J.M. MONDELO and the definitions ∆1 ai (θ) =a(ϕT /m (Ki (θ))) − a(Ki+1 (θ + mω )) Z 1 ω = Da(Ki+1 (θ + m ) + sEi (θ))Ei (θ) ds 0 and ∆2 ai (θ) =a(ϕT /m (Ki (θ))) − a(Ki+1 (θ + mω )) − Da(Ki+1 (θ + m ω ))Ei (θ) Z 1 = (1 − s) D2 a(Ki+1 (θ + m ω ) + sEi (θ))[Ei (θ), Ei (θ)] ds. 0 Hence, using that ω DϕT /m (Ki (θ))DKi (θ) − DKi+1 (θ + m ) = DEi (θ) and that m−1 X m−1 X ω > hDKi+1 (θ + m ) a(Ki+1 (θ + ω m ))i = hDKi (θ)> a(Ki (θ))i i=0 i=0 we have: m−1 X m−1 X  31 hηi (θ)i = hDKi (θ)> DϕT /m (Ki (θ))> a(ϕT /m (Ki (θ))) − a(Ki (θ)) i i=0 i=0 m−1 X − hDEi (θ)> ∆1 ai (θ) + DKi+1 (θ + ω > 2 m ) ∆ ai (θ)i i=0 m−1 X > = hD pT /m (Ki (θ)) i i=0 m−1 X − hDEi (θ)> ∆1 ai (θ) + DKi+1 (θ + ω > 2 m ) ∆ ai (θ)i i=0 m−1 X =− hDEi (θ)> ∆1 ai (θ) + DKi+1 (θ + ω > 2 m ) ∆ ai (θ)i, i=0 which is quadratically small. We have used the exactness of the Hamil- tonian flow, being pt is the primitive function of ϕt . References [1] R.A. Adomaitis, I.G. Kevrekidis, and R. de la Llave. A computer-assisted study of global dynamic transitions for a noninvertible system. Internat. J. Bifur. Chaos Appl. Sci. Engrg., 17(4):1305–1321, 2007. [2] Eugene L. Allgower and Kurt Georg. Numerical continuation methods, vol- ume 13 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1990. An introduction. [3] Nicola Baresi, Zubin P. Olikara, and Daniel J. Scheeres. Fully numerical meth- ods for continuing families of quasi-periodic invariant tori in astrodynamics. The Journal of the Astronautical Sciences, 2018. COMPUTATION OF INVARIANT TORI 51 [4] E. Oran Brigham. The Fast Fourier Transform and its applications. Prentice– Hall, 1988. [5] X. Cabré, E. Fontich, and R. de la Llave. The parameterization method for invariant manifolds. I. Manifolds associated to non-resonant subspaces. Indiana Univ. Math. J., 52(2):283–328, 2003. [6] X. Cabré, E. Fontich, and R. de la Llave. The parameterization method for invariant manifolds. II. Regularity with respect to parameters. Indiana Univ. Math. J., 52(2):329–360, 2003. [7] X. Cabré, E. Fontich, and R. de la Llave. The parameterization method for invariant manifolds. III. Overview and applications. J. Differential Equations, 218(2):444–515, 2005. [8] R. Calleja, A. Celletti, and R. de la Llave. A KAM theory for conformally symplectic systems: efficient algorithms and their validation. J. Differential Equations, 255(5):978–1049, 2013. [9] R. Calleja and R. de la Llave. Fast numerical computation of quasi-periodic equilibrium states in 1D statistical mechanics, including twist maps. Nonlin- earity, 22(6):1311–1336, 2009. [10] R. Calleja and J.-Ll. Figueras. Collision of invariant bundles of quasi-periodic attractors in the dissipative standard map. Chaos: An Interdisciplinary Jour- nal of Nonlinear Science, 22(3):033114, 2012. [11] M. Canadell and A. Haro. Parameterization method for computing quasi- periodic reducible normally hyperbolic invariant tori. In F. Casas, V. Martı́nez (eds.), Advances in Differential Equations and Applications, volume 4 of SEMA SIMAI Springer Series. Springer, 2014. [12] M. Canadell and À. Haro. Computation of Quasi-Periodic Normally Hyper- bolic Invariant Tori: Algorithms, Numerical Explorations and Mechanisms of Breakdown. J. Nonlinear Sci., 27(6):1829–1868, 2017. [13] M. Canadell and À. Haro. Computation of Quasiperiodic Normally Hyperbolic Invariant Tori: Rigorous Results. J. Nonlinear Sci., 27(6):1869–1904, 2017. [14] Enric Castellà and Àngel Jorba. On the vertical families of two-dimensional tori near the triangular points of the bicircular problem. Celestial Mech. Dynam. Astronom., 76(1):35–54, 2000. [15] T.N. Chan. Numerical Bifurcation Analysis of Simple Dynamical Systems. PhD thesis, Concordia University, Montreal, Canada, September 1983. [16] B.V. Chirikov. A universal instability of many-dimensional oscillator systems. Phys. Rep., 52(5):264–379, 1979. [17] R. de la Llave. A tutorial on KAM theory. In Smooth ergodic theory and its applications (Seattle, WA, 1999), volume 69 of Proc. Sympos. Pure Math., pages 175–292. Amer. Math. Soc., Providence, RI, 2001. [18] R. de la Llave, A. González, À. Jorba, and J. Villanueva. KAM theory without action-angle variables. Nonlinearity, 18(2):855–895, 2005. [19] R. de la Llave and A. Olvera. The obstruction criterion for non-existence of invarian circles and renormalization. Nonlinearity, 19(8):1907–1937, 2006. [20] D.J. Dichmann, E.J. Doedel, and R.C. Paffenroth. The computation of periodic solutions of the 3-body problem using the numerical continuation software AUTO. In G. Gómez, M. W. Lo, and J. J. Masdemont, editors, Libration Point Orbits and Applications. World Scientific, 2003. [21] D. W. Dunham and R. W. Farquhar. Libration point missions, 1978–2002. In G. Gómez, M. W. Lo, and J. J. Masdemont, editors, Libration Point Orbits and Applications, Singapore, 2003. World Scientific. Proceedings of the conference Libration Point Orbits and Applications, Aiguablava (Girona, Spain), June 10–14, 2002. 52 ALEX HARO AND J.M. MONDELO [22] J.-Ll. Figueras, A. Haro, and A. Luque. Rigorous Computer-Assisted Ap- plication of KAM Theory: A Modern Approach. Found. Comput. Math., 17(5):1123–1193, 2017. [23] Jordi-Lluı́s Figueras and Àlex Haro. Different scenarios for hyperbolicity break- down in quasiperiodic area preserving twist maps. Chaos, 25(12):123119, 16, 2015. [24] E. Fontich, R. de la Llave, and Y. Sire. Construction of invariant whiskered tori by a parameterization method. I. Maps and flows in finite dimensions. J. Differential Equations, 246(8):3136–3213, 2009. [25] C. Gasquet and P. Witomski. Fourier analysis and applications, volume 30 of Texts in Applied Mathematics. Springer-Verlag, New York, 1999. Filtering, nu- merical computation, wavelets, Translated from the French and with a preface by R. Ryan. [26] G. Gómez, À. Jorba, C. Simó, and J. Masdemont. Dynamics and mission design near libration points. Vol. III, volume 4 of World Scientific Monograph Series in Mathematics. World Scientific Publishing Co. Inc., River Edge, NJ, 2001. Advanced methods for collinear points. [27] G. Gómez and J.M. Mondelo. The dynamics around the collinear equilibrium points of the RTBP. Phys. D, 157(4):283–321, 2001. [28] A. González, A. Haro, and R. de la Llave. Singularity theory for non-twist KAM tori. Mem. Amer. Math. Soc., 227(1067):vi+115, 2014. [29] À. Haro, M. Canadell, J.-Ll. Figueras, A. Luque, and J.-M. Mondelo. The pa- rameterization method for invariant manifolds, volume 195 of Applied Math- ematical Sciences. Springer, [Cham], 2016. From rigorous results to effective computations. [30] A. Haro and R. de la Llave. Manifolds on the verge of a hyperbolicity break- down. Chaos, 16(1):013120, 8, 2006. [31] A. Haro and R. de la Llave. A parameterization method for the computation of invariant tori and their whiskers in quasi-periodic maps: numerical algorithms. Discrete Contin. Dyn. Syst. Ser. B, 6(6):1261–1300, 2006. [32] A. Haro and R. de la Llave. A parameterization method for the computation of invariant tori and their whiskers in quasi-periodic maps: rigorous results. J. Differential Equations, 228(2):530–579, 2006. [33] A. Haro and R. de la Llave. A parameterization method for the computation of invariant tori and their whiskers in quasi-periodic maps: explorations and mechanisms for the breakdown of hyperbolicity. SIAM J. Appl. Dyn. Syst., 6(1):142–207 (electronic), 2007. [34] Peter Henrici. Fast fourier methods in computational complex analysis. SIAM Review, 21(4):481–527, 1979. [35] M. Huang, T. Küpper, and N. Masbaum. Computation of invariant tori by the Fourier methods. SIAM J. Sci. Comput., 18(3):918–942, 1997. [36] G. Huguet. The role of hyperbolic invariant objects: from Arnold difussion to biological clocks. PhD thesis, Departament de Matemàtica Aplicada I, Uni- versitat Politècnica de Catalunya, 2008. http://www.ma.utexas.edu/users/ ghuguet/tesi/. [37] G. Huguet, R. de la Llave, and Y. Sire. Computation of whiskered invariant tori and their associated manifolds: new fast algorithms. Discrete Contin. Dyn. Syst., 32(4):1309–1353, 2012. [38] À. Jorba. A methodology for the numerical computation of normal forms, centre manifolds and first integrals of Hamiltonian systems. Experiment. Math., 8(2):155–195, 1999. COMPUTATION OF INVARIANT TORI 53 [39] À. Jorba and E. Olmedo. On the computation of reducible invariant tori on a parallel computer. SIAM J. Appl. Dyn. Syst., 8(4):1382–1404, 2009. [40] Àngel Jorba. Numerical computation of the normal behaviour of invariant curves of n-dimensional maps. Nonlinearity, 14(5):943–976, 2001. [41] Àngel Jorba and Josep J. Masdemont. Dynamics in the center manifold of the Restricted Three–Body Problem. Physica D, 132:189–213, 1999. [42] A. Luque and J. Villanueva. A KAM theorem without action-angle variables for elliptic lower dimensional tori. Nonlinearity, 24(4):1033–1080, 2011. [43] K. R. Meyer, G. R. Hall, and D. Offin. Introduction to Hamiltonian Dynamical Systems and the N –Body Problem. Springer–Verlag, 2nd edition, 2009. [44] J.-M. Mondelo, E. Barrrabés, G. Gómez, and M. Ollé. Fast numerical computa- tion of Lissajous and quasi–halo libration point trajectories and their invariant manifolds. Paper IAC-12,C1,6,9,x14982. Proceedings of the 63rd International Astronautical Congress, 1-5 October 2012, Naples, Italy. [45] A. Olvera and C. Simó. An obstruction method for the destruction of invariant curves. Phys. D, 26(1-3):181–192, 1987. [46] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numeri- cal Recipes in C++: The Art of Scientific Computing. Cambridge University Press, 2nd edition, 2002. [47] F. Schilder, W. Vogt, S. Schreiber, and H.M. Osinga. Fourier methods for quasi-periodic oscillations. Internat. J. Numer. Methods Engrg., 67(5):629– 671, 2006. [48] C. L. Siegel and J. K. Moser. Lectures on celestial mechanics. Classics in Math- ematics. Springer-Verlag, Berlin, 1995. Translated from the German by C. I. Kalme, Reprint of the 1971 translation. [49] C. Simó. Effective computations in celestial mechanics and astrodynamics. In Modern methods of analytical mechanics and their applications (Udine, 1997), volume 387 of CISM Courses and Lectures, pages 55–102. Springer, Vienna, 1998. [50] J. Stoer and R. Bulirsch. Introduction to numerical analysis. Springer, Berlin, second edition, 1993. [51] Victor Szebehely. Theory of orbits. The Restricted Problem of Three Bodies. Academic Press, 1967. [52] Robert L. Warnock. Close approximations to invariant tori in nonlinear me- chanics. Phys. Rev. Lett., 66(14):1803–1806, 1991. Departament de Matemàtiques i Informàtica and BGSMath, Univer- sitat de Barcelona, Gran Via 585, 08007 Barcelona, Spain. Email address: alex@maia.ub.es Departament de Matemàtiques & CERES-IEEC & BGSMath, Univer- sitat Autònoma de Barcelona. Av. de l’Eix Central, Edifici C, 08193 Bellaterra (Barcelona), Spain. Email address: jmm@mat.uab.cat