Elliptic tori in FPU non-linear chains with a small number of nodes
aa r X i v : . [ m a t h - ph ] F e b Elliptic tori in FPU non-linear chainswith a small number of nodes ∗ CHIARA CARACCIOLO
Dipartimento di Matematica, Universit`a degli Studi di Roma “Tor Vergata”,via della Ricerca Scientifica 1, 00133 — Roma (Italy).
UGO LOCATELLI
Dipartimento di Matematica, Universit`a degli Studi di Roma “Tor Vergata”,via della Ricerca Scientifica 1, 00133 — Roma (Italy).e-mails: [email protected], [email protected],
Abstract
We revisit an algorithm constructing elliptic tori, that was originally designed forapplications to planetary hamiltonian systems. The scheme is adapted to properlywork with models of chains of N + 1 particles interacting via anharmonic potentials,thus covering also the case of FPU chains. After having preliminarily settled theHamiltonian in a suitable way, we perform a sequence of canonical transformationsremoving the undesired perturbative terms by an iterative procedure. This is doneby using the Lie series approach, that is explicitly implemented in a programmingcode with the help of a software package, which is especially designed for computeralgebra manipulations. In the cases of FPU chains with N = 4 ,
8, we successfullyapply our new algorithm to the construction of elliptic tori for wide sets of theparameter ruling the size of the perturbation, i.e., the total energy of the system.Moreover, we explore the stability regions surrounding 1D elliptic tori. We compareour semi-analytical results with those provided by numerical explorations of theFPU-model dynamics, where the latter ones are obtained by using techniques basedon the so called frequency analysis. We find that our procedure works up to valuesof the total energy that are of the same order of magnitude with respect to themaximal ones, for which elliptic tori are detected by numerical methods. ∗ Primary: 37J40; Secondary: 37N05, 65P40, 70H08, 70H15.
Key words and phrases:
FPU problem, lower dimensional invariant tori, KAM theory, normal formmethods, perturbation theory for Hamiltonian systems.
C. Caracciolo, U. Locatelli
In the 50’s, Fermi, Pasta and Ulam designed one of the very first numerical experimentsin the history of science (see [15]), realized thanks to the remarkable implementation ofM. Tsingou on the computer Maniac I in the Los Alamos laboratories. It is well knownthat the authors were surprised by the results. In the so called FPU non-linear chains,when just a few slow frequency normal modes are initially excited with relatively lowamplitudes of oscillation, the system does not seem to tend toward the equipartition ofthe energy among all the normal modes themselves. Such a behavior was observed topersist also for very long times of numerical integration and was unexpected, because itis in contrast with the predictions of the Statistical Mechanics. This phenomenon hasbeen explained in terms of metastability (see, e.g., the reviews [1] and [2]). Moreover, thethermalization was observed for larger values of the initial amplitudes of oscillation. Thecomplete understanding of the eventual transition to ergodicity is still an open problem,in particular for what concerns the “thermodinamic limit”, i.e., by keeping fixed the valueof the specific energy E S = E/N with N → ∞ , E being the total energy of the systemand N + 1 the number of nodes making part of the non-linear FPU chain (see, e.g., thosesections in [3] that are devoted to review part of the wide scientific literature on thistopic).The present work is in the spirit of the approach described in [10] and [11], where thepersistence of (meta)stable configurations has been interpreted in terms of the so called q − tori. They are invariant tori of lower dimension, i.e., smaller than the maximal onethat is equal to the number of degrees of freedom N −
1. In the limit of the total energygoing to zero, the motion on q − tori is given by the composition of small oscillations fora reduced number of normal modes. For relatively low values of the initial oscillationsamplitudes, the sharing of the energy with the remaining modes is very limited. For whatconcerns such a phenomenon of (so called) localization, the behavior of the motions on q − tori is very similar to those considered in the pioneering numerical simulations of theFPU chains. Therefore, a special attention is devoted to the linear stability of q − tori,because the concept of similarity between those orbits would be obviously reinforced,if the initial conditions with just a few modes excited would be in a sort of stabilityregion around the q − tori. In [10] and [11], q − tori are determined by using the method ofPoincar´e–Lindstedt series, that is able to compute the shape of invariant manifolds andthe motion on them; however, in order to have information about the local dynamics inthe neighborhood of those tori, some additional work is required (see., e.g., [8]). Let usalso recall that q − tori are a generalization of q − breathers, that are periodic orbits whichcan be obtained as a continuation of the small oscillations involving one single normalmode. Their role in the dynamics of FPU chains has been widely investigated startingby a sequence of articles written by S. Flach and some of his coworkers (see [44], [17] andtherein references).In KAM theory, linearly stable q − tori are commonly known as elliptic tori, that arecompact invariant manifold of dimension smaller than the number of degrees of freedom.Moreover, they are such that the transverse motions are given by a superposition ofsmall oscillations (as in the case of the elliptic equilibrium points). Our main goal is to lliptic tori in FPU non-linear chains . . . q − breathers.Let us emphasize that our approach, being based on a normal form method, is in avery good position for this type of investigation, because it directly gives the equations ofmotion ruling the dynamics in the neighborhood of the elliptic tori without any additionalwork. In the present paper we mainly focus on the comparisons between the resultsprovided by the explicit construction of normal forms with those given by a suitablestudy of the numerical integrations. The accuracy of the former computational method(that is often said of semi-analytic type) is strongly limited when the number of degreesof freedom is increased. This is the reason why we prefer to consider non-linear chainswith small numbers of nodes, in order to obtain sharp results eventually giving a clearvalidation of our approach. In this respect, we have to admit that the spirit of our workis evidently far from the original motivations making so interesting the study of the FPUmodel, because they mainly concern with the pathway to the statistical equilibrium andthe eventual occurrence of obstacles represented by metastable configurations.The paper is organised as follows. In section 2, we introduce the basic settings in aHamiltonian framework that are needed to properly define the FPU model. The algorithmconstructing the normal form for elliptic tori is described in section 3. As a furtherdifference with respect to the works on q − tori, where the numerical integrations wereanalyzed by using the Generalized ALignment Indices method (GALI for short, see [50]),here the results provided by the semi-analytic method are compared with those obtainedby applying Frequency Analysis (see, e.g., [35]). The adaptation of this latter technique ofnumerical investigation to the present context is described in section 4. Most of the resultstesting both the validity of our approach and its performances are collected in section 5. C. Caracciolo, U. Locatelli
Section 6 is devoted to the introduction of suitable Birkhoff normal forms centered aroundthe invariant elliptic tori; this allows to study the dynamics in the vicinity of them. Boththe conclusions and the main perspectives are discussed in section 7.
The Hamiltonian describing a FPU chain of N + 1 particles is the following: H ( y , x ) = 12 N − X j =0 (cid:2) y j + ( x j +1 − x j ) ) (cid:3) + α N − X j =0 ( x j +1 − x j ) + β N − X j =0 ( x j +1 − x j ) , (1)with α, β ∈ R , where x k is the displacement of the k -th particle with respect to theequilibrium position and y k its canonically conjugate momentum. Usually, N = 2 p andboundary conditions are fixed so that x = x N = 0 , y = y N = 0. The case α = 0 , β = 0is called α − model; when α = 0 , β = 0 it is called β − model.Using the variables ( Y , X ), that are related to the normal modes of oscillation and aredefined by the following canonical transformation ∀ l = 1 , . . . , N − x l = r N N − X j =1 X j √ ν j sin (cid:18) jlπN (cid:19) , y l = r N N − X j =1 Y j √ ν j sin (cid:18) jlπN (cid:19) , (2)with ν j = 2 sin (cid:18) jπ N (cid:19) , (3)the Hamiltonian can be rewritten as follows: H ( Y , X ) = 12 N − X j =1 h ν j (cid:0) X j + Y j (cid:1)i + H ∗ ( X ) , (4)with H ∗ containing only cubic and quartic terms in X .Now, it is convenient to temporarily introduce action-angle coordinates adapted to afixed number of n ≤ N − Y j = p I j cos ϕ j , X j = p I j sin ϕ j , ∀ j = 1 , . . . , n . This leads Hamiltonian (4) to the following form: H ( I , ξ , ϕ , η ) = n X j =1 ω j I j + n X j =1 Ω j ( ξ j + η j )2 + H ∗ ( I , ξ , ϕ , η ) , (5)where we have split ν as ν = ( ω , Ω ) ∈ R n × R n , with n = N − − n , and we haveintroduced ξ j = Y j + n , η j = X j + n , ∀ j = 1 , . . . , n . lliptic tori in FPU non-linear chains . . . I ⋆ of the n –dimensional action vector I and then we perform a translation, by defining p j = I j − I ⋆j ∀ j = 1 , . . . , n , (6)and leaving the other coordinates unchanged. Finally, we expand the Hamiltonian inpower series of p , ξ , η and, renaming the angles ϕ as q , in Fourier series in the angles q ,so that the Hamiltonian will assume the following form: H ( p , q , ξ , η ) = ω (0) · p + n X j =1 Ω (0) j ξ j + η j ) + F ( p , q , ξ , η ) . (7)Let us stress that in the expression of the Hamiltonian above, we have omitted the para-metric dependency with respect to the initial translation vector I ⋆ , in spite of its crucialimportance, that we are going to discuss in this final part of the present section. Asa first remark, let us recall that in the phase space the distance from the equilibriumpoint corresponding to the origin of the canonical coordinates I , ξ and η is essentiallygoverned by the initial translation vector I ⋆ and so also for what concerns the width ofthe oscillations in the normal mode variables ( Y , X ). From an heuristic point of view,this allows us to see that the size of the perturbations is ruled by I ⋆ .Hereafter, in the notation we will stress the parametric dependency on the initialtranslation I ⋆ (actually, on the angular velocity vector ω (0) that can be put in a bijectivecorrespondence with I ⋆ ) just when it is really necessary for the understanding. In partic-ular, this will be mandatory in section 3.4, where the formal statement of theorem 3.1 isreported. In fact, that rigorous result claims that our algorithm constructing elliptic tori(which is described in the next section) is converging just in a subset of an open ball of theinitial vectors I ⋆ (cid:0) ω (0) (cid:1) . Under suitable conditions, such a subset has positive Lebesguemeasure which is larger and larger (in a probabilistic sense) for perturbations shrinkingto zero. We present the formal algorithm making reference to a generic Hamiltonian with n + n degrees of freedom, where the canonical coordinates ( p , q , ξ , η ) can be naturally split intwo parts, that are ( p , q ) ∈ R n × T n , action-angle variables in the neighborhood of thewanted torus, and ( ξ , η ) ∈ R n × R n , in the neighborhood of the origin.Let us introduce some notation. For some fixed positive integer K we introduce thedistinguished classes of functions b P sK ˆ m, ˆ ℓ , with integers ˆ m, ˆ ℓ, s ≥ Actually, when the whole procedure constructing the invariant elliptic torus is successful, on thatsame torus the motion law of the action vector I will be such that k I ( t ) − I ⋆ k = o ( k I ⋆ k ). Reasonableways to practically choose the initial translation I ⋆ will be described in section 5. C. Caracciolo, U. Locatelli g ∈ b P sK ˆ m, ˆ ℓ can be written as g ( p , q , ξ , η ) = X m ∈ N n | m | = ˆ m X ( ℓ , ¯ ℓ ) ∈ N n | ℓ | + | ¯ ℓ | =ˆ ℓ X k ∈ Z n | k |≤ sK p m ξ ℓ η ¯ ℓ (cid:2) c m , ℓ , ¯ ℓ , k cos( k · q ) + d m , ℓ , ¯ ℓ , k sin( k · q ) (cid:3) , (8)with the coefficients c m , ℓ , ¯ ℓ , k , d m , ℓ , ¯ ℓ , k ∈ R ; in the previous formula, we have introducedthe symbol | · | to denote the ℓ -norm and we have adopted the multi-index notation,i.e., p m = Q n j =1 p m j j . We stress that the positive integer K is chosen in such a way toexploit the Fourier decay of the coefficients in the analytic functions. As it is detailedin subsection 4.4 of [24], this is made in order to split all the Hamiltonian in blocks ofdifferent orders of magnitude with respect the natural small parameter related to thequasi-integrable problem that is the object of our study. To fix the ideas, a suitablechoice of the parameter K allows to obtain initial estimates of the type described into thefundamental assumption (c) of theorem 3.1. We say that g ∈ P sKℓ if g ∈ [ ˆ m ≥ , ˆ ℓ ≥
02 ˆ m +ˆ ℓ = ℓ b P sK ˆ m, ˆ ℓ . (9)Hereafter, we will omit the dependence of the function from the variables, unless it hassome special meaning. Moreover, we will adopt the usual notation for the average of afunction g with respect to the generic angles ϑ ∈ T n , i.e., h g i ϑ = R T n d ϑ . . . d ϑ n g/ (2 π ) n .Let us consider a Hamiltonian that can be written in the following way: H (0) ( p , q , ξ , η ) = ω (0) · p + n X j =1 Ω (0) j ξ j + η j ) + X s ≥ X ℓ ≥ f (0 , s ) ℓ ( p , q , ξ , η )+ X s ≥ X ℓ =0 f (0 , s ) ℓ ( p , q , ξ , η ) , (10)where f (0 , s ) ℓ ∈ P sKℓ , while the first upper index is related to the normalization step.Starting from the Hamiltonian described in (1) and following the prescriptions givenin the previous section, one can easily check that a FPU chain of N + 1 particles canbe brought to the form above. In other words, the Hamiltonian in formula (7) can beexpanded as H (0) in (10), where f (0 , s ) ℓ = 0 when s ≥ f (0 , ℓ ∈ P Kℓ , f (0 , ℓ ∈ P Kℓ ∀ ℓ ≥
0, with K = 2. This holds true, both for the α –model and the β one.Our main purpose is to eliminate from the Hamiltonian all the terms having totaldegree less than three in the square root of the actions; by referring to the paradigmatic Setting K = 2 is quite natural for Hamiltonian systems close to stable equilibria, see, e.g., [27]. Thereason of such a choice can be understood by referring to the discussion in Section 2. The size of theperturbing terms making part of the Hamiltonian term H ∗ in (5) is O (cid:0) k I ⋆ k s/ (cid:1) , s being the degree inthe square root of the actions I . Each monomial term appearing in the perturbation will have a maximaltrigonometric degree in ϕ that is not greater than s . The subsequent introduction of the actions p modifies the dependency on the actions, but does not change the relation between the order of magnitudein the (small) parameter k I ⋆ k and the degree in the angles, the latter being at most twice bigger withrespect to the former. lliptic tori in FPU non-linear chains . . . H ( ∞ ) ( P , Q , Ξ , Θ ) = ω ( ∞ ) · P + n X j =1 Ω ( ∞ ) j j + Θ j ) + X s ≥ X ℓ ≥ f ( ∞ , s ) ℓ ( P , Q , Ξ , Θ ) + E ( ∞ ) , (11)with f ( ∞ , s ) ℓ ∈ P sKℓ and E ( ∞ ) ∈ P , i.e., it is a constant. It is easy to check that( P ( t ) , Q ( t ) , Ξ ( t ) , Θ ( t )) = ( , Q + ω ( ∞ ) t, , ) (12)is a solution of Hamilton’s equations related to (11), since that Hamiltonian, except forits first part, contains terms of type O ( k P k ), O ( k P kk ( Ξ , Θ ) k ) and O ( k ( Ξ , Θ ) k ) only.Therefore, the motion law (12), that is generated by the initial condition ( , Q , , ), isquasi-periodic with an angular velocity vector equal to ω ( ∞ ) and the corresponding orbitlies on the n − dimensional invariant torus P = , Ξ = Θ = . The energy level ofsuch a manifold is H ( ∞ ) ( , Q , , ) = E ( ∞ ) . Moreover, it is elliptic in the sense that thetransverse dynamics in a neighborhood of the invariant torus itself is given by oscillationswhose corresponding angular velocity vector is approaching Ω ( ∞ ) in the limit of k ( Ξ , Θ ) k going to zero.The formal algorithm for the construction of the normal form is composed by a se-quence of canonical transformations, defined using the formalism of Lie series. We cansummarize the r –th normalization step, by giving the formula defining the canonicalchange of coordinates that transforms the intermediate Hamiltonian H ( r − into the sub-sequent H ( r ) , i.e., exp (cid:0) L χ ( r )0 (cid:1) ◦ exp (cid:0) L χ ( r )1 (cid:1) ◦ exp (cid:0) L χ ( r )2 (cid:1) ◦ D ( r ) , (13)where L f · = {· , f } is the Lie derivative operator related to the Poisson bracket {· , ·} andthe Lie series operator exp (cid:0) L χ ( r ) j (cid:1) · = ∞ X i =0 i ! L iχ ( r ) j · (14)removes the Hamiltonian terms with total degree in the square root of the actions equalto j and with trigonometric degree in the angles q up to rK . Moreover, by a linearcanonical transformation D ( r ) , the terms that are quadratic in ( ξ , η ) and do not dependon the angles q are brought to a diagonal form. At the end of this r –th normalization Because of the so called “exchange theorem” (see [29]), the new Hamiltonian H ( r ) is obtained fromthe old one, by applying the Lie series to H ( r − in reverse order with respect to what is written in (13).This is consistent with the order of the discussion in the following subsections: the first stage of the r –thnormalization step deals with the canonical transformation generated by χ ( r )0 , the second one with χ ( r )1 and the last one with both χ ( r )2 and D ( r ) . A rather self-consistent introduction to the Lie series formalismin the Hamiltonian framework can be found, e.g., in [22], where the method is described in an explicitlyalgorithmic way, going back to the seminal articles written by Hori, Deprit and Henrard (see [32], [13]and [31], resp.). C. Caracciolo, U. Locatelli step, the ineliminable terms that are independent on the angles q and either linear in p or quadratic in ( ξ , η ) are added to the normal form part. This requires to update theangular velocities from (cid:0) ω ( r − , Ω ( r − (cid:1) to (cid:0) ω ( r ) , Ω ( r ) (cid:1) , that is why in (11) the Hamiltonianin normal form has new frequency vectors ω ( ∞ ) and Ω ( ∞ ) .Now, let us suppose that the Hamiltonian at step r − H ( r − ( p , q , ξ , η ) = ω ( r − · p + n X j =1 Ω ( r − j ξ j + η j ) + X s ≥ r X ℓ =0 f ( r − , s ) ℓ ( p , q , ξ , η )+ X s ≥ X ℓ ≥ f ( r − , s ) ℓ ( p , q , ξ , η ) + E ( r − , (15)where f ( r − , s ) ℓ ∈ P sKℓ and E ( r − ∈ P , i.e., it is a constant gathering all the terms whichcannot be eliminated by the first homological equations that have already been solvedfor all the normalization steps 1 , . . . , r − H (0) written in formula (10) is in form (15), with r = 1 and E (0) = 0. In the following subsections we are going to detail how the algorithmactually works. r –th normalization step In the context of the r -th normalization step, the first stage aims to remove the termsdepending just on the angles q up to the trigonometrical degree rK , i.e., the termscollected in f ( r − ,r )0 . We determine the generating function χ ( r )0 by solving the homologicalequation n χ ( r )0 , ω ( r − · p o + f ( r − , r )0 ( q ) − h f ( r − ,r )0 ( q ) i q = 0 . (16)Let us remark that, since f ( r − , r )0 depends on q only, h f ( r − , r )0 i q is a constant term thatgives its contribution to the energy level of the final invariant elliptic torus; thus, weintroduce E ( r ) = E ( r − + h f ( r − , r )0 i q . (17)Starting from the Fourier expansion f ( r − ,r )0 ( q ) = X | k |≤ rK (cid:2) c k cos( k · q ) + d k sin( k · q ) (cid:3) , (18)with c k , d k real known coefficients, we readily get the following expression for the gener-ating function χ ( r )0 ( q ) = X < | k |≤ rK (cid:20) − c k sin( k · q ) k · ω ( r − + d k cos( k · q ) k · ω ( r − (cid:21) . (19)Therefore, the homological equation (16) admits a solution for the generating function χ ( r )0 , provided that the frequency vector ω ( r − is non-resonant up to the trigonometricdegree rK . To fix the ideas, we assume that it fulfills the following Diophantine condition:min < | k |≤ rK | k · ω ( r − | ≥ γ ( rK ) τ . (20) lliptic tori in FPU non-linear chains . . . L χ ( r )0 to the Hamiltonian at step r −
1, renaming the new variables as the old ones (by abuse of notation), the transformedHamiltonian reads as H (I; r ) = exp (cid:16) L χ ( r )0 (cid:17) H ( r − = ω ( r − · p + n X j =1 Ω ( r − j ( ξ j + η j )2 + X s ≥ r X ℓ =0 f (I; r, s ) ℓ + X s ≥ X ℓ ≥ f (I; r, s ) ℓ + E ( r ) , (21)where the functions f (I; r, s ) ℓ ∈ P sKℓ contributing to the definition of the new Hamiltonianare recursively defined as follows: f (I; r, r )0 = 0 ,f (I; r, r + m )0 = f ( r − , r + m )0 for 0 < m < r ,f (I; r, s ) ℓ = ⌊ s/r ⌋ X j =0 i ! L iχ ( r )0 f ( r − , s − ir ) ℓ +2 i for ℓ = 0 , s ≥ r or ℓ = 1 , , s ≥ r or ℓ ≥ , s ≥ . (22)Indeed, in order to practically implement the definitions included in formula (22), wehave found convenient to proceed in the following way. First, we define the new termsas the old ones, i.e., f (I; r, s ) ℓ = f ( r − , s ) ℓ ; then, each term generated by Lie derivatives withrespect to the generating function is added to the corresponding class. This is made by asequence of redefinitions, for each of them an abuse of notation has to be tolerated. Sincethe Lie derivative of the generating function χ ( r )0 decreases by 1 the degree in p , while thetrigonometrical degree in the angles q is increased by rK , we set f (I; r, s + ir ) ℓ − i ← ֓ i ! L iχ ( r )0 f ( r − , s ) ℓ ∀ ℓ ≥ , ≤ i ≤ ⌊ ℓ/ ⌋ , s ≥ , (23)where with the notation a ← ֓ b we mean that the quantity a is redefined so as to be equalto a + b . Let us recall that f (I; r, r )0 = 0, because of the homological equation (16). r –th normalization step The second stage of the r -th normalization step acts on the Hamiltonian that is initiallyexpanded as in (21), in order to remove the perturbing terms linear in ( ξ , η ) and inde-pendent of p , that are those collected in f (I; r, r )1 . Thus, we have to solve the followinghomological equation: ( χ ( r )1 , ω ( r − · p + n X j =1 Ω ( r − j (cid:0) ξ j + η j (cid:1)) + f (I; r, r )1 ( q , ξ , η ) = 0 . (24)In order to solve such an equation, it is convenient to temporarily replace the transversevariables ( ξ , η ) with action-angle canonical coordinates, by defining ξ j = p J j cos( ϕ j )0 C. Caracciolo, U. Locatelli and η j = p J j sin( ϕ j ). Let us write the expansion of f (I; r, r )1 ( q , J , ϕ ) as follows: f (I; r, r )1 ( q , J , ϕ ) = X ≤ k ≤ rK n X j =1 p J j h c ( ± ) k , j cos( k · q ± ϕ j ) + d ( ± ) k , j sin( k · q ± ϕ j ) i , (25)where c ( ± ) k , j and d ( ± ) k , j are some real coefficients. Therefore, the generating function χ ( r )1 solving equation (24) is determined in such a way that χ ( r )1 ( q , J , ϕ ) = X ≤ k ≤ rK n X j =1 p J j " − c ( ± ) k , j sin( k · q ± ϕ j ) k · ω ( r − ± Ω ( r − j + d ( ± ) k , j cos( k · q ± ϕ j ) k · ω ( r − ± Ω ( r − j . (26)This expression is well-defined, provided that the frequency vector ω ( r − satisfies theso-called first Melnikov non-resonance condition up to order rK (see [38]), i.e.,min < | k |≤ rK, | ℓ | =1 (cid:12)(cid:12) k · ω ( r − + ℓ · Ω ( r − (cid:12)(cid:12) ≥ γ ( rK ) τ and min | ℓ | =1 (cid:12)(cid:12) ℓ · Ω ( r − (cid:12)(cid:12) ≥ γ , (27)for some fixed values of both γ > τ > n −
1. By applying the Lie series exp( L χ ( r )1 ) tothe old Hamiltonian H (I; r ) , we have a new one, i.e., H (II; r ) = exp (cid:0) L χ ( r )1 (cid:1) H (I; r ) , where nowthe generating function χ ( r )1 is reconsidered to be dependent on the variables ( q , ξ , η ). Theexpansion of the new Hamiltonian will have exactly the same structure as that describedin (21), with f (II; r, s ) ℓ in place of f (I; r, s ) ℓ and f (II; r, r )1 = 0. The new terms f (II; r, s ) ℓ thatcompose the Hamiltonian can be determined with calculations similar to those listedduring the description of the first stage of normalization. For what concerns the explicitcomputer algebra manipulations, it is enough to remark that the Lie derivative withrespect to χ ( r )1 decreases by 1 the total degree in the square root of the actions. Thisexplains why it is convenient to set the sequence of redefinitions of the new terms f (II; r, s ) ℓ in such a way that f (II; r, s + ir ) ℓ − i ← ֓ i ! L iχ ( r )1 f (I; r, s ) ℓ ∀ ℓ ≥ , ≤ i ≤ ℓ , s ≥ . (28)This also means that f (II; r, r )1 = 0, due to the homological equation (24). r –th normalization step The third and last stage of normalization is more elaborated. It aims to remove terms be-longing to two different classes: first, those linear in p and independent of ( ξ , η ), moreover,other terms that are quadratic in ( ξ , η ) and independent of p . Such a part of the perturb-ation is removed by the composition of two canonical transformations expressed by Lieseries, with generating functions X ( r )2 ( p , q ) ∈ b P rK , and Y ( r )2 ( q , ξ , η ) ∈ b P rK , , respectively.Moreover, the third stage is ended by a linear canonical transformation D ( r )2 that leavesthe pair ( p , q ) unchanged and it aims to diagonalize the terms that are quadratic in ( ξ , η ) lliptic tori in FPU non-linear chains . . . q . Let us detail all these changes of coordinates, so thatthe algorithm is unambiguously defined.The generating functions X ( r )2 is in charge to remove terms that are linear in p anddo depend on the angles q up to the trigonometric degree rK . Therefore, it is a solutionof the following homological equation: n X ( r )2 , ω ( r − · p o + f (II; r, r )2 ( p , q ) − h f (II; r, r )2 ( p , q ) i q = 0 . (29)Let us recall that f (II; r, r )2 ∈ P rK = b P rK , ∪ b P rK , ; in general, such a function really dependson all the canonical variables, i.e., f (II; r, r )2 = f (II; r, r )2 ( p , q , ξ , η ). Therefore, we denote with f (II; r, r )2 ( p , q ) the subpart of f (II; r, r )2 that is depending just on ( p , q ). Analogously, in thefollowing f (II; r, r )2 ( q , ξ , η ) will denote the subpart of f (II; r, r )2 that does depend on all thecanonical variables but the actions p and so on also for what concerns f (II; r, r )2 ( ξ , η ). Thishighly non-standard notation will be adopted all along the present subsection. Let us hereemphasize that the term h f (II; r, r )2 ( p , q ) i q will be added to the part in normal form, with achange of the angular velocity vector ω ( r − . We can write the expansion of f (II; r, r )2 ( p , q )as f (II; r, r )2 ( p , q ) = X < k ≤ rK n X j =1 p j [ c k , j cos( k · q ) + d k , j sin( k · q )] , (30)thus, the corresponding solution of the homological equation (29) is given by X ( r )2 ( p , q ) = X < k ≤ rK n X j =1 p j (cid:20) − c k , j sin( k · q ) k · ω ( r − + d k , j cos( k · q ) k · ω ( r − (cid:21) . (31)Let us remark that X ( r )2 is well-defined provided that the frequency vector ω ( r − satisfiesthe non-resonance condition already reported in (20).The generating function Y ( r )2 aims to remove the part of the term of f (II; r, r )2 that isquadratic in ( ξ , η ) and does depend on the angles q . Therefore, Y ( r )2 has to solve thefollowing homological equation: ( Y ( r )2 , ω ( r − · p + n X j =1 Ω ( r ) j (cid:0) ξ j + η j (cid:1)) + f (II; r, r )2 ( q , ξ , η ) − h f (II; r, r )2 ( q , ξ , η ) i q = 0 . (32)It is convenient to defer the discussion about the treatment of the term h f (II; r, r )2 ( q , ξ , η ) i q ,because it is rather peculiar. Now, let us temporarily replace the transverse canonicalvariables ( ξ , η ) with the action-angle coordinates ( J , ϕ ) (exactly in the same way as inthe previous section 3.2), so as to write the expansion of the perturbing term we nowwant to remove in the following way: f (II; r, r )2 ( q , J , ϕ ) = X < k ≤ rK n X i, j =1 p J i J j h c ( ± , ± ) k , i, j cos( k · q ± ϕ i ± ϕ j )+ d ( ± , ± ) k , i, j sin( k · q ± ϕ i ± ϕ j ) i . (33)2 C. Caracciolo, U. Locatelli
Thus, the generating function Y ( r )2 is determined by equation (32) in such a way that Y ( r )2 ( q , J , ϕ ) = X < k ≤ rK n X i, j =1 p J i J j " − c ( ± , ± ) k , i, j sin( k · q ± ϕ i ± ϕ j ) k · ω ( r − ± Ω ( r − i ± Ω ( r − j + d ( ± , ± ) k , i, j cos( k · q ± ϕ i ± ϕ j ) k · ω ( r − ± Ω ( r − i ± Ω ( r − j , (34)that is well defined provided that the angular velocity vector ω ( r − satisfies both thealready mentioned Diophantine inequality (20) and the so-called second Melnikov non-resonance condition up to order rK (see [38]), i.e.,min < | k |≤ rK, | ℓ | =2 (cid:12)(cid:12) k · ω ( r − + ℓ · Ω ( r − (cid:12)(cid:12) ≥ γ ( rK ) τ (35)with fixed values of both the parameters γ > τ > n − p or quadratic in ( ξ , η ) and do not depend on q . The former ones can bedirectly added to the part in normal form, whereas the latter have to be preliminarily putin diagonal form. This can be done by a linear canonical transformation D ( r ) , that canbe computed in a direct way so that n X j =1 Ω ( r − j ( ξ j + η j )2 + f (II; r, r )2 ( ξ , η ) ! (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( ξ , η )= D ( r ) ( ¯ ξ , ¯ η ) = n X j =1 Ω ( r ) j ( ¯ ξ j + ¯ η j )2 , (36)as it is explained, e.g., in section 7 of [23]. Such a direct calculation can be performedprovided that min | ℓ | =2 (cid:12)(cid:12) ℓ · Ω ( r − (cid:12)(cid:12) ≥ γ . (37)Finally, we need to understand how these generating functions affect the Hamiltonianterms. Let us recall that the Poisson brackets with the generating functions X ( r )2 and Y ( r )2 do not change the total degree in the square root of the actions. In order to describe thedefinitions of the new Hamiltonian terms it is convenient to introduce the intermediatefunctions g ( r, s ) ℓ , g ′ ( r, s ) ℓ in the following way. First, we define g ( r, s ) ℓ = f (II; r, s ) ℓ for all non-negative values of the indexes ℓ , r and s ; then, we consider the effects induced by theapplication of the Lie series with generating function X ( r )2 to the Hamiltonian. In orderto do that, we gather the new terms in g ( r, s ) ℓ according to both their total degree in the In practical implementations, such a change of coordinates D ( r ) can be conveniently defined bycomposing a subsequence of Lie series, each of them being related to a quadratic generating function D ( r ; m )2 ( ξ , η ). By the way, this allows to automatically determine D ( r ) in such a way it is a near to theidentity canonical transformation. Moreover, as a further alternative method, when one is dealing withthe estimates needed to prove the convergence of the algorithm, in [26] the use of the Lie transforms (thatare equivalent to the composition of infinite sequences of Lie series) has been found to be very suitable. lliptic tori in FPU non-linear chains . . . g ( r, s + ir ) ℓ ← ֓ i ! L iX ( r )2 f (II; r, s ) ℓ ∀ i ≥ , ℓ ≥ , s ≥ . (38)In the same way, we gather in g ′ ( r, s ) ℓ all the Hamiltonian terms created by the applicationof the Lie series having Y ( r )2 as generating function. This is summarized by the followingredefinitions: g ′ ( r, s ) ℓ = g ( r, s ) ℓ and g ′ ( r, s + ir ) ℓ ← ֓ i ! L iY ( r )2 g ( r, s ) ℓ ∀ i ≥ , ℓ ≥ , s ≥ . (39)Let us remark that each class of type b P sℓ ,ℓ or P sℓ is preserved by the diagonalizationlinear operator D ( r ) , for all non-negative values of the indexes ℓ , ℓ , ℓ and s . Therefore,it is natural to put f ( r,s ) ℓ = g ′ ( r, s ) ℓ ◦ D ( r ) . (40)for all indexes ℓ ≥ s ≥ r –th normalization step, it is convenient that the terms linearlydepending just on p or J are included in the main part of the Hamiltonian, because allof them belong to the same class of functions, i.e. P . For this purpose, we introduce thenew angular velocity vector ω ( r ) , in such a way that ω ( r ) · p = ω ( r − · p + f (II; r, ( p ) , (41)while the new values of the components of Ω ( r ) are defined by equation (36), that alsoallows us to put f ( r,r )2 = 0.The Hamiltonian at the end of the r -th normalization step H ( r ) has the same structureof H ( r − in (15), but the new perturbative terms f ( r,s ) ℓ with ℓ = 0 , , rK . The non-resonance conditions we have assumed in (20), (27), (35) and (37) can be sum-marized in the following way:min < | k |≤ rK, ≤| ℓ |≤ (cid:12)(cid:12) k · ω ( r − + ℓ · Ω ( r − (cid:12)(cid:12) ≥ γ ( rK ) τ and min < | ℓ |≤ (cid:12)(cid:12) ℓ · Ω ( r − (cid:12)(cid:12) ≥ γ , (42)with γ > τ > n −
1. Let us here resume the parametric dependence of all theHamiltonian terms on the initial value of the angular velocity vector ω (0) , that is in abijective correspondence with the initial translation vector I ⋆ , as it has been discussedat the end of sect. 2. In particular, in the Diophantine inequalities reported in (42)4 C. Caracciolo, U. Locatelli the angular velocity vectors at the r -th normalization step are functions of ω (0) , i.e., ω ( r − = ω ( r − ( ω (0) ) and Ω ( r − = Ω ( r − ( ω (0) ). Let us recall that we do not try tokeep a full control on the way ω ( r − ( ω (0) ) and Ω ( r − ( ω (0) ) are modified passing from the r –th normalization step to the next one. Such an approach is in contrast with what isusually done to construct the Komogorov normal form for maximal invariant tori, wherethe angular velocities are kept fixed (see [33] or, e.g., [24]), but it is somehow unavoidablebecause of the occurrence of the transversal angular velocities Ω ( r − ( ω (0) ) that in generalcannot remain constant along the normalization procedure. This seems to prevent thecomplete construction of the normal form and, therefore, the proof of the existence of anelliptic torus. Nevertheless, following the approach designed by P¨oschel in [45], it can beproved that the Lebesgue measure of the resonant regions where the Melnikov conditionsare not satisfied shrinks to zero with the size of the perturbation. Therefore, an analysismade with a suitable scheme of estimates allows to prove the convergence of the algorithm,as it is summarized by the statement below. Theorem 3.1
Consider the following family of real Hamiltonians: H (0) ( p , q , ξ , η ; ω (0) ) = ω (0) · p + n X j =1 Ω (0) j ( ω (0) )2 ( ξ j + η j )++ X s ≥ X ℓ ≥ f (0 , s ) ℓ ( p , q , ξ , η ; ω (0) ) + X s ≥ X ℓ =0 f (0 , s ) ℓ ( p , q , ξ , η ; ω (0) ) , (43) where ( p , q , ξ , η ) ∈ O × T n ×O with O and O open neighborhoods of the origin in R n and R n , respectively, while ω (0) ∈ U , U being an open subset of R n , and f (0 , s ) ℓ ∈ P sKℓ .Moreover, we assume that:(a) all the functions Ω j : U 7→ R and f (0 , s ) ℓ : O × T n × O × U 7→ R , appearing in (43) ,are analytic functions with respect to ω (0) ∈ U ;(b) Ω (0) i ( ω (0) ) = Ω (0) j ( ω (0) ) and Ω (0) i ( ω (0) ) = 0 for ω (0) ∈ U and ≤ i < j ≤ n , ≤ i ≤ n ;(c) for some fixed and positive values of ε and E , one has sup ( p , q , ξ , η ; ω (0) ) ∈O × T n ×O ×U (cid:12)(cid:12)(cid:12) f (0 ,s ) ℓ ( p , q , ξ , η ; ω (0) ) (cid:12)(cid:12)(cid:12) ≤ ε s E (44) ∀ s ≥ , ℓ ≥ and ∀ ℓ ≥ s = 0 .Then, there is a positive ε ⋆ such that for ≤ ε < ε ⋆ the following statement holds true:there exists a non-resonant set U ( ∞ ) ⊂ U of positive Lebesgue measure and with themeasure of U \ U ( ∞ ) tending to zero for ε → for bounded U , such that for each ω (0) ∈U ( ∞ ) there exists an analytic canonical transformation ( p , q , ξ , η ) = ψ ( ∞ ) ω (0) ( P , Q , Ξ , Θ ) leading the Hamiltonian to the normal form H ( ∞ ) ( P , Q , Ξ , Θ ; ω (0) ) = ω ( ∞ ) · P + n X j =1 Ω ( ∞ ) j (cid:0) Ξ j + Θ j (cid:1) X s ≥ X ℓ ≥ f ( ∞ , s ) ℓ ( P , Q , Ξ , Θ ; ω (0) ) + E ( ∞ ) , (45) lliptic tori in FPU non-linear chains . . . where ω ( ∞ ) = ω ( ∞ ) ( ω (0) ) , Ω ( ∞ ) = Ω ( ∞ ) ( ω (0) ) , f ( ∞ , s ) ℓ ∈ P sKℓ and E ( ∞ ) = E ( ∞ ) ( ω (0) ) isa finite real value fixing the constant energy level that corresponds to the invariant torus (cid:8) ( P = , Q ∈ T n , Ξ = , Θ = ) (cid:9) . Actually, the statement above does not add anything new with respect to that reportedin [45]. In spite of the fact that, here, we do not aim to deal with another purely analyticaldemonstration of the existence of elliptic tori, what really matters is in the proof oftheorem 3.1, because it is based on a formal algorithm that is substantially the same (apartvery minor differences) with respect to that described in subsections 3.1–3.3. Therefore,the theorem above rigorously ensures the convergence of the whole procedure, providedthe perturbation is small enough. Such an approach, that goes back to the proof of KAMtheorem using classical series expansions, has been showed to be in a very good positionfor the applications to realistic models (see, e.g., [24], [19], [27]), while the proof adoptedin [45] is based on a fast convergence scheme of quadratic type (a so called Newton-likemethod).Hereafter, we abandon mainly theoretical arguments; in particular, in section 5, ournormal form algorithm will be used in practice, for explicitly constructing elliptic tori inthe FPU model. In particular, since we are not going to apply the theorem, we do notcare about the threshold value ε ⋆ , that is extremely low. Actually, it is well known thatcomputer assisted proofs can ensure the existence of invariant KAM surfaces for valuesof the small parameter that are very close to its upper limit; moreover, there are versionsof this kind of rigorous results that are suitable for elliptic tori (see [16] and [37]). Sincewe just want to extract from the proof of theorem 3.1 a key element, i.e., the widelydescribed constructive algorithm, and use it in practice , we will limit ourselves to checkfor its convergence with a numerical criterion that will be shown later. Several powerful numerical approaches have been designed in order to explore the dynam-ics of Hamiltonian systems. In the present context, we think that the so called FrequencyAnalysis Method (hereafter FA, see [35] for an introduction) is in a better position, be-cause it can be efficiently adapted to the search for elliptic tori (see, e.g., [42] and [12] forapplications to problems of interest in Astronomy).
Following the usual prescriptions provided in the framework of FA, we study complexsignals, that are formed by pairs of real canonical variables of cartesian type, namely weconsider the motion laws t Y j ( t ) + i X j ( t ), where j = 1 , . . . , N − Y , X )are related to the normal modes of oscillation (recall formulæ (2)–(4)). Actually, we focus A complete proof is included in C. Caracciolo:
On the stability in the neighborhood of invariantelliptic tori , Ph.D. thesis, Univ. of Rome “Tor Vergata” (2021), which is available on request to theauthor. C. Caracciolo, U. Locatelli on the following decomposition of the signals: Y j ( t ) + i X j ( t ) ≃ N C X s =1 A j,s e i( υ j,s t + φ j,s ) , (46)where N C is the number of the components we want to consider, while A j,s ∈ R + and φ j,s ∈ [0 , π ) are the amplitude and the initial phase of the s -th component of the j -thsignal, respectively. Moreover, υ j,s is a local maximum point of the “tuning function”related to the fundamental integral of the frequency analysis method, i.e. υ
7→ T j ( υ ) = 1 T (cid:12)(cid:12)(cid:12)(cid:12) Z T d t [ Y j ( t ) + i X j ( t )] e − i υt W ( t ) (cid:12)(cid:12)(cid:12)(cid:12) . (47)Here, W is a suitable weight function such that R T d t W ( t ) = T . In all the computations,we have adopted a so-called “Hanning filter” adapted to the interval of time [0 , T ], i.e.,we have set W ( t ) = 1 + cos[ π (2 t/T − t Y j ( t ) + i X j ( t ), in practice we obviously limit ourselves to deal with a discretization ofthe signal made by sets of finite points computed on a regular grid in the interval [0 , T ],i.e., for t = i ∆, where i = 0 , . . . , N P and the sampling time is ∆ = T / N P .In order to produce the discretized signals (cid:8)(cid:0) Y j ( i ∆) , X j ( i ∆) (cid:1)(cid:9) N P i =0 to be analysed, wehave preliminarly used the symplectic integrator of type SBAB C , which is describedin [36]. The splitting of the initial Hamiltonian (1) in two integrable parts A and B (thatis requested for implementing such a kind of numerical integrator) is made so that thequadratic part is stored in A , i.e., A = P N − j =0 (cid:2) y j + ( x j +1 − x j ) ) (cid:3) , and B = H − A .In order to reduce the effects due to the accumulation of the round-off errors, we havefound convenient to perform the numerical integrations by implementing the long double floating-point arithmetics in our C code, while we used the standard double precisionvariables in all the C programs devoted to the FA.As a first application of FA, let us consider a very classical type of numerical explor-ation about the dynamical behavior of the FPU chains. Indeed, we are going to considerinitial conditions where just the first oscillation mode is excited, in agreement with whatwas done in the original experiments about the FPU model; therefore, we study thechaoticity of the motions starting from those initial data. Actually, we set x l (0) = r N A √ ν sin (cid:18) lπN (cid:19) , y l (0) = 0 , (48) ∀ l = 1 , . . . , N − E = H ( y (0) , x (0)) dependsuniquely on the amplitude A of the initial excitation of the first mode, since we keep α and β as fixed parameters (recall definition (1)). We numerically determine the globalmaximum point ω f ;1 of the function υ
7→ T ( υ ) defined in (47). Of course, the cal-culation of the integral in (47) is made by a quadrature method that is applied afterhaving replaced the function t (cid:0) Y ( t ) , X ( t ) (cid:1) (where t ∈ [0 , t ]) with the finite sequence (cid:8)(cid:0) Y ( i ∆) , X ( i ∆) (cid:1)(cid:9) T/ ∆ i =0 . Therefore, this same procedure is repeated so as to determine lliptic tori in FPU non-linear chains . . . | ∆ ω f; | E S = E/ N N=4N=8N=16 | ∆ ω f; | E S = E/ N N=4N=8N=16
Figure 1:
Variation of the first angular velocity as a function of the specific energy in FPUchains. All the considered motions start from initial conditions where just the first normal modeis excited. The results are plotted in red, green and blue when they refer to N = 4, 8 and 16,resp. On the left, the α − model is studied with α = 1 /
4; on the right, the β − model is analysedwith β = 1 / the global maximum point of the function defined in (47) when the integral appearing inthat formula is performed on the interval [ T , T ] instead of [0 , T ]. In Figure 1, the vari-ation | ∆ ω f ;1 | of the maximum points over those two subsequent time windows is reportedas a function of the specific energy E S = E/N for different number of particles, namely N = 4 , ,
16, for both a α –model and a β one (left and right boxes, resp.). All the signalsanalysed to produce the plots in Figure 1 are such that the sampling time is ∆ = 0 . T = 2 = 131072. Allthe preliminary numerical integrations (that have been necessary to create the signals)have been performed by setting the timestep at 0 . | ∆ ω f ;1 | are increasing with respect thenumber N of nodes: on one hand, it is known that the variation of the local maxima ofthe function (47) is O (1 /T ) when quasi-periodic motions are studied and the Hanningfilter is adopted (see [35]), on the other hand, ω f ;1 ≃ ν that is defined in (3), then thefundamental period corresponding to the main oscillation (of the first normal mode) is O ( N ). Therefore, in relative terms, when N is increased, larger times are needed to letthe angular velocities relax to their final values; this explains why the variation of | ∆ ω f ;1 | gets bigger with N , by keeping fixed the total time 2 T . By comparison with Figure 6of [7], we see that our results are in good agreement with those obtained by computingthe maximal Lyapunov exponent for the α -model: the transition to the instability occursfor values of the specific energy E S larger than 0 . .
02, when N = 8 and N = 16,respectively. However, the FA highlights also the weak chaoticity regions likely due tothe crossing of resonances. This occurs for values of the specific energy that are smaller8 C. Caracciolo, U. Locatelli than the previous ones and, then, such that the diffusion in the phase space is not veryremarkable. By comparing the plots in the right box of Figure 1 with those in the leftone, it is quite evident that both the crossing of resonances and the ultimate transitionto chaos is sharper in the case of the β -model. This was easy to expect, due to the factthat the perturbing terms are quartic instead of being cubic. Computational algorithms that are commonly used into the framework of FA can beefficiently adapted for the localization of elliptic tori. For such a purpose, we basicallyfollow the approach described in [43], where an attempt to provide a theoretical foundationof such a method is also made. The technique that iteratively refines the approximationof a generic n –dimensional elliptic torus is the keystone of this procedure. It can besummarized as follows, by referring to the specific case of the FPU chains, just to fix theideas. Before entering the details of the procedure, we emphasize that it is tailored forlower dimensional tori with a purely elliptic character in their transverse oscillations; themethod is expected to not be working in the case of hyperbolic dynamics in at least onedegree of freedom.(a) First, let us consider an initial condition ( Y (0) , X (0)) eventually close enough toan invariant elliptic torus.(b) By using a suitable numerical integrator, we produce N − (cid:8)(cid:0) Y j ( i ∆) , X j ( i ∆) (cid:1)(cid:9) T/ ∆ i =0 , N being the number of nodes, ∆ the sampling time, T thewidth of the total time interval of integration and j = 1 , . . . , N − Y j ( t ) + i X j ( t ) ≃ P N C s =1 A j,s e i( υ j,s t + φ j,s ) , N C being the fixed num-ber of summands. By eventually reordering the components, we can assume that theamplitudes are decreasing with respect to the index s , i.e., A j, ≥ A j, ≥ . . . ≥ A j, N C ∀ j = 1 , . . . , N − n –dimensionalvector ω f , whose components are the so called fundamental angular velocities.(d1) As in the previous subsection, we set ω f ;1 = υ , that is the global maximumpoint of the function T ( υ ) defined in (47).(d2) For each j = 2 , . . . , N −
1, we select the largest summand (corresponding to avalue, say, ¯ s j of the index s ) that cannot be explained as a linear combinationwith integer coefficients of the current components of the fundamental angularvelocities. This means that ∀ ≤ s ≤ ¯ s j − k ∈ Z j − suchthat (cid:12)(cid:12) υ j,s − P j − i =1 ¯ k i ω i (cid:12)(cid:12) ≤ ε tol with | ¯ k | ≤ K M , while (cid:12)(cid:12) υ j, ¯ s j − P j − i =1 k i ω i (cid:12)(cid:12) > ε tol ∀ k ∈ Z j − such that | k | ≤ K M . The parameters ε tol and K M are convenientlyfixed so as to let the whole procedure converge (see the discussion at the endof the present section). lliptic tori in FPU non-linear chains . . . j = 2 , . . . , N −
1, we extract the new j -th component of the funda-mental angular velocities vector from the value υ j, ¯ s j selected as explained atthe previous point (d2).(d4) When each j -th component of the fundamental angular velocities is expected tobe approximately equal to a value ν j that is known a priori (for instance, referto equation (3) in the case of FPU chains), then the choice can be adapted onsuch a purpose. In such a situation, first, ∀ j = 1 , . . . , N − k , . . . , ¯ k j such that (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ¯ k j υ j, ¯ s j + j − X i =1 ¯ k i ω f ; i − ν j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = min k ∈ Z j | k |≤ KM , kj =0 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k j υ j, ¯ s j + j − X i =1 k i ω f ; i − ν j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , therefore, we put ω f ; j = ¯ k j υ j, ¯ s j + P j − i =1 ¯ k i ω f ; i . In words, by linear combinationswith integer coefficents, we iteratively define a set of fundamental angularvelocities that are close as much as possible to the expected values.(e) We modify the initial conditions as follows. Let us introduce ( ˜ Y (0) , ˜ X (0)) so that˜ Y j (0) + i ˜ X j (0) = P N C s =1 ˜ A j,s e i φ j,s ∀ j = 1 , . . . , N −
1, where ˜ A j,s = 0 when (cid:12)(cid:12) υ j,s − k · ω f (cid:12)(cid:12) > ε tol ∀ k ∈ Z n such that | k | ≤ K M , otherwise we simply keep ˜ A j,s = A j,s .Once again, the previous formulæ can be rephrased in a more clear way as follows:the initial conditions are chosen in order to suppress all the components related toangular velocities that are not a linear combination (by integer coefficients) of thefundamental vector ω f , which has been determined as explained at the previouspoints (d1)–(d4).(f) We redefine the initial conditions so that ( Y (0) , X (0)) = ( ˜ Y (0) , ˜ X (0)). The exe-cution of the algorithm is positively concluded when the signals are close enoughto their quasi-periodic approximations based on the fundamental vector of angularvelocities ω f , i.e., we check if the following condition is satisfied ∀ j = 1 , . . . , N − { t i : t i = i ∆ , ∀ ≤ i ≤ T/ ∆ } (cid:12)(cid:12)(cid:12) Y j ( t i ) + i X j ( t i ) − P N C s =1 ˜ A j,s e i( υ j,s t i + φ j,s ) (cid:12)(cid:12)(cid:12) | A , | ≤ µ tol ;the execution of the algorithm is ended by setting a flag variable on an error statusif we reached a prefixed maximum value on the number of times the points (b)–(f)are iterated; otherwise, we restart from point (b).A good choice of the parameters ε tol , µ tol , K M and N C is crucial for avoiding infiniteloops. On one hand, by adopting a small value of the tolerance errors ε tol and µ tol , thealgorithm is forced to look for very accurate quasi-periodic approximations. On the otherhand, there are so many calculations involved (during both the numerical integrations To fix the ideas, in our numerical experiments such an upper bound has been set so as to be equalto 20. C. Caracciolo, U. Locatelli and the FA) that the accumulation of the round-off errors is relevant; therefore, if weask for too much precision, for instance, the stop condition described at point (f) abovecould never be fulfilled. We are obviously interested in keeping a large enough valueof the upper bound K M on the l –norm of the Fourier harmonics k . However, becauseof the Fourier decay of the coefficients in the Hamiltonian expansions, terms related toharmonics having a big norm | k | are expected to play a negligible role; therefore, it isessential to adopt a small enough value for the upper bound K M in order to prevent thewrong choice of spurious terms in the definitions introduced at the previous points (d2),(d4) and (e). The pros and cons concerning the choice of the number of components aresimilar to those regarding the tolerance errors: small values of N C decrease the accuracyof the selected approximations at the end of the whole procedure, but too large valuesof N C increase dramatically the computational time and could prevent its convergencebecause of the fact that components related to very small amplitudes are strongly affectedby numerical errors.The applications of this method, based on FA and designed for the localization ofelliptic tori, are discussed in the next section. For all those numerical computations, wehave found that it is convenient to set ε tol = 10 − , µ tol = 2 × − , K M = 20 and N C = 25. Moreover, we will keep fixed the parameters concerning also the preliminarysymplectic integration (that are the timestep, ∆ and T ), whose values are reported in thecomments to Figure 1 in the previous subsection 4.1. The approach described in section 3 can be implemented in a programming code. Forsuch a purpose, we have used X ̺ ´o νoς , a software package especially designed for doingcomputer algebra manipulations into the framework of Hamiltonian perturbation theory(see [28] for an introduction to its main concepts). In this section we will discuss the resultsof the application of the algorithm constructing elliptic tori for the FPU chains; moreover,we will compare these results with the technique for doing numerical explorations thathas been previously described in section 4. Let us start by studying a non-trivial case that is as simple as possible. We focus ona FPU chain having N + 1 particles, with N = 4. This means that the system has N − N = 2 m with m ∈ N . First, we begin studying 2D elliptic tori.As it has been widely stressed in section 2, all the procedure constructing the normalform for a n –dimensional elliptic torus basically depends on the preliminar translationvector I ⋆ ∈ R n + , which also rules the size of the perturbation. Moreover, looking at the Since the pioneering work [15], it is usual to set N equal to an integer power of 2. This allows tocompute very efficiently the passage from the original coordinates ( y , x ) to those related to the normalmodes, that are ( Y , X ), and vice versa (see formula (2)). lliptic tori in FPU non-linear chains . . . Loga r i t h m o f t he no r m Normalization step: r 1e-251e-201e-151e-101e-051e+000 2 4 6 8 10 12
Loga r i t h m o f t he no r m Normalization step: r
Figure 2:
Numerical check of the convergence of the algorithm constructing the normal formrelated to an elliptic torus for the FPU model: we plot the norm of the generating functions ateach normalization step r . The symbols +, × , △ , ∗ and (cid:3) refer to the norm of the generatingfunctions χ ( r )0 , χ ( r )1 , X ( r )2 , Y ( r )2 and the diagonalizing transformation D ( r ) . The norm is calcu-lated by adding up the absolute value of all the coefficients appearing in the expansion of eachgenerating function or the linear change of coordinates defined by D ( r ) , for what concerns thelatter case. description of the algorithm made in section 3, it is expected to converge if and only ifthe generating functions decrease in a geometrical way. Actually, an automatic criterionto decide about the practical usefulness of the normal form expansions can be given bymonitoring the behavior of the norms of the generating functions (see, e.g., [51]). In orderto fix the ideas, it is convenient to see how such an approach can be adapted to the presentframework.The results summarized in Figure 2 refer to the α − model with α = 0 .
25, where thecomponents of the initial translation ( I ⋆ , I ⋆ ) are fixed as (0 . , . ,
1) forwhat concerns the left box and the right one, respectively. Actually, the construction ofthe normal forms related to both those translation vectors are reported also in Figure 3:they correspond to the symbols plotted in the centre of the left box and in its top-right corner, resp. For the sake of simplicity, we have decided to truncate the initialexpansions of H (0) in (10) so that the total degree ℓ of the classes of functions P sKℓ isnot greater than 8. Such an upper bound is automatically satisfied by any Hamiltonian H ( r ) , because the Poisson brackets with the generating functions do not increase thetotal degree in the square root of the actions. In order to deal with all the computationsrequired by the algorithm in a reasonably short amount of time, we have limited ourselvesto perform just the first 12 normalization steps. With these parameters, the computationof a single 2D elliptic torus requires 77 seconds of CPU time on a computer equippedwith a Xeon 18-Core 6154 (3.0 Ghz) processor.
Accordingly, we have decided totruncate the Fourier expansions so as to contain all the terms of the final generatingfunctions X (12)2 , Y (12)2 ∈ P K . With these last settings, the computational algorithmdescribed in section 3 is completely defined. In the semi-log plot reported in the left box Let us recall that in section 3 we decided to set K = 2. C. Caracciolo, U. Locatelli E E N=4, α = 0.25, β = 0 E / E E N=4, α = 0.25, β = 0 Figure 3:
On the left, success/failure of our procedure constructing 2–dimensional elliptic torifor different translation vectors, the symbol is plotted when the algorithm is considered tobe convergent, because it satisfies the rules (A)–(B), the △ otherwise. On the right, plot ofthe ratio between the energy of the first two normal modes (see the text for a more precisedefinition) for motions of the same type as those studied in the first work about FPU chains.The curve appearing in the right box has been reported also on the left as a solid line. of Figure 2, the decrease of the norms of χ ( r )0 , χ ( r )1 , Y ( r )2 , X ( r )2 is pleasantly regular apart afew unavoidable fluctuations. Such generating functions are listed according to their sizein an increasing order, that is coherent with the expected accumulation of small divisors,as it has been described in [26] and in the reference reported in footnote 5 . Since thelinear canonical transformation D ( r ) is determined with the aim to remove the same kindof terms for any normalization step r , there is not any small divisor entering the solutionof equation (36). Therefore, it is natural to expect that the norm corresponding to D ( r ) is smaller than any other. This explains also why its plot in the left box of Figure 2is the most irregular, since it is very exposed to the numerical effects due to round-offerrors, because of its smallness. In short, the generating functions introduced at the 12-th normalization step are 6 or 7 order of magnitude smaller than those defined at thebeginning of the algorithm and their decrease is rather regular. This strongly suggeststhat the sequence of the corresponding canonical transformations is more and more closeto the identity and their composition is converging to the change of coordinates bringingthe Hamiltonian to the wanted normal form. On the other hand, the behavior of the normsappearing in the right box of Figure 2 is not converging at all. Looking at the differentvalues of the vectors I ⋆ that are corresponding to the left and right boxes, respectively,and keeping in mind the discussion in section 2, the following explanation of the observedphenomena looks quite natural: in the case of the right box, the initial translation is solarge that the induced perturbation prevent the existence of the wanted elliptic torus.This does not occur when I ⋆ and the corresponding perturbation are small enough, as forthe plot in the left box.We are now ready to formulate some practical rules, in order to establish when we canassume that the procedure is properly working or not. We consider that the normalizationalgorithm is converging if both the following conditions are satisfied: lliptic tori in FPU non-linear chains . . . k X ( r )2 kk X (1)2 k + k X (2)2 k < (0 . r − for r > k X (¯ r )2 kk X (1)2 k < , where ¯ r = 12 is the last normalization step performed.Let us stress that the geometrical decrease of the generating functions (when it occurs)is usually so sharp that a reasonable modification of the values appearing in the r.h.s.of the inequalities included in the rules (A)–(B) above does not considerably affect theresults. We have decided to refer the tests to the sequence X ( r )2 because their normsare usually greater than the corresponding ones for all other generating functions, asalready remarked in the comments about Figure 2. Let us recall that the functionalnorm is given by the sum of the absolute values of all the coefficients appearing in the(finite) expansion. We have considered many different initial translations, so that thecomponents of every vector I ⋆ are located in correspondence with the points of a regular2D grid of log–log type. For each initial translation I ⋆ making part of such a finite set,we have performed the normalization procedure, at the end of which we have applied thepractical criterion described above so as to discriminate between the convergence and thedivergence of that single application of the constructive algorithm. The successes andfailures of this explicit method are summarized in the left box of Figure 3. Actually, thelocation of the points making part of the regular grid is rescaled in terms of the energyof the quadratic approximation based on the normal modes oscillations; this means thatwe put the values of ν I ⋆ and ν I ⋆ in abscissas and ordinates, resp. This has been donefor an easier comparison with the results reported in the right box of the same figure.The plot in the left box of Figure 3 has been realized by considering values of theenergy of the first two normal modes in a wide range from 10 − to 5; the normalizationalgorithm looks to be convergent often, even when both those values are approximatelyequal to 1. This is certainly a remarkably good performance, because it is not far fromthe upper bound on the specific energy over the which quasi-periodic motions are no moredetectable, by exploring the dynamics with the FA method (see Figure 1). However, amore careful analysis of the results highlights that the motions studied in the presentsection are quite different with respect to those of the previous one. Let us now focuson the same initial conditions of “semi-sinusoidal” type (48) that have been consideredin [15]; then, we can numerically integrate the equations of motion and calculate thequasi-periodic decomposition (46). In the right box of Figure 3, we report the ratio( ν , A , ) / ( ν , A , ) of the energies of the maximal quasiperiodic components related tothe second and the first normal mode of oscillation, respectively. The plot of such aquantity, as a function of the total energy, clearly shows that the oscillations with angularvelocities close to that of the second normal mode get considerably excited just for largeenough initial amplitudes of the first mode. On the other hand, the more the energy islow, the more the sharing between normal modes is inhibited (as it is well known, see,e.g., [4]). For values of the total energy that are relatively small, it is reasonable to assumethat the two different ways we adopted to define the energies related to the normal modesare approximately equal. Therefore, it makes sense to draw the same curve that appearsin the right box of Figure 3 also on the left, after having reported the numerical data fromthe former to the latter in the corresponding units. Such a graphical operation makes4 C. Caracciolo, U. Locatelli evident that the motions generated by initial conditions of “semi-sinusoidal” type are notin the domain of applicability of our algorithm when it is used to explicitly construct 2Delliptic tori. Furthermore, the curve we have drawn in the left box looks parallel andwell separated with respect to a sort of natural border of the domain of applicability ofthe method. This remark and the inhibition effect on the sharing of the energy (that isinitially exciting just the first mode) strongly suggest that the kind of motions originallystudied in the FPU chains are substantially more similar to those on 1D elliptic tori.This is the reason why, hereafter, such invariant objects will play a prominent role in ourinvestigations, in spite of the fact that they are simply linearly stable periodic orbits.
One-dimensional elliptic tori in FPU chains can be found, by using a continuation methodbased on the algorithm explained in subsection 4.2. Let us recall that in the scientificliterature one can find several different numerical techniques that have been shown tobe very efficient in searching for periodic orbits. However, we prefer a method basedon FA that is designed to work properly just for (quasi)periodic motions lying on lowerdimensional tori having an elliptic character in the transverse degrees of freedom. Thisis perfectly coherent with the type of normal forms we aim to construct. In this con-text, let us also mention that there are also methods that are mainly analytic by theirnature. For instance, an approach based on the Reeb’s theorem requires to perform apreliminary averaging procedure on the angular part of the canonical coordinates; then,non-degenerate equilibrium points of the reduced (and averaged) system are proved tocorrespond to periodic solutions of the original one. Moreover, they also inherit theireventual linear stability from the corresponding equilibrium points (see [46] and [34] fora recent application to rotating H´enon–Heiles systems).In order to describe in a more detailed way the continuation method based on FA, theprocedure can be summarized as it follows. • The numerical procedure of subsection 4.2 is performed from point (a) to (f), startingfrom an initial condition ( Y (0) , X (0)) of type (48), where A is so small that suchan algorithm is able to easily find the corresponding 1D elliptic torus. At the endof this stage, we compute the value of the energy ˜ E = H ( ˜ Y (0) , ˜ X (0)), where theHamiltonian H is defined in (4) and ( ˜ Y (0) , ˜ X (0)) is defined at the last executionof point (e). • Every time the conclusion at point (f) is positive, the corresponding value of theenergy ˜ E is updated in the same way we described just above. Therefore, the wholeprocedure of subsection 4.2 is restarted from point (a), right after having redefined( Y (0) , X (0)) = (1 + ζ )( ˜ Y (0) , ˜ X (0)), where ζ is a small parameter (see below) and( ˜ Y (0) , ˜ X (0)) is the initial conditions vector corresponding to the last computationof ˜ E , which is also the last numerical determination of the elliptic torus. • If the prescription at point (f) returns a negative conclusion and ζ is still larger thana prefixed minimum value ¯ ζ , then ζ is reduced by half. Therefore, the procedure of lliptic tori in FPU non-linear chains . . . • If the small parameter ζ ≤ ¯ ζ , then the whole continuation method is finally stopped. ω E S = E / N N=4, α = 0.25, β = 0 ω E S = E / N N=4, α = 0, β = 0.25 ω E S = E / N N=8, α = 0.25, β = 0 ω E S = E / N N=8, α = 0, β = 0.25 Figure 4:
Comparisons between 1D tori constructed by our normal form algorithm (dots denotedby • ) and those detected by using the continuation method (solid lines), for both the α − modeland the β one (on the left / right, resp.) in the cases with N = 4 , The solid lines appearing in Figure 4 are obtained by applying this continuationmethod. In such numerical experiments, the value of ζ has been initially set at 0 . ζ = 0 . E S instead of the total one. By comparing Figures 4and 1, one can appreciate that for the β –model the breakdown threshold with respectto E s for the invariant manifold corresponding to 1D elliptic tori is approximately thesame as that for the motions starting from the same type of initial conditions origin-ally considered in [15] and for which the transition to chaoticity looks very sharp. Forwhat concerns the α –model, 1D elliptic tori appear to be more robust also because theyseem not being affected from the transitions through resonances, that apparently exertsa strong impact in the numerical experiments summarized in Figure 1.6 C. Caracciolo, U. Locatelli -3.0e+00-2.0e+00-1.0e+000.0e+001.0e+002.0e+003.0e+000.00001 0.00010 0.00100 0.01000 0.10000 1.00000 θ E S = E / N | ∆ ω f; | E S = E/ N Figure 5:
Study of the transition to chaoticity in the β –model of the FPU chain with N = 4(and β = 1 / λ ± j = e ± i θ j related to any monodromymatrix, that has been computed for each of the 1D elliptic tori considered in the top-right boxof Figure 4. More precisely, the values of ± θ j for j = 1 , , In Figure 4, there are also some points that have been plotted with bold dots, so thatthey can be easily seen. They refer to the applications of our normalization algorithm,which has been described in section 3 and has been adapted to construct 1D elliptic tori.Actually, we have selected a finite set of values for the initial translation I ⋆ ∈ R + , insuch a way that the sequence log I ⋆ makes a regular grid. For each single element of sucha set, we have performed our whole normalization procedure starting from that value ofthe initial translation I ⋆ . Let us now focus on the comparison between expansions (11)and (15). When the algorithm seems to be converging from a numerical point of view,the limit values ω ( ∞ )1 and E ( ∞ ) , which are quantities characterizing the wanted 1D elliptictorus, are expected to be well approximated by ω (¯ r )1 and E (¯ r ) , resp., which are explicitlycomputed during each construction of the normal form up to the step ¯ r . Therefore, alsoour semi-analytic algorithm provides the value of both the angular velocity and the energylevel which are directly comparable with those given by the continuation method basedon FA, that has been previously described in this same section. Looking at Figure 4, onecan appreciate that the agreement between these two computational methods is excellent,because the superposition between the plots is nearly perfect in all the four cases we haveconsidered.The extension of the range for relatively large values of the energy is definitely thehardest challenge for what concerns the applicability of our algorithm. However, let uscheck that the numerical method we used for locating elliptic tori does not introduceany artificial simplification. In other words, we want to ensure that the family of periodiccurves cannot be easily prolonged by continuation over what we call breakdown threshold.This is because of a very concrete dynamical obstruction: the transition to a resonant lliptic tori in FPU non-linear chains . . . β –model with N = 4. Letus remark that the monodromy matrix related to a linearly stable periodic orbit can beeasily calculated (see, e.g., [30] for a gentle introduction to these concepts). This canbe done by referring to the corresponding normal form approximation of a 1D elliptictorus, that we can write as ω P + P j =1 Ω j (Ξ j + Θ j ) /
2, in analogy with the statementof theorem 3.1. In fact, its six eigenvalues are such that λ ± j = exp( ± i2 π Ω j /ω ) for j = 1 , λ ± = 1, being the period equal to 2 π/ω . For what concerns the graphsreported on the left of Figure 5, let us emphasize that we focus on all the 1D elliptictori that have been numerically found, in order to draw the solid curve appearing inthe top-right box of Figure 4. One can easily deduce that the ratio Ω /ω is going toconverge to 2 in correspondence with the breakdown threshold, whose characteristicvalue can be expressed in terms of the specific energy as E S ≃ .
25. In the box onthe right of Figure 5 we consider a slight modification of every initial condition relatedto a 1D elliptic torus, whose the corresponding monodromy matrix has been studied inthe plot on the left. In fact, the values of the canonical coordinates are multiplied by acommon factor 1 .
01 for all these initial conditions; then, the corresponding variation of thefirst angular velocity | ∆ ω f ;1 | is computed (according to the description in subsection 4.1)and it is plotted in the graph included in the box on the right. The sudden jump byeight orders of magnitude that is corresponding to the last evaluation about the variationof the fundamental frequency, jointly with the comments concerning the eigenvalues ofthe monodromy matrix, allows us to draw the following conclusion. 1D elliptic tori areentering the chaotic zone surrounding the 2:1 resonance between the angular velocities ω and Ω , when the specific energy is converging to the breakdown value E s ≃ .
25. Thisexplains why the one-parameter family of 1D elliptic tori cannot be further prolonged forlarger values of E S , by preserving its continuity. The same situation has been observedin the cases with the α –model instead of the β –one and/or N = 8, where the patternsof the graphs representing the eigenvalues of the monodromy matrix are eventually morecomplicated. For the sake of brevity, the figures corresponding to these other cases areomitted.Looking at the comparison between the boxes on top of Figure 4 with those on bottom,it is evident that our algorithm works better in the case N = 4 with respect to N = 8. Ingeneral, it is very natural to expect that, for our semi-analytic procedure, the situationis more difficult when the number of degrees of freedom is increased. In particular, wehad to reduce a few parameters ruling the truncations of the Hamiltonians, in orderto make reasonable both the computational time and the usage of memory in the casewith N = 8. In fact, while for N = 4 these parameters are exactly the same as thosedescribed in subsection 5.1, when N = 8 we have decided to truncate the expansions of H ( r ) so that the total degree in the square root of the actions is not greater than 4, for Let us consider the dashed line approaching θ ± = 0 from below. Looking at such a plot, it isevident that the ratio θ j / (2 π ) = Ω j /ω → m , for some integer numbers j and m . By following sucha curve backward for energies smaller and smaller, we can see that it is going to converge to the value2 πν /ν − π ≃ −
1, where ν and ν are given in (3). They are the angular velocities of the first andsecond normal mode in the case of small oscillations, so they are equal to the limit values of ω and Ω for E S →
0, respectively. Therefore, we can conclude that j = 1 and m = 2. C. Caracciolo, U. Locatelli any normalization step r ranging from 0 to ¯ r , where this upper bound on r has beendecreased to ¯ r = 8 (in the case with N = 4, it was fixed to ¯ r = 12). Since we havekept the same practical criterion to decide if the algorithm is converging or not, thismeans that rule (B) of subsection 5.1 has become more restrictive, because the decreaseof the generating function must be faster to compensate the reduction on the numberof normalization steps. For what concerns the CPU usage, the computation of a single1D elliptic torus requires around 33 seconds when N = 4 and 112 seconds in the case with N = 8. Of course, such a comparison is rather rough, because the truncation parametershave been chosen in a different way with the aim of accomplishing all the computationaltasks in a reasonable total time, as we have discussed above. Both the truncation rulesand the computational resources are exactly the same as in the case of the 2D elliptictori that has been described in the previous subsection. Here, the time needed by thecomputation of a single invariant torus is about 43 % with respect to the amount that hasbeen reported during the discussion of Figure 2. The reason of such a gain can be easilyunderstood by looking at the expansion of a generating function; to fix the ideas, let usfocus on equation (26). For instance, since the number of summands making part of thegenerating function χ ( r )1 is increased by a factor r passing from n = 1 and n = 2 to n = 2 and n = 1 (let us recall that we fixed K = 2 when dealing with FPU models), oneimmediately realizes that the whole normalization procedure is based on Poisson bracketsthat involve many more terms in the case of the explicit construction of 2D elliptic tori.This explains why more computational time is needed in that case.Apparently, our approach works better when it is applied to the β –models instead of α ones; this is highlighted by the comparison between the boxes on the left of Figure 4with those on the right. In our opinion, this phenomenon is due to two main reasons.First, the transition to the chaoticity is much sharper in the case of a quartic perturbationinstead of a cubic one; therefore, it is natural to expect that a normal form method canapproach better the breakdown threshold for any invariant manifold, when the β –model isconsidered. Moreover, there is an accurate integrable approximation of the α –model: theToda lattice, because their expansions coincide up to terms of third order (see [18] and,e.g., the introduction of [47] in [20]). We think that this fact should play an important rolein order to stabilize the dynamics when relatively high energies are considered, but ourapproach is completely blind to such an accurate approximation provided by a Toda–likeHamiltonian. This further explains why we have achieved better results in the case of the β –model.By comparing the cases reported in Figure 4 all together, the best performance is in the β –model with N = 4. In that case, our algorithm is able to construct the normal formfor a 1D elliptic torus, when the energy is about 93 % of the value corresponding toits breakdown threshold E S ≃ .
25, as it has been numerically computed by using thecontinuation method based on FA. The results for the β –model are still good for N = 8,where our algorithm fulfills the construction of tori up to 64 % of the threshold energy.On the other hand, for the α –model the performance is good only for N = 4, since for N = 8 we are really far from the numerical threshold. lliptic tori in FPU non-linear chains . . . According to the discussion at the end of subsection 5.1, the orbits generated by initialconditions of “semi-sinusoidal” type are close to be 1D elliptic tori. In the present section,we enforce this concept, trying to explain the apparent stability (i.e., quasi-periodicity)of the former motions in terms of their vicinity to the latter ones. We emphasize that1D elliptic tori are expected to be quite robust, because either the construction of the cor-responding normal form is stopped (due to the eventual occurrence of resonant terms) or,if it is converging, arbitrarily small divisors cannot be produced. In order to understandthis argument, it is convenient to reconsider the non-resonance condition (42), in the casethe sequence of the angular velocity vector ω ( r − is one-dimensional. Here, we aim todescribe how the elliptic tori influence the dynamics in their neighborhood, by imple-menting a suitable Birkhoff normalization algorithm. Thus, we basically follow the ideasdescribed in [39] with an approach similar to that fully developed in [25] and [27], wherethe main perturbative terms are removed from a Hamiltonian in Kolmogorov’s normalform by an iterative procedure. Moreover, for those dynamical systems such a methodhas been also complemented with suitable final estimates about the time of effective sta-bility. We adapt that approach in such a way to remove the terms eventually dependingon the angles which appear in the expansion of the normal form Hamiltonian (11) andare o ( k P k + k ( Ξ , Θ ) k ).Let us first rewrite the expansion (11) of Hamiltonian H ( ∞ ) , that is the normal formcorresponding to an elliptic torus, in the following way: H (0) ( P , Q , J , ϕ ) = n X j =1 ω ( ∞ ) j P j + n X j =1 Ω ( ∞ ) j J j + X ℓ> f (0) ℓ ( P , Q , J , ϕ ) , (49)where ( J , ϕ ) are the action–angle variables related to the canonical coordinates ( Ξ , Θ )that are transverse to the elliptic torus; moreover, the perturbing terms f (0) ℓ are of totalorder ℓ in the square root of the actions ( P , J ), while they are of trigonometrical degreenot greater than ℓ with respect to the angles ϕ . Indeed, such functions f (0) ℓ are given bythe series P s ≥ f ( ∞ , s ) ℓ , where the summands f ( ∞ , s ) ℓ appear in (11). Therefore, the angulardependence of f (0) ℓ is analytic with respect to Q .By performing a finite sequence of canonical transformations (expressed, once again,in terms of the Lie series formalism), we aim to lead the Hamiltonian to the Birkhoffnormal form up to order r , that is H ( r ) ( P , Q , J , ϕ ) = r X s =0 Z s ( P , J ) + X ℓ>r f ( r ) ℓ ( P , Q , J , ϕ ) , (50)where the first part Z ( r ) = P rs =0 Z s is the integrable part, because it depends on theactions only, while R ( r ) = P ℓ>r f ( r ) ℓ = O (cid:0) k ( P , J ) k ( r +3) / (cid:1) is the remainder.Let us briefly describe the r -th normalization step in a fully generic way. For such apurpose, it is convenient to recollect the action–angle variables in ( I , ϑ ) ∈ R n × T n , where I = ( J , P ) and ϑ = ( Q , ϕ ). We focus on the expansion of the Hamiltonian produced by0 C. Caracciolo, U. Locatelli the first r − H ( r − ( I , ϑ ) = P r − s =0 Z s ( I ) + P ℓ ≥ r f ( r ) ℓ ( I , ϑ ),where the main integrable term Z can be rewritten as ν · I with ν = ( ω ( ∞ ) , Ω ( ∞ ) ). Weaim to remove those perturbative terms having a total degree in the square root of theactions that is equal to r + 2. Therefore, the r -th generating function is determined asthe solution of the homological equation (cid:8) χ ( r ) , ν · I (cid:9) + f ( r − r = Z r . (51)This means that the new summand making part of the normal form is given by Z r ( I ) = h f ( r − r i ϑ , while the new generating function can be written as χ ( r ) ( I , ϑ ) = X | m | = r +2 X k =( k , k ∈ Z n × Z n k = , | k |≤ r +2 I m / (cid:20) − c k , m sin( k · ϑ ) k · ν + d k , m cos( k · ϑ ) k · ν (cid:21) , (52)where the coefficients c k , m and d k , m appear in the following expansion of the perturbativeterm to be averaged: f ( r − r ( I , ϑ ) = X | m | = r +2 X k =( k , k ∈ Z n × Z n | k |≤ r +2 I m / (cid:2) c k , m cos( k · ϑ ) + d k , m sin( k · ϑ ) (cid:3) . (53)As it is usual, the iterative definition of the functions f ( r ) ℓ ∀ ℓ > r appearing in formula (50)can be easily obtained. In fact, it is just matter to gather separately all the summandshaving a total degree in the square root of the actions equal to ℓ +2 among those appearingin the expansion of the Lie series which defines the Hamiltonian in Birkhoff normal formup to order r , i.e., H ( r ) = exp (cid:0) L χ ( r ) (cid:1) H ( r − . (54)After having replaced r with r + 1, the next normalization step could be exhaustivelydefined by using again all the formulæ (51)–(54). This ends the description of the iterativealgorithm constructing the Birkhoff normal form up to a finite order in the neighborhoodof an elliptic torus.In order to ensure that the homological equation (51) can be solved by the generatingfunction χ ( r ) written in (52), once again we have to assume a suitable non-resonancecondition on the angular velocity vector, for instance, the Diophantine one, i.e.,min k =( k , k ∈ Z n × Z n k = , | k |≤ r +2 | k · ν | ≥ γ ( r + 2) τ , (55)for some fixed value of the parameters γ > τ ≥ n −
1, with n = n + n . Let usrecall that the algorithm constructing the Birkhoff normal form is not convergent. Thismeans that we cannot make the Hamiltonian integrable, by iterating the procedure up toinfinity. Indeed, it is convenient to perform the r -th normalization step until the sup–normof the remainder R ( r ) = P ℓ>r f ( r ) ℓ is decreasing (with respect to the previous one, thatis R ( r − ) on the set { ( I , ϑ ) ∈ B ̺ ( ) × T n } , where ̺ is the maximal distance in actionfrom the elliptic torus in its neighborhood. The so called estimates `a la Nekhoroshev lliptic tori in FPU non-linear chains . . . r ∝ ̺ − / [2(1+ τ )] ;therefore, it is easy to prove that sup B ̺ ( ) × T n |R (¯ r ) | is bounded from top by a quantity ∼ exp (cid:0) − ( ̺ ∗ /̺ ) / [2(1+ τ )] (cid:1) that is exponentially small with respect to the inverse of afractional power of the distance ̺ from the torus, ̺ ∗ being a suitable positive constant.In words, this phenomenon is due to the fact that the distance ̺ from the torus rules theperturbation: in its vicinity, the integrable approximation is better and the perturbation(that is given by the remainder) is smaller; otherwise, the approximation provided by theBirkhoff normal form is not so accurate, the perturbative terms play a remarkable roleand just a few normalization steps can be profitably performed. Much more details canbe found in the introductory lectures [22], while a quantitative analysis in the context ofsome computer-assisted proofs is described in [5].Birkhoff normal forms are commonly used in order to provide lower bounds on thediffusion time that is needed by a Hamiltonian system to escape from a local region of thephase space. For the sake of simplicity, let us rewrite the Birkhoff’s normal form up toorder r , that is described in (50), in the following more concise way: H ( r ) ( I , ϑ ) = Z ( r ) ( I )+ R ( r ) ( I , ϑ ). Since ˙ I = (cid:8) I , H ( r ) (cid:9) = (cid:8) I , R ( r ) (cid:9) , the smaller is the remainder, the longerwill be the escaping time from the neighborhood B ̺ ( ) × T n (around the elliptic torus)for any motion starting from an initial condition with k I (0) k ≤ ̺ < ̺ . One can easilyverify that also the lower bound for the diffusion time is of order ∼ exp (cid:0) ( ̺ ∗ /̺ ) / [2(1+ τ )] (cid:1) for ̺ → + . Therefore, for some realistic models of physical interest, this can be usedto ensure the effective stability, because the diffusion time can be made larger than theexpected lifetime of the system, by choosing a value of ̺ that is small enough. In thepresent work, we do not provide any evaluation of the drift speed ˙ I , also because thereis not a natural lifetime of the FPU model to be compared with the estimates about thediffusion time. Therefore, we prefer to show the evidence that the quasi-periodic motions,observed in the numerical experiments reported in subsection 4.1, are generated by someinitial conditions, which are so close to an elliptic torus that the remainder of the Birkhoffnormal form can be made effectively small.In order to produce the plots reported in Fig. 6, we reconsider all the normal formsrelated to the one-dimensional elliptic tori we have explicitly constructed for both the α − model and the β one in the cases with N = 4 ,
8. More precisely, we study separatelyevery Kolmogorov-like normal form of type (11), which corresponds to a dot (denotedby • ) that has been plotted in Fig. 4. Let us recall that we have preliminarly producedsuch a Hamiltonian by several algebraic manipulations (using the libraries of utilitiesprovided by the software package X ̺ ´o νoς ). In order to not increase dramatically thecomputational complexity of the normalization algorithm, we prefer to avoid the searchfor the optimal step. Indeed, for each of those Kolmogorov-like normal forms, we performa further Birkhoff normalization up to order r = 5 for N = 4 and r = 1 for N = 8.Moreover, ∀ s = 0 , . . . , r we truncate the Hamiltonians H ( s ) , defined by the iterativeprocedure summarized by formulæ (51)–(54), up to the maximal trigonometric degree24 (or 16) in the angles ϑ and up to order 8 (or 4) in the square root of the actions I when N = 4 ( N = 8, respectively). Let us stress that these limits force us to manipulateseveral millions of terms in the case with N = 8, because the Hamiltonians depend on2 C. Caracciolo, U. Locatelli seven pairs of canonically conjugate variables. For each of the Birkhoff normal forms H ( r ) ( I , ϑ ) = Z ( r ) ( I ) + R ( r ) ( I , ϑ ) expanded in the vicinity of a 1D elliptic torus, we selectamong the initial conditions of “semi-sinusoidal” type the configuration minimizing theabsolute value of the remainder |R ( r ) | . In view of the truncation parameters, such aremainder R ( r ) = P ℓ>r f ( r ) ℓ is of order 8 (or 4) with respect to k I k / . In Fig. 6, for eachof the Birkhoff normal forms H ( r ) we have computed, we plot the value of the specificenergy for such an initial condition and the corresponding minimum of the absolute valuefor what concerns the remainder |R ( r ) | .It is very meaningful to compare every case reported in Fig. 6 with the correspondingone in Fig. 1. The simpler behavior is shown in the case of the β − model with N = 4: withthe exception of some fluctuations for very small amplitudes of the initial conditions of“semi-sinusoidal” type that are probably due to round-off errors, the size of the remaindergrows linearly as a function of the specific energy in the semi-log plot (i.e., with theprescribed polynomial law) up to values of E S ≃
1. This is in agreement with thecorresponding plot of the variation of the first angular velocity in Fig. 1, that is of orderof the roundoff error until a sharp transition to larger values (that we associate to a chaoticregime) occurs for E S &
1. In the case of the α − model with N = 4, the chaotic regimecan be appreciated in Fig. 1 in the region where E S > . E S that is slightly larger than 1. Moreover, the slope of the curve looks lower inthe region with E S ∈ (0 . ,
1) with respect to the trend we can appreciate for smallerenergies. Let us recall that such a change of the slope in the behavior of the remindersize is usually connected with a decrease of the optimal step, when the normalizationalgorithm `a la
Birkhoff is performed (see, e.g., the discussion in subsection 4.3 of [27]).This explains why such a procedure is less effective when E S ∈ (0 . , α − model is considered with N = 8, in Fig. 1 a few of very narrow chaotic regions canbe observed in the region with E S ∈ (0 . , . E S is approximately in the interval (0 . , .
2) and the graph stopsthere. In the last case, that is the β − model with N = 8, the variation of the first angularvelocity in Fig. 1 shows a sudden increase in the interval of values of E S ∈ (0 . , . E S & .
2. There is a smalltrack of the intermediate chaotic region in the behavior of the size of the remainder inthe bottom-right box of Fig. 6 (see the rather irregular plot of three consecutive pointswith specific energy E S . . E S ≃ .
2. In short, the From a conceptual point of view, it would be more elegant to fix the initial conditions in advance,in such a way to look for the Birkhoff normal form minimizing the remainder. We stress that such aprocedure would require to handle many huge expansions to be conveniently stored at the same time.Instead, the computational algorithm described in the text above (that prescribes to consider a Birkhoffnormal form at once) allows us to enormously reduce the requirements about the memory storage. lliptic tori in FPU non-linear chains . . . -5 m i n | R ( ) | E S = E / N N=4, α = 0.25, β = 0. m i n | R ( ) | E S = E / N N=4, α = 0., β = 0.25 -5 m i n | R ( ) | E S = E / N N=8, α = 0.25, β = 0. -5 m i n | R ( ) | E S = E / N N=8, α = 0., β = 0.25 Figure 6:
Energy level |R ( r ) | for the remainder of the Birkhoff normal form around an elliptictorus (after having performed r = 5 , N = 4 ,
8, respect-ively). Such an evaluation is done for the initial conditions of “semi-sinusoidal” type minimizingthe remainder itself. The corresponding value of the specific energy for such initial conditionsis reported on the axis of abscissas. The α -model and the β one are considered on the left andon the right, resp.; the cases with N = 4 , • ) are connectedby dashed segments. All the other graphs refer to 1D elliptic tori. coherence between the data reported in Figs. 6 and 1 reinforces the idea that a motionstarting from initial conditions of “semi-sinusoidal” type is quasi-periodic, when it is in thestability region surrounding an elliptic torus, where the size of the remainder is negligiblewith respect to the value of the total energy of the Hamiltonian.In the top-left box of Fig. 6, the size of the remainder for the Birkhoff normal formsaround 2D elliptic tori is also included (the corresponding plot is connected by dashedsegments). In more detail, this further process performing the elimination of the per-turbative terms starts from the normal forms represented by the dots appearing on thediagonal edge of the region marked with many symbols • , that are plotted in the left boxof Fig. 3. Let us recall that they are expected to correspond to the 2D elliptic tori thatare the closest ones with respect to the initial conditions of “semi-sinusoidal” type. By4 C. Caracciolo, U. Locatelli comparing those plots, one can appreciate that the remainders of Birkhoff normal formscentered about 2D elliptic tori are much larger than the corresponding ones that are re-lated to linearly stable periodic orbits. This fact unambiguously shows that the latterare much more effective, for what concerns the approximation of the orbits originating bythe initial conditions of “semi-sinusoidal” type. Such a remark is in agreement with thediscussion at the end of subsection 5.1 and it can be seen as a natural completion of that.
In this work we have shown the effectiveness of our procedure explicitly constructing thenormal form for elliptic tori. It is based on an adaptation of the algorithm introducedin [49] and proved to be convergent in [26] under rather mild hypotheses. There is amajor difference with respect to the framework which has been assumed in those works,because it is typical into the study of planetary systems. There, the dynamics that istransverse to the elliptic tori is secular with respect to the revolution periods of the orbitslying on the elliptic tori; here, there is not any separation between the degrees of freedomaccording to whether they are fast or slow. An exhaustive analysis of the convergence ofour algorithm has been deferred to another work (see the reference reported in footnote 5 ),while here we prefer to focus on the semi-analytic applications to FPU chains with a smallnumber of nodes. This means that the formal procedure constructing the normal formsrelated to elliptic tori is explicitly implemented with the help of a specific computeralgebra software, specialized in the implementation of perturbation methods based oncanonical transformations expressed by Lie series. The results given by such a semi-analytic approach are compared with the purely numerical ones, which are produced byother technical tools that are used often in the field of Astronomy and are commonlyknown with the name of frequency analysis (FA) methods. FA is used here also in aversion that is extremely performing in the search of elliptic tori. We have devoted aspecial care in the comparisons between the semi-analytic results and the numerical ones,for what concerns the periodic orbits that are linearly stable in the transverse dynamics(describing the limit of small oscillations around the periodic orbits themselves). We haveconsidered both the α − model and the β one in the cases with N = 4 ,
8. As it was rathernatural to expect, we have obtained our best result when the perturbation is given bythe quartic terms of the potential energy (i.e., the β − model) with the lowest significantnumber of particles ( N = 4). In such a specific case, our semi-analytic algorithm hasbeen able to construct 1-D elliptic tori up to values of the energy that are about 93 %of their breakdown threshold, according to the numerical exploration performed by usingFA, which has been widely described in subsection 5.2. Such a level of performanceis comparable with that of the best computer-assisted proof (hereafter, CAP) for whatconcerns the existence of KAM tori for a Hamiltonian system, i.e., the so called forcedpendulum problem (see [9]). Let us also recall that, in some special models describedby symplectic maps, there are CAPs working up to values of the small parameter thatare extremely close to the numerical breakdown threshold for the wanted invariant tori(see [16]). In this context, it is natural to raise the following couple of interesting questions. lliptic tori in FPU non-linear chains . . . N going to infinity is beyond the scopes of the present work, it would be interesting toextend our approach to values of N larger than 4 or 8. In order to do that, in our semi-analytic method based on the construction of suitable normal forms we should includesome substantial improvements; let us sketch a few of them. It is well known that wheninitial conditions of “semi-sinusoidal” type and rather low values of the total energy areconsidered, the oscillation amplitudes of the high frequency normal modes are extremelysmall (see, e.g., [4]). Therefore, in our codes doing computer algebra manipulations weshould introduce a hierarchy in the truncation rules of the expansions (that are nowcompletely homogeneous), in such a way to approximate more carefully the dynamics ofthe low frequency normal modes with respect to those oscillating faster. Such a sourceof improvement is very technical, but it is essential in order to implement a normalform approach that is able to properly describe the behavior of the FPU chains alsofor a not so small number of nodes. Moreover, our semi-analytic method has shownbetter performances in the applications to β − models rather than to α − ones. In order tocompensate this imbalance, the introduction of action–angle variables well adapted to theToda lattice is probably unavoidable. Producing the preliminary expansions of the FPUHamiltonian in such coordinates is a task far from being trivial.In a recent work, we have improved very significantly some previous results aboutthe secular dynamics of the υ –Andromedæ planetary system, by constructing the normalforms for an elliptic torus and, therefore, for the orbit lying on a (KAM) invariant mani-fold. In order to achieve such a goal, the application of the same algorithm described andtested in the present paper has played a prominent role. Such an implementation has been6 C. Caracciolo, U. Locatelli so successful also because of the two following facts: in a purely secular approximationof a planetary system, there is not any more distinction between fast and slow dynamics(as in the framework considered here); moreover, the problem is described by a Hamilto-nian depending on the canonical variables as an even polynomial (i.e., the same situationholding true for the β − model). We think that the back transfer of our approach into thefield of Astronomy could provide interesting results. There are possible applications thatlook rather promising, for what concerns the partially unknown 3D architecture of theextrasolar planetary systems detected by the radial velocity method. This holds particu-larly true for the orbits of the exoplanets that are in a so called mean-motion resonance(see [48] as an example of a recent work about such a subject). Acknowledgments
This work was partially supported by the MIUR-PRIN project 20178CJA2B – “NewFrontiers of Celestial Mechanics: theory and Applications”, the “Beyond Borders” pro-gramme of the University of Rome Tor Vergata through the project ASTRID (CUPE84I19002250005) and by the “Progetto Giovani 2019” programme of the National Groupof Mathematical Physics (GNFM–INdAM) through the project “Low-dimensional Invari-ant Tori in FPU–like Lattices via Normal Forms”. The authors acknowledge also theMIUR Excellence Department Project awarded to the Department of Mathematics of theUniversity of Rome “Tor Vergata” (CUP E83C18000100006).
References [1] D. Bambusi, A. Ponno:
Resonance, Metastability and Blow up in FPU , in G. Galla-votti, man. ed., see [20], 191–205 (2008).[2] G. Benettin, A. Carati, L. Galgani, A. Giorgilli:
The Fermi–Pasta–Ulam Problemand the Metastability Perspective , in G. Gallavotti, man. ed., see [20], 151–189 (2008).[3] G. Benettin, A. Ponno:
Understanding the FPU state in FPU–like models , , doi:10.3934/mine.2021025 , 1–22 (2020).[4] L. Berchialla, L. Galgani, A. Giorgilli: Localization of energy in FPU chains , Discr.& Cont. Dyn. Syst., , 855–866 (2004).[5] C. Caracciolo, U. Locatelli: Computer-assisted estimates for Birkhoff normal forms ,J. of Comput. Dyn., , 425–460. (2020).[6] C. Caracciolo, U. Locatelli, M. Sansottera, M. Volpi: “Librational KAM tori in thesecular dynamics of the υ –Andromedæ planetary system”, preprint.[7] L. Casetti, M. Cerruti-Sola, M. Pettini, E.G.D. Cohen: The Fermi-Pasta-Ulamproblem revisited: Stochasticity thresholds in nonlinear Hamiltonian systems , Phys.Rev. E, , 6566–6574 (1997). lliptic tori in FPU non-linear chains . . . Local Behavior Near Quasi-Periodic Solu-tions of Conformally Symplectic Systems , J. of Dyn. & Diff. Eqs., , 821–841 (2013).[9] A. Celletti, A. Giorgilli, U. Locatelli: Improved estimates on the existence of invarianttori for Hamiltonian systems , Nonlinearity, , 397–412 (2000).[10] H. Christodoulidi, C. Efthymiopoulos, T. Bountis: Energy localization on q-tori,long-term stability, and the interpretation of Fermi-Pasta-Ulam recurrences , Phys.Rev. E, , 016210/1–16 (2010).[11] H. Christodoulidi, C. Efthymiopoulos: Low-dimensional q-Tori in FPU lattices: Dy-namics and localization properties , Physica D, , 92–113 (2013).[12] J. Couetdic, J. Laskar, A.C.M. Correia, M. Mayor, S. Udry:
Dynamical stabilityanalysis of the HD 202206 system and constraints to the planetary orbits , Astron. &Astroph., , A10 (2010).[13] A. Deprit:
Canonical transformations depending on a small parameter , CelestialMechanics, , 12–30 (1969).[14] L.H. Eliasson: Perturbations of stable invariant tori for Hamiltonian systems , Ann.Scuola Norm. Sup. Pisa, Cl. Sci., IV Serie, , 115–147 (1988).[15] E. Fermi, J. Pasta, S. Ulam: Studies of Nonlinear Problems , Los Alamos Report LA-1940 (1955); reprinted in both [20] and “Collected Papers of Enrico Fermi”, editedby E. Segr´e, University of Chicago, Vol. 2, p. 978, Chicago (1965).[16] J.-L. Figueras., A. Haro, A. Luque:
Rigorous computer-assisted application of KAMtheory: A modern approach , Found. Comput. Math., , n. 5, 1123–1193 (2017).[17] S. Flach, A. Ponno: The Fermi-Pasta-Ulam problem: Periodic orbits, normal formsand resonance overlap criteria , , 908–917 (2008).[18] H. Flaschka: The Toda lattice. II. Existence of integrals , Phys. Rev., , 1924–1925(1974).[19] F. Gabern, A. Jorba, U. Locatelli: On the construction of the Kolmogorov normalform for the Trojan asteroids , Nonlinearity, , n.4, 1705–1734 (2005).[20] G. Gallavotti, man. editor: The Fermi-Pasta-Ulam problem. A status report , LectureNotes in Physics, , Springer, Berlin (2008).[21] A. Giorgilli:
Quantitative methods in classical perturbation theory , proceedings ofthe Nato ASI school “From Newton to chaos: modern techniques for understandingand coping with chaos in N–body dynamical systems”, A.E. Roy e B.D. Steves eds.,Plenum Press, New York (1995).8
C. Caracciolo, U. Locatelli [22] A. Giorgilli:
Notes on exponential stability of Hamiltonian systems , in
DynamicalSystems, Part I . Pubbl. Cent. Ric. Mat. Ennio De Giorgi, Sc. Norm. Sup. Pisa,87–198 (2003).[23] A. Giorgilli, A. Delshams, E. Fontich, L. Galgani, C. Sim´o:
Effective stability fora Hamiltonian system near an elliptic equilibrium point, with an application to therestricted three body problem , J. Diff. Equations, , 167–198 (1989).[24] A. Giorgilli, U. Locatelli: Kolmogorov theorem and classical perturbation theory , J.of App. Math. and Phys. (ZAMP), , 220–261 (1997).[25] A. Giorgilli, U. Locatelli, M. Sansottera: Kolmogorov and Nekhoroshev theory for theproblem of three bodies , Cel. Mech. & Dyn. Astr., , 159–173 (2009).[26] A. Giorgilli, U. Locatelli, M. Sansottera:
On the convergence of an algorithm con-structing the normal form for lower dimensional elliptic tori in planetary systems ,Cel. Mech. & Dyn. Astr., , 397–424 (2014).[27] A. Giorgilli, U. Locatelli, M. Sansottera:
Secular dynamics of a planar model of theSun-Jupiter-Saturn-Uranus system; effective stability in the light of Kolmogorov andNekhoroshev theories , Regular and Chaotic Dynamics, , 54–77 (2017).[28] A. Giorgilli, M. Sansottera: Methods of algebraic manipulation in perturbation theory ,in “Chaos, Diffusion and Non-integrability in Hamiltonian Systems – Applications toAstronomy”, Proceedings of the Third La Plata International School on Astronomyand Geophysics, P.M. Cincotta, C.M. Giordano and C. Efthymiopoulos eds., Univer-sidad Nacional de La Plata and Asociaci´on Argentina de Astronom´ıa Publishers, LaPlata (2012).[29] W. Gr¨obner:
Die Lie-Reihen und Ihre Anwendungen , Springer Verlag, Berlin (1960);Italian transl.:
Le serie di Lie e le loro applicazioni , Cremonese, Roma (1973).[30] J.D. Hadjedemetriou:
Periodic orbits in gravitational systems , 43–79, in “ChaoticWorlds: from Order to Disorder in Gravitational N-Body Dynamical Systems”,B.A. Steves, A.J. Maciejewski, M. Hendry (eds)., Nato Science Series II, ,Springer Netherlands (2006).[31] J. Henrard:
The algorithm of the inverse for Lie transform , in “Recent advances inDynamical Astronomy”, B.D: Tapley and V. Szebehely, V. (eds.), Reidel, Dortrecht(1973).[32] G. Hori:
Theory of general perturbations with unspecified canonical variables , Publ.Astron. Soc. Japan, , 287–296 (1966).[33] A.N. Kolmogorov: Preservation of conditionally periodic movements with smallchange in the Hamilton function , Dokl. Akad. Nauk SSSR, , 527 (1954). Engl.transl. in: Los Alamos Scientific Laboratory translation LA-TR-71-67; reprinted in:Lecture Notes in Physics, . lliptic tori in FPU non-linear chains . . . Reeb’sTheorem and Periodic Orbits for a Rotating H´enon–Heiles Potential , J. Dyn. Diff.Eqs., https://doi.org/10.1007/s10884-019-09814-6 (2019).[35] J. Laskar:
Frequency Map analysis and quasi periodic decompositions , in Benest, D.,Froeschl´e, C. and Lega E. (eds.):
Hamiltonian systems and Fourier analysis , Taylorand Francis (2003).[36] J. Laskar, P. Robutel:
High order symplectic integrators for perturbed Hamiltoniansystems , Cel. Mech. & Dyn. Astr., , 39-62 (2001).[37] A. Luque, J. Villanueva: A KAM theorem without action-angle variables for ellipticlower-dimensional tori , Nonlinearity, , 1033-1080 (2011).[38] V.K. Melnikov: On some cases of conservation of almost periodic motions with asmall change of the Hamiltonian function , Dokl. Akad. Nauk SSSR, , 1245–1248(1965).[39] A. Morbidelli, A. Giorgilli:
Superexponential stability of KAM tori , J. Stat. Phys., , 1607–1617 (1995).[40] N.N. Nekhoroshev: Exponential estimates of the stability time of near–integrableHamiltonian systems.
Russ. Math. Surveys, , 1 (1977).[41] N.N. Nekhoroshev: Exponential estimates of the stability time of near–integrableHamiltonian systems, 2.
Trudy Sem. Petrovs., , 5 (1979).[42] B. Noyelles: Expression of Cassini’s third law for Callisto, and theory of its rotation ,Icarus, , 225–239 (2009).[43] B. Noyelles, N. Delsate, T. Carletti:
Equilibrium search algorithm of a perturbedquasi-integrable system , preprint (2012); https://arxiv.org/pdf/1101.2138.pdf [44] T. Penati, S. Flach:
Tail resonances of Fermi-Pasta-Ulam-breathers and their impacton the pathway to equipartition , Chaos, , 023102 (2007).[45] J. P¨oschel: On elliptic lower dimensional tori in Hamiltonian sytems , Math. Z., ,559–608 (1989).[46] G. Reeb:
Sur certaines propri´et´es topologiques des trajectoires des syst´emes dy-namiques , Acad. Roy. Belgique. Cl. Sci. M´em. Coll. in 8 ◦ . 27 article 9 (1952).[47] B. Rink: An Integrable Approximation for the Fermi–Pasta–Ulam Lattice , in G. Gal-lavotti, man. ed., see [20], 283–301 (2008).[48] M. Sansottera, A.-S. Libert: “Resonant Laplace-Lagrange theory for extrasolar sys-tems in mean-motion resonance”, Cel. Mech. & Dyn. Astr., , 38 (2019).0
C. Caracciolo, U. Locatelli [49] M. Sansottera, U. Locatelli, A. Giorgilli:
A semi-analytic algorithm for constructinglower dimensional elliptic tori in planetary systems , Cel. Mech. & Dyn. Astr., ,337–361 (2011).[50] C. Skokos, T. Bountis, C. Antonopoulos:
Geometrical properties of local dynamicsin Hamiltonian systems: The Generalized Alignment Index (GALI) method , PhysicaD, , 30–54 (2007).[51] M. Volpi, U. Locatelli, M. Sansottera:
A reverse KAM method to estimate unknownmutual inclinations in exoplanetary systems , Cel. Mech. & Dyn. Astr.,130